### Data sources

Three sources of data were used. Incidence of dementia and relative risk of death among demented subjects were estimated from the Paquid cohort. Paquid is a cohort of 3,777 subjects aged 65 years or over at inclusion representative from two French departments^{34}. Subjects randomly selected from the electoral rolls who agreed to participate were interviewed at home by trained neuropsychologists at baseline in 1989 and subsequently every two or three years over 27 years. Diagnosis of dementia was assessed using DSM IIIR criteria in a two-phase procedure including a screening by the neuropsychologist and a clinical examination at home by a neurologist. Vital status and exact date of death were collected all along the follow-up. Mortality in the Paquid cohort was shown to be very close from national mortality rate in France for the same period^{24}.

Prevalence and incidence of chronic use of BZD from age 65 were estimated from the Echantillon Généraliste des Bénéficiaires (EGB) using data from 2010 to 2015. EGB was established in 2005. This is a 1/97th nationwide representative random sample of the French health insurance beneficiaries whether they are receiving healthcare or not^{35}. In France, every citizen is covered by health insurance and EGB includes most insurance systems except a few systems for civil servants and student. Participants are followed until their death, emigration or insurance cessation while newly insured patients are added on a quarterly basis. This database includes age and gender and information on all out-hospital reimbursements of drugs for more than 600,000 subjects representative of the French population. In this work, relying on reimbursement data, chronic use of BZD was defined as 6 months of continuous use as in most published studies about the association between BZD use and dementia^{7,10,36}.

Finally the population size by gender at age 65 and the estimations (till 2015) and projections (from 2016) for the overall mortality rates by age and gender were provided by the French National Institute of Statistics (INSEE)^{37}, for each year from 1950 to 2070.

### “Illness-death” model

The method relies on a non-homogeneous Markov illness-death model depicted in Supplementary Fig S1 online where the three-states are “non-demented”, “demented” and “dead”. The transition intensities between the states that represent respectively the incidence of dementia ((alpha_{01})), the mortality of healthy subjects ((alpha_{02})) and the mortality of demented subjects ((alpha_{12})) were modeled by the following proportional intensity models:

$$ alpha_{01} left( {a,t,zleft( {a,t} right)} right) = alpha_{01}^{0} left( {a,t} right)theta_{01} left( a right)^{{zleft( {a,t} right)}} ; $$

(1)

$$ alpha_{02} left( {a,t,zleft( {a,t} right)} right) = alpha_{02}^{0} left( {a,t} right)theta_{02} left( a right)^{z(a,t)} ; $$

(2)

$$ alpha_{12} left( {a,t,zleft( {a,t} right)} right) = alpha_{12}^{0} left( {a,t} right)theta_{12} left( a right)^{{zleft( {a,t} right)}} = alpha_{02}^{0} left( {a,t} right)RRleft( a right)theta_{12} left( a right)^{{zleft( {a,t} right)}} ; $$

(3)

where (a) is the age, (t) is the calendar year and (zleft( {a,t} right)) is the time and age dependent BZD exposure: (zleft( {a,t} right)) equals 1 if the subject is exposed at time *t* and age *a* and 0 if not. Dementia incidence and mortality rates depend on both age and calendar time. The mortality among demented is proportional to the mortality of healthy subjects with a relative risk depending on age ((RR(a))) while (theta_{01} left( a right)), (theta_{02} left( a right)) and (theta_{12} left( a right)) are the relative risks associated with the exposure to BZD, respectively for dementia, mortality among healthy subjects and mortality among demented subjects.

### Estimations of the incidence and prevalence at age 65 of chronic use of BZD

Subjects from EGB were eligible for this study if they: (i) had been continuously affiliated to the main French Health Insurance scheme (so-called General scheme) at least from 2009/01/01 to 2010/02/01, (ii) had a complete follow-up from 2009 until death, date of general scheme exit, or until 31 December 2015, whichever came first, (iii) were aged 65 or over on 2010/01/01. Subjects were considered as chronic BZD users if they received reimbursement for BZD drugs for at least 6 months continuously. Prevalence of chronic use at age 65 was estimated from this sample. Incidence of chronic use was estimated by eliminating prevalent users of benzodiazepine at 2010/01/01.

### Modeling assumptions regarding exposure

A subject was considered as exposed as soon as he/she had BZD dispensing for at least 6 consecutive months ((z(a,t) = 1)). Then he/she remains exposed for any later time *t* and age *a* because we hypothesized that chronic BZD use increased the risk of dementia for the rest of the life. The relative risk of dementia for chronic BZD users was fixed at (theta_{01} = 1.6) according to Billioti de Gage et al.^{9}. We assumed that the relative risks of death for BZD exposure were identical among demented and non-demented subjects ((theta_{02} (a) = theta_{12} (a))) and constant by 5-years intervals. Based on Mathieu et al.’s results^{38}, we took (theta_{02} (a) = theta_{12} (a) = 2.45) for (65 le a < 70), (theta_{02} left( a right) = theta_{12} left( a right) = 1.69) for (70 le a < 75), (theta_{02} left( a right) = theta_{12} left( a right) = 1.3) for (75 le a < 80), (theta_{02} left( a right) = theta_{12} left( a right) = 1.1) for (80 le a < 85) and (theta_{02} left( a right) = theta_{12} left( a right) = 1) for (a ge 85). A sensitivity analysis was conducted assuming that BZD use did not increase mortality (theta_{02} left( a right) = theta_{12} left( a right) = 1.)

### Estimations of transition intensities in the illness-death models

The illness-death model stratified by gender was estimated from the Paquid cohort without accounting for BZD consumption. Indeed, in the Paquid cohort, we do not have the information to identify chronic BZD users as defined above: BZD consumption at each visit (every 2 or 3 years) was collected but without information on duration of use. Estimates were obtained by a semi-parametric penalized likelihood approach^{39}, using the R package SmoothHazard^{40}. This method avoids parametric assumptions on the risk functions and account for both the competing risk of death and the interval censoring of age at dementia between the visit of diagnosis and the previous one. The age-specific relative risk of death for demented subjects, (RRleft( a right)), was estimated as the ratio of the estimated mortalities among demented (alpha_{12} left( a right)) and non-demented subjects (alpha_{02} left( a right)) without any parametric assumptions.

As several studies have suggested a decreasing calendar time trend for dementia incidence^{25}, our main analysis assumed that dementia incidence decreased by 1% each year. This decreasing incidence may reflect the improvement of education or a better care of medical conditions associated with the risk of dementia (hypertension, depression,…). The estimated incidence from Paquid was considered as the incidence for the generations who were 75 years old in 1990 (the median generation in Paquid) and the reduction was 1% each year for the following generations. A sensitivity analysis was conducted assuming a constant age-specific incidence of dementia (alpha_{01}).

Then, using the gender-specific estimates of the incidence of dementia, (alpha_{01} left( {a,t} right)), and of the relative risk of death for demented subjects, (RRleft( a right)), as well as the nation-wide projections of age-and gender-specific mortality, we computed the age-, gender- and calendar year-specific projections of mortality for demented ((alpha_{12} (a,t))) and non-demented subjects ((alpha_{02} (a,t))) by solving a differential equation as described in Wanneveich et al.^{31}.

The age-, gender- and calendar year-specific incidence of dementia among non-BZD users (alpha_{01}^{0} left( {a,t} right)) was computed from the incidence of dementia in the overall population, the age-, gender- and calendar year-specific prevalence of exposure among non-demented subjects ((P_{Z/ND} left( {a,t} right))) and the relative risk (theta_{01} left( a right)) by solving the following equation:

$$ alpha_{01} left( {a,t} right) = P_{Z/ND} left( {a,t} right)alpha_{01}^{0} left( {a,t} right)theta_{01} left( a right) + (1 – P_{Z/ND} left( {a,t} right))alpha_{01}^{0} left( {a,t} right). $$

(4)

Similarly, mortality among exposed and unexposed demented subjects were computed by

$$ alpha_{12} left( {a,t} right) = P_{Z/D} left( {a,t} right)alpha_{12}^{0} left( {a,t} right)theta_{12} left( a right) + (1 – P_{Z/D} left( {a,t} right))alpha_{12}^{0} left( {a,t} right); $$

(5)

and mortality among exposed and unexposed non-demented subjects were computed by Eq. (4) replacing (alpha_{01}^{0} left( {a,t} right)) by (alpha_{02}^{0} left( {a,t} right)) and (theta_{01} left( a right)) by (theta_{02} left( a right)).

### Monte Carlo algorithm for projections in 2040

To provide projections for dementia for subjects older than 65 in 2040, the algorithm was the following:

For each cohort born in year (b) with (2010 – 105 le b le 2040 – 65) ((1935 le b le 1975)), we generated a sample of 10,000 subjects alive and non-demented aged 65 in (b + 65).

BZD use at age 65 was generated for each subject according to prevalence of BZD use at 65

For each subject and each year from (b + 66) to (b + 105) or death, we simulated successively the occurrence of BZD use, the death and the onset of dementia according to their respective annual incidences.

Then empirical estimates of the following epidemiological indicators were computed using data from the generated samples with formulas detailed in the Supplementary Online Resource 1: prevalence of exposure (current or past chronic BZD use), prevalence of dementia, life-expectancy without dementia, life-long probability of dementia, mean age at dementia, mean number of years spent with dementia.

A first run of the algorithm was used to estimate only the prevalence of exposure (past or current chronic BZD consumption) that was used to update the computation of the three transition intensities among unexposed subjects ((alpha_{kj}^{0} ; k = 0, 1; j = 1, 2)) by solving Eqs. (4) and (5). Using these updated values a second run was used to provide central estimations of epidemiological indicators.

Variances of the estimates of epidemiological indicators were computed with 100 runs of the algorithm using 100 different sets of values generated from the posterior distribution of all the parameters estimated in the preliminary steps: (alpha_{01} left( a right)), (RRleft( a right)), BZD chronic use prevalence and BZD chronic use incidence.

All computations were performed separately for women and men for 3 differents scenarios:

**Scenario 0**: With constant BZD consumption**Scenario 1**: Incidence after 65 of chronic BZD use divided by two from 2020**Scenario 2**: Incidence after 65 of chronic BZD use considered as null from 2020.

Details for scenarios 1 and 2 are given in Supplementary Online Resource 2. The Monte Carlo algorithm was implemented in the R package MCSPCD available on GitHub.