Mathematical Methods for Hippocampus Cell Turnover Carbon DNA dating in neuronal and nonneuronal cells from the hippocampus shows turnover. Turnover rate estimates were based on systematic exploration of mathematical models for different scenarios of lifetime cell turnover. average age = ( 1 − f ) t + f ( 1 − exp ( − r t ) ) r .

The best-fitted models suggest that only a subset of the neurons is renewed during lifetime and that a majority of the neurons do not renew significantly after birth. About one third (35%) of the neurons renew regularly. Within this fraction, the turnover rate varies among subjects, with a median turnover rate of 1.75%. The formula (based on scenario 2POP), for a subject age t, with turnover rate r and fraction renewing f, is Thus, for a 50-year-old subject (t = 50 years, f = 0.35, and r = 0.0175 per year) the average age of neurons would be 44 years. t 1 / 2 = ln 2 r .

Other scenarios suggest that there might be a decrease in turnover with age. For these scenarios, it is better to compute the half-life of a cell, rather than the average age. For scenarios with constant turnover rates (A, 2POP, XPOP), the half-life is ( g 1 g 1 + t 1 / 2 ) g 1 g 0 = 1 2 .

For scenarios E and 2POPE, the half-life is given by solving the equation ( g 1 g 1 + t 1 / 2 ) g 1 g 0 exp ( − d t 1 / 2 ) = 1 2 .

For scenario 2POPEd, the half-life is given by solving the equation Because scenarios E, 2POPE and 2POPEd predict a preferential renewal of young cells, the half-lives are smaller than for scenarios A, 2POP, and XPOP. Nonneuronal cells also form a heterogeneous population in terms of turnover rates. Best fits suggest that a slight majority of nonneuronal cells turnover at few percentages per year. Hippocampal neuronal (NeuN+) and nonneuronal (NeuN-) samples were analyzed. From a total number of 156 samples (NeuN+ n = 74; NeuN- n = 82), 55 NeuN+ and 65 NeuN- samples were analyzed. Samples were excluded when the amount of 14C was too small to yield precise measurements, or when missing data prevented turnover rate analysis. In a few cases (4 for neurons and 1 for nonneuronal cells), it was not possible to fit the measured 14C level using scenario 2POP (see details in Individual turnover analysis below). For the individual turnover analysis based on scenario 2POP, the number of samples used was n = 51 for neurons, and n = 64 for nonneuronal cells.

Mathematical Modeling of the Turnover Rates 14C data modeling was based on an approach developed previously ( Bergmann et al., 2009 Bergmann O.

Bhardwaj R.D.

Bernard S.

Zdunek S.

Barnabé-Heider F.

Walsh S.

Zupicich J.

Alkass K.

Buchholz B.A.

Druid H.

et al. Evidence for cardiomyocyte renewal in humans. C data modeling was based on an approach developed previously (). Briefly, a continuous birth-and-death model (von Foerster-McKendrick equations) represented by a linear transport equation, where cells are born at age 0 and age at speed 1 per year (cells get older by one year each year). A cell can die according to a mortality curve determined by the age of the cell and the age of the individual. Cells are born either to replace a dead cell, or according to a predefined birth rate that depends on the age of the individual. 14C levels between birth and death of the individual, it is possible to reproduce the 14C levels a sample would have if it followed the birth-and-death process. By fitting the models to the data, it is possible to infer how much cell renewal is needed to reproduce the observed 14C levels. Details on mathematical modeling can be found in ( Bergmann et al., 2009 Bergmann O.

Bhardwaj R.D.

Bernard S.

Zdunek S.

Barnabé-Heider F.

Walsh S.

Zupicich J.

Alkass K.

Buchholz B.A.

Druid H.

et al. Evidence for cardiomyocyte renewal in humans. Bernard et al., 2010 Bernard S.

Frisén J.

Spalding K.L. A mathematical model for the interpretation of nuclear bomb test derived 14C incorporation in biological systems. The solution of the birth-and-death process yields the distribution of ages of all the cells alive in an individual at any given age. By integrating along the atmosphericC levels between birth and death of the individual, it is possible to reproduce theC levels a sample would have if it followed the birth-and-death process. By fitting the models to the data, it is possible to infer how much cell renewal is needed to reproduce the observedC levels. Details on mathematical modeling can be found in (). To account for the food chain delay between integration of 14C into the biomass and the integration into replicating DNA, we introduced food lag of 1 year, that is, the atmospheric bomb curve was shifted to the right by one year, so that a cell born in 1987 would have a 14C level corresponding to 1986. A large proportion of neurons generated during development are produced before birth. To account for that, we considered that the beginning of development was 6 months before birth. All estimates were calculated with these parameters. 14C level. The equation reads ( PDE ) ∂ n ( t , a ) ∂ t ︸ change in person age + ∂ n ( t , a ) ∂ a ︸ change in cell age = − γ ( t , a ) n ( t , a ) ︸ loss of cells , ( Initial Condition ) n ( − 0.5 , a ) = δ ( a ) ︸ all cells aged 0 initially , ( Boundary Condition ) n ( t , 0 ) = b ( t ) ︸ cell birth .

The equation for the birth-and-death process is a linear partial differential equation (PDE) with an age structure that will be associated with the atmosphericC level. The equation reads The function n ( t , a ) is the cell density of age a in a person of age t. The term density is used because n is not a cell number; it has units in cells/age (exactly in N 0 cells/year, where N 0 is the initial number of cells). The PDE represents the fact that cells age at the same time as the person, and that they die at a rate γ ( t , a ) (units: per year). The death rate can depend either on the age of the cell itself (the variable a), or on the age of the person (the variable t). In the different scenarios described below, many types of death rates are presented. Cell birth is not included directly in the PDE but is supplied as a boundary condition, since we assume that cells are always born at age a = 0. This assumption relies on the biological assumption that new cells in the hippocampus (neurons or other cell types) come from the proliferation of progenitors and that the DNA is contemporary to their birth. The boundary condition describes how many cells are born in a person of age t. The boundary condition does not depend explicitly on cell age, because of the assumption that all cells are born at age a = 0. The form the birth rate b(t) can be relatively complex and will be described in the scenarios below. The initial condition describes what was the cell age distribution six months (0.5 years) before birth. We assumed that most of the neurons were formed before birth; this explains the choice of the timing for the initial condition. The right-hand-side of the initial condition state that all cells were initially of age a = 0. With an appropriate choice of the death rate γ and birth rate b, any renewal dynamics can be modeled, from the simplest one where cells are replaced at a constant rate, to more complex ones where cell number can vary in time and birth and death also vary in time. Although it is possible, we do not explore population feedback mechanisms for regulating the number of cells, i.e., the death and birth rates do not depend on the state of the system.

Description of the Scenarios Bergmann et al., 2009 Bergmann O.

Bhardwaj R.D.

Bernard S.

Zdunek S.

Barnabé-Heider F.

Walsh S.

Zupicich J.

Alkass K.

Buchholz B.A.

Druid H.

et al. Evidence for cardiomyocyte renewal in humans. All models used here are similar, in that they describe birth and death of cells. We refer to each flavor of the model as a scenario. A scenario represents a simple hypothesis on how cells are formed and replaced each year. (The letters are based on the notation introduced in). The simplest scenario, scenario A, assumes that all cells are formed at birth and that every cell is replaced when it dies. The probability for a cell to die during a time interval is constant. The rate at which cells are replaced is termed turnover rate. Scenario A A single parameter describes the annual turnover rate of all cells. Scenario A offers a minimalistic model against which more realistic scenarios can be compared. The single parameter is the death rate, which also corresponds, under steady state (constant cell number), to the turnover, or the amount of cells replaced each year. For that reason, we present the results for scenario A even though it might not be the best to describe the data. Death rate, γ ( t , a ) = r constant; birth rate, b ( t ) = r . Scenario B The two parameters, the birth and death rates, are independent. In particular, cell number increases or decreases with time when birth and death rates are different. Death rate, γ ( t , a ) = r constant; birth rate, b ( t ) = b constant. Scenario C The turnover rate is decreasing exponentially with the age of the subject: r 0 exp ( − r 1 t ) . The two parameters are the initial turnover rate r 0 and the turnover decrease exponent r 1 . All dying cells are replaced. Death rate, γ ( t , a ) = γ ( t ) = r 0 exp ( − r 1 t ) ; birth rate, b ( t ) = γ ( t ) . Scenario E g 0 g 1 / ( g 1 + a ) . The two parameters are the initial death rate and the half-death constant, the age at which death rate is halved. All dying cells are replaced. The birth rate must be determined by solving a Volterra equation of type II, also called renewal equation. Death rate, γ ( t , a ) = γ ( a ) = g 0 g 1 / ( g 1 + a ) ; birth rate is obtained by solving the Volterra equation: b ( t ) ︸ number of new cells needed = γ ( t ) exp ( − ∫ 0 t γ ( s ) d s ) ︸ cells from birth dying today + ∫ 0 t b ( t − a ) γ ( a ) exp ( − ∫ 0 a γ ( s ) d s ) ︸ cells born a years ago dying today d a .

The death rate of the cells decreases with their age a according to the function. The two parameters are the initial death rate and the half-death constant, the age at which death rate is halved. All dying cells are replaced. The birth rate must be determined by solving a Volterra equation of type II, also called renewal equation. Death rate,; birth rate is obtained by solving the Volterra equation: Scenario 2POP The population consists of two populations, one renewing and one not renewing. The renewing fraction follows scenario A. The two parameters are the fraction of renewing cells and the turnover rate within this fraction. Death rate, γ ( t , a ) = r constant; birth rate, b ( t ) = r ; fraction of renewing subset, f constant. Scenario XPOP r ( x ) of population x, with x between 0 and 1. The variable x represents the ‘turnover class’, with x = 0 the class with the highest turnover and x = 1 the class with the lowest turnover rate. The turnover rate function is r ( x ) = r 0 x h k h + x h .

The population consists of a continuum of renewing populations, each with a different turnover rate. The three parameters determine the turnover rateof population x, with x between 0 and 1. The variable x represents the ‘turnover class’, with x = 0 the class with the highest turnover and x = 1 the class with the lowest turnover rate. The turnover rate function is The three parameters are, r 0 the highest turnover rate (corresponding to class x = 0), k the fraction of cells with turnover rates higher than r 0 / 2 and h a measure of heterogeneity (a large value of h means that there is a clear separation between a renewing and a nonrenewing fraction). ( PDE ) ∂ n ( t , a , x ) ∂ t ︸ change in person age + ∂ n ( t , a , x ) ∂ a ︸ change in cell age = − r ( x ) n ( t , a , x ) ︸ loss of cells turnover class x , ( Initial Condition ) n ( − 0.5 , a , x ) = { δ ( a ) if x ∈ [ 0 , 1 ] , 0 otherwise , ︸ all cells aged 0 initially ( Boundary Condition ) n ( t , 0 , x ) = r ( x ) ︸ cell birth .

For this scenario, the variable x acts as a population structure variable, like the age a does. Therefore, formally we need to introduce a new PDE, along with an initial condition and a boundary condition to describe the scenario. However, because cells always stay in the same turnover class, the PDE associated to scenario XPOP is similar to the PDE for the other scenarios. Here, the initial cell density is uniform over the interval [0,1], i.e., each turnover class contains the same ‘number’ of cells (it is actually a density). Because all cells are replaced and cells cannot move from one class to another, the cell density will stay uniform over [0,1] at any age of the person. The parameter h control which proportion of the cells have a high-turnover rate by setting a threshold on x over which the turnover rate function r ( x ) falls off rapidly. In practice, to solve numerically the PDE for this scenario, it is necessary to discretize the variable t, a and x, so x can be considered as a discrete class. We have used 17 turnover classes linearly spaced between 0 and 1: x = i / 16 , i = 0 , … , 16 . This way, scenario XPOP can be handled by much of the same numerical methods developed for the other scenarios. Scenario 2POPE γ ( t , a ) = γ ( a ) = g 0 g 1 ( g 1 + a ) .

This is a combination of scenario 2POP with scenario E. The three parameters are the fraction of renewing cells, the initial death rate and the half-death constant. Death rate of the renewing subset: Birth rate of the renewing subset is given by the Volterra equation introduced in scenario E. Fraction of renewing subset, f constant. Scenario 2POPEd γ ( t , a ) = γ ( a ) = g 0 g 1 ( g 1 + a ) .

Like 2POPE, but nonrenewing cells die at a constant rate. The four parameters are the fraction of renewing cells (at birth, since the fraction is changing), the death rate within the nonrenewing fraction, the initial death rate within the renewing fraction and the half-death constant. This is the most complete scenario, used for Figure 6 in the main text. Death rate of the renewing subset: Birth rate of the renewing subset is given by the Volterra equation introduced in scenario E. Death rate of the nonrenewing subset, d constant; fraction of renewing subset initially, f constant. Because the number of nonrenewing cell is declining with age of the person, the actual fraction of renewing cells is increasing. Other Scenarios Bergmann et al. (2009) Bergmann O.

Bhardwaj R.D.

Bernard S.

Zdunek S.

Barnabé-Heider F.

Walsh S.

Zupicich J.

Alkass K.

Buchholz B.A.

Druid H.

et al. Evidence for cardiomyocyte renewal in humans. Other scenarios, either from, or new scenarios were also explored. Different death rates functions were implemented and tested for scenarios C and E. scenario D, were the death rate increases with the age of the cells, was also tested, as well as scenarios 2POPB and 2POPC, which were built the same way as scenario 2POPE. More complex scenarios (F, G, H) with varying cell numbers and complex birth-and-death dynamics were also tested. All these scenarios either gave inconsistent results (cell number too high or too low) or gave the same results as a simpler scenario. For instance scenario D consistently predicted no increase in cell death rate, and was therefore not more informative than scenario A.

Calculating the 14C Level Associated to a Scenario n ( t , a ) can be numerically computed by direct integration of the PDE, given the initial condition and boundary condition. (All the methods apply also to scenario XPOP). The average 14C level C ˜ of a DNA sample of a person born at calendar year t b and collected at calendar year t d (here this is the date of death of the person), is C ˜ = ∫ 0 t = t d − t b K l a g ( t d − a ) n ( t , a ) ︷ contribution of cells aged a to 14 C d a N ( t ) ︸ Total number of cells .

The cell densitycan be numerically computed by direct integration of the PDE, given the initial condition and boundary condition. (All the methods apply also to scenario XPOP). The averageC levelof a DNA sample of a person born at calendar year tand collected at calendar year t(here this is the date of death of the person), is K l a g is the food-lag atmospheric 14C level curve. When evaluated at calendar year y, K l a g ( y ) represents the actual 14C concentration that will integrate newly synthesized DNA, and will therefore correspond to an average of past atmospheric concentration, to account for the food supply chain, from photosynthesis to the plate. Here we have considered two different forms of K l a g functions. The simplest form, and the one that was used for the numerical computation is K l a g ( y ) = K ( y − t l a g ) , i.e., delaying the atmospheric curve by a certain amount of time t l a g . After exploring different values of t l a g , we settled to t l a g = 1 year, given that any values between 0 and 2 years did not change numerical results significantly. To test whether a smoother lag representing the fact that food comes from different sources with different delays would be more appropriate, we also tested convoluting the atmospheric 14C levels with a delay kernel K l a g ( y ) = ∫ 0 τ m a x K ( y − τ ) g ( τ ; k , q ) d τ

The functionis the food-lag atmosphericC level curve. When evaluated at calendar year y,represents the actualC concentration that will integrate newly synthesized DNA, and will therefore correspond to an average of past atmospheric concentration, to account for the food supply chain, from photosynthesis to the plate. Here we have considered two different forms offunctions. The simplest form, and the one that was used for the numerical computation is, i.e., delaying the atmospheric curve by a certain amount of time. After exploring different values of, we settled toyear, given that any values between 0 and 2 years did not change numerical results significantly. To test whether a smoother lag representing the fact that food comes from different sources with different delays would be more appropriate, we also tested convoluting the atmosphericC levels with a delay kernel g ( τ ; k , q ) = k q Γ ( q ) τ q − 1 exp ( − k τ ) .

We choose as a delay kernel the gamma function, which is used to describe the time needed to go through a fixed number of steps, each characterized by an exponential waiting time. The gamma distribution has two parameters, the speed k of each step and the number of steps q, For integers, Γ ( q ) = ( q − 1 ) ! is the factorial function. Because of the additional integral that had to be computed for the delay kernel, numerical simulations were significantly slowed down, so we tested only a few scenarios. For delay kernel restricted on less than three years, results were similar to that of the delayed 14C levels.

Numerical Simulations of the Scenarios All simulations were performed with MATLAB (version R2012b), with parallel codes developed specifically for the task. Briefly, all equations were integrated using appropriate numerical integral algorithms (nonadaptive Simpson quadrature and adaptive Simpson quadrature). The atmospheric 14C level curve was sampled at midpoint each year (1993.5, 1994.5, and so on) and linearly interpolated to convert it to a continuous function for use in the numerical integral functions. Simulations were performed on an iMac with Intel i5-2500s quad-core CPU at 2.70GHz. An entire set of simulations takes around 16 days to run.

Individual Turnover Analysis Individual turnover rates can be estimated in individual subjects for scenario A. Other scenarios have too many parameters to obtain a unique estimate. For scenario 2POP, we used a fraction of renewing cells of f = 0.35 for the neuronal cells, and f = 0.55 for the nonneuronal cells. The other scenarios were not fitted to individual samples because there were too many parameters to fit. Scenario A: Individual Turnover Rates Neuronal cells (NeuN+, n = 55) had mean and median individual turnover rates (per year) of 0.0077 and 0.0048, respectively. Nonneuronal cells (NeuN-, n = 65) had mean and median individual turnover rates of 0.0174 and 0.0118, respectively. The p value for neuronal versus nonneuronal cells (Wilcoxon rank sum test) was 5 × 10−8. Scenario 2POP: Individual Turnover Rates (4 neuronal cell outliers and 1 nonneuronal cell outliers were removed). Neuronal cells (NeuN+, n = 51) had a median individual turnover rate (per year) of 0.0175. Nonneuronal cells (NeuN-, n = 64) had a median individual turnover rate of 0.0350. The p value for neuronal versus nonneuronal cells (Wilcoxon rank sum test) was 0.005. scenario 2POP assumes that only a fraction of the cells are renewing. The size of the renewing fraction was set to 0.35 for neurons and 0.55 for nonneuronal cells, based on globally estimated values (see Global turnover rate analysis for details on global analysis). In a few cases (4 for neurons and 1 for nonneuronal cells), it was not possible to fit the measured 14C level using these renewal fractions, indicating that in those cases, the renewing fraction would be larger. However, for the comparison of turnover rates, is was important to keep the renewing fractions equal for all subjects, so the sample that could not be fitted were excluded. Pairwise Analysis between Neuronal and Nonneuronal Cells Paired neuronal cell analysis (NeuN+, n = 40) had a median individual turnover rate per year of 0.0047 (scenario A) and 0.0150 (2POP). Paired nonneuronal cell analysis (NeuN-, n = 40) had a median individual turnover rate of 0.0120 (scenario A) and 0.0350 (2POP). The p value for paired neuronal versus nonneuronal cells (Wilcoxon signed-rank test) was 2e-8 for scenario A and 0.054 for scenario 2POP.

Global Turnover Rate Analysis To fit the birth-and-death models to the data, we have assumed that the calculated 14C level C could be compared directly to the measured 14C levels in a DNA sample. This entails two biological assumptions we thought were realistic: i) most of neurons and other hippocampal cells are postmitotic cells, and their DNA is stable after the last division, ii) new cells come from a small but actively proliferating pool of progenitor, so that DNA synthesis corresponds to the birth of the cell. Old postmitotic cells migrating from a slowly depleting reservoir was not considered. With these assumptions in mind, we used many approaches to try to find the best parameters values for each scenario. This means finding within a parameter space which parameters values will optimize a predefined goodness-of-fit function. The dimension of the parameter space in each scenario is equal to the number of parameters to estimate. This varies from 1 for scenario A to 4 for scenario 2POPEd. The higher the dimension, the harder it is to find an optimal parameter set. We explain below how the problem of selecting best fits and comparing scenarios was tackled. ( C ˜ ) and the measured ones (C). The SSE is S S E = ∑ i = 1 n ( C ˜ i − C i ) 2 w i 2 .

To estimate the goodness-of-fit of a parameter set, we looked at a weighted sum of squares of the errors (SSE) between the calculated (predicted) carbon levelsand the measured ones (C). The SSE is The weight w i corresponds to the confidence we had in a 14C sample, based on the amount of carbon measured. Small samples (below 10 μg) were given a weight of 1, while samples of mass over 10 μg where assigned a weight of 2. For any given scenario, we found, using a direct search method for minimizing a multivariate function (fminsearch in MATLAB), a minimum of the SSE. The SSE was summed up over all neuronal cell samples and over all nonneuronal cell samples (in the SSE formula, n is the number of samples). The index i refers to a single sample. For the minimum SSE search, an initial guess for the parameter set must be provided. The main problem with minimizing the SSE is that in general, a minimum is found locally to that initial guess, and there is no guarantee that there is no other parameter set that will have a lower SSE in some other part of the parameter space. The only practical way around the local minimum problem is to explore the parameter space in a systematic way. First we have used different initial guesses to make sure that fminsearch would always find the same minimum, so that we could be confident that the minimum is not only local but also global. Then we used a Markov Chain Monte Carlo method to explore the space in a systematic way. This allowed us not only to confirm that the parameter sets that minimized the SSE were most probably the best fits (that is, with a SSE globally minimal), but also, in certain case, to estimate the uncertainty on the parameter values. The algorithm used is the Metropolis-Hastings algorithm. Briefly, successive parameter sets are tried and compared to the current parameter set. If the new parameter set is better (larger likelihood, or smaller SSE), it is accepted and replaces the current parameter set. If it is not better, it replaces the current parameter set with a probability related to the likelihood ratio between the two sets. This process is iterated a fixed number of times. Accepted and rejected parameter sets, along with their likelihood, give a global picture of the parameter space, and allows the determination of parameter ranges that are compatible with the data. At each iteration, a new parameter set was determined randomly within an interval around the current parameter set. The size of the interval was adapted so that the probability of accepting the new parameter set is around 25%. New parameter sets falling outside a pre-set interval were automatically rejected (preset intervals are indicated for each scenario in the corresponding figure legend). For the analysis, the first 10,000 iterations were discarded, to remove the transients (the burn-in). Bergmann et al., 2009 Bergmann O.

Bhardwaj R.D.

Bernard S.

Zdunek S.

Barnabé-Heider F.

Walsh S.

Zupicich J.

Alkass K.

Buchholz B.A.

Druid H.

et al. Evidence for cardiomyocyte renewal in humans. A I C = p + log 2 S S E . Although the AIC cannot be rigorously applied to the types of models presented here, it provides a rule of thumb to select the best scenario. With this criterion, scenario E (2 parameters) must have a SSE of less than half scenario A (1 parameter) to be considered better. We expect that more complex scenarios fit the data better, because they have more free parameters to adjust. It does not mean that they are better than a simple scenario. We had already used before a simple criterion, the AIC, to compare how well each scenario fitted the data (). The AIC states that each additional parameter in the model should make the SSE decrease by a factor 2, otherwise the additional parameter is not useful:Although the AIC cannot be rigorously applied to the types of models presented here, it provides a rule of thumb to select the best scenario. With this criterion, scenario E (2 parameters) must have a SSE of less than half scenario A (1 parameter) to be considered better.

Toward a “Preferential Loss of Younger Neurons” Scenario? When evaluated by their goodness-of-fit, the best scenarios all predict that young neurons die younger. This ‘preferential loss of younger neurons’ can be achieved through two distinct mechanisms. First, the best overall scenario (scenario 2POP) suggests that a subset of hippocampal neurons renew at a high rate, while another subset does not renew at a significant rate. This leads to a bimodal population consisting of young, renewing neurons, and older, nonrenewing neurons. Second, there might be a continuum in the renewing rates with respect to neuron ages. Younger neurons would be more likely to be replaced than older neurons, as in scenario E. Available data suggest that scenario 2POP explains the data best. This is mainly due to the contribution of the subjects born during the years 1955-1965, around the sharp increase and decrease of the atmospheric 14C levels. 14C concentrations in neurons samples are quite low during that time, indicating that a significant proportion of the cells are not renewed. Without these subjects, it was not possible to distinguish between scenario 2POP and a homogeneous decreasing turnover such as scenario C or E.

Markov Chain Monte Carlo Simulations for Nonneuronal Cells Six scenarios were tested using MCMC simulations: scenario A, scenario E ( Figure S2 A), scenario 2POP ( Figure S2 B), scenario XPOP ( Figure S2 C), scenario 2POPE ( Figure S2 D), scenario 2POPEd ( Figure S2 E). Scenario A: Nonneuronal Cells Markov Chain Monte Carlo simulations for scenario A confirmed there is only one best fit, with a turnover rate of 1.426% per year, consistent with the nonlinear least-square fit. The simulation explored turnover rate values between 0 and 1 per year, with 1 × 105 iterations. Scenario E: Nonneuronal Cells 5 ( Markov Chain Monte Carlo parameter estimation for scenario E. The initial death rate was bounded between 0 and 1 per year and the half-death rate constant between 0 and 80 years. The best fit (3.727% per year, 14.35 years) is consistent with the nonlinear least-square fit. Number of Monte Carlo iterations: 1 × 10 Figure S2 A). Scenario 2POP: Nonneuronal Cells 5 ( Markov Chain Monte Carlo parameter estimation for scenario 2POP. The best fit (0.5596 fraction renewing cells, 0.04923 per year) is consistent with the nonlinear least-square fit. The fraction of renewing cells was bounded between 0 and 1 and the turnover rate between 0 and 0.4 per year. Number of iterations: 1 × 10 Figure S2 B). Scenario XPOP: Nonneuronal Cells 0 , the maximal turnover rate; g 1 the fraction of cells having a renewal rate of at least g 0 / 2 , and the hill coefficient, which controls the steepness of the slope around g 1 . According to the MCMC simulations, the acceptable range for g 0 is g 0 > 0.05 per year, g 1 between 0.1 and 0.5 (interpreted as the fraction with high turnover) and a hill coefficient above 3. The results of the MCMC simulations are consistent with scenario 2POP. Number of iteration: 3 × 105 ( Markov Chain Monte Carlo parameter estimation for scenario XPOP. A three-dimensional space parameter was sampled. The three parameters are: g, the maximal turnover rate; gthe fraction of cells having a renewal rate of at least, and the hill coefficient, which controls the steepness of the slope around g. According to the MCMC simulations, the acceptable range for gis g> 0.05 per year, gbetween 0.1 and 0.5 (interpreted as the fraction with high turnover) and a hill coefficient above 3. The results of the MCMC simulations are consistent with scenario 2POP. Number of iteration: 3 × 10 Figure S2 C). 0 . Values of g 0 were divided into 100 bins ranging from 0.0 to 1.0 per year, and the hill function corresponding to the parameter set with the highest likelihood was plotted ( Best fits with the MCMC simulation were analyzed. All parameter values with a likelihood L > 0.55 were selected and sorted according to the parameter g. Values of gwere divided into 100 bins ranging from 0.0 to 1.0 per year, and the hill function corresponding to the parameter set with the highest likelihood was plotted ( Figure 4 B). Scenario 2POPE: Nonneuronal Cells 5 ( Markov Chain Monte Carlo simulation for the combined scenario 2POPE. A three-dimensional space parameter has been sampled. The three parameters are the fraction of renewing cells (range [0, 1]), the initial death rate (range [0, 0.2] per year) and the half-death rate constant (range [0, 80] years). The fraction renewing is most likely to be between 0.5 and 0.6. The initial death rate is likely to around 0.05 per year. The half-death rate constant (indicating at what age the cell have divided their death rate by two) is most likely around 20 years or very large (this would indicate no decrease in turnover rate). Number of iterations: 3 × 10 Figure S2 D). Scenario 2POPEd: Nonneuronal Cells 5 ( Markov Chain Monte Carlo simulation with scenario 2POPEd. A four-dimensional space parameter has been sampled. Nonrenewing cells die at a slow rate, and renewing cells replace young cells preferentially. The four parameters are the fraction of renewing cell (range [0, 1]), the death rate of nonrenewing cells (range [0, 0.01] per year), the initial death rate (range [0, 0.2] per year) and the half-death rate constant (range [0, 80] years). Number of iterations: 3 × 10 Figure S2 E).