 Research
 Open Access
 Published:
Multitype BellmanHarris branching model provides biological predictors of early stages of adult hippocampal neurogenesis
BMC Systems Biology volume 11, Article number: 90 (2017)
Abstract
Background
Adult hippocampal neurogenesis, the process of formation of new neurons, occurs throughout life in the hippocampus. New neurons have been associated with learning and memory as well as mood control, and impaired neurogenesis has been linked to depression, schizophrenia, autism and cognitive decline during aging. Thus, understanding the biological properties of adult neurogenesis has important implications for human health. Computational models of neurogenesis have attempted to derive biologically relevant knowledge, hard to achieve using experimentation. However, the majority of the computational studies have predominantly focused on the late stages of neurogenesis, when newborn neurons integrate into hippocampal circuitry. Little is known about the early stages that regulate proliferation, differentiation, and survival of neural stem cells and their immediate progeny.
Results
Here, based on the branching process theory and biological evidence, we developed a computational model that represents the early stage hippocampal neurogenic cascade and allows prediction of the overall efficiency of neurogenesis in both normal and diseased conditions. Using this stochastic model with a simulation program, we derived the equilibrium distribution of cell population and simulated the progression of the neurogenic cascade. Using BrdU pulseandchase experiment to label proliferating cells and their progeny in vivo, we quantified labeled newborn cells and fit the model on the experimental data. Our simulation results reveal unknown but meaningful biological parameters, among which the most critical ones are apoptotic rates at different stages of the neurogenic cascade: apoptotic rates reach maximum at the stage of neuroblasts; the probability of neuroprogenitor cell renewal is low; the neuroblast stage has the highest temporal variance within the cell types of the neurogenic cascade, while the apoptotic stage is short.
Conclusion
At a practical level, the stochastic model and simulation framework we developed will enable us to predict overall efficiency of hippocampal neurogenesis in both normal and diseased conditions. It can also generate predictions of the behavior of the neurogenic system under perturbations such as increase or decrease of apoptosis due to disease or treatment.
Background
Adult neurogenesis generates new neurons throughout life in two distinct regions of the mammalian brain: the subventricular zone, involved in olfactory processes, and the subgranular zone (SGZ) of the dentate gyrus [1–4], where new neurons have been associated with learning and memory as well as mood control [5–7]. The addition of new neurons is not merely static or restorative; it constitutes an adaptive response to the animal’s environment and/or its internal state. For example, hippocampal neurogenesis can be regulated both positively and negatively by external stimuli, such as learning [8], exercise [9], environment [10] and stress [11], as well as pathological states such as epilepsy [12–16], drug addiction [17–19], depression [20–22] and schizophrenia [23, 24]. Thus, adult neurogenesis represents another means, apart from molecular, synaptic, or morphological changes of an individual cell, to alter the functional circuitry depending on the demand. However, despite a significant functional relevance of this form of wholecell plasticity, little is known about the processes that regulate it.
During physiological conditions, adult neurogenesis maintains a steadystate. At any given moment, neural stem and progenitor cells (NPCs) may undergo one of three possible fates – they proliferate, producing more of identical daughter cells; they differentiate, giving rise to immature neurons called neuroblasts; or they die [25–28] (Fig. 1). It is believed that the basal rate of neurogenesis is genetically determined [29], but the mechanisms that govern it under various physiological and pathological stimuli are poorly understood. Most research on the neurogenic regulatory mechanisms has centered on the factors that regulate NPC proliferation and integration of newborn neurons into the dentate circuitry during the late stages of neurogenic cascade [30, 31]. However, early stages of neurogenesis are very complex, as mechanisms that determine cell proliferation, differentiation, and death are active at the same time. Further, the influence of the newborn cell death on adult neurogenesis is not known, even though it has been established that the majority of newborn cells in the adult dentate SGZ die during the first week of life, presumably undergoing apoptosis [32]. Newborn cell apoptosis may also be important for spatial learning [33], a hippocampaldependent task suggested to require neurogenesis [34]. Thus, understanding the early stages of neurogenesis is critical if we are to manipulate this process to enhance the number of viable newborn neurons as treatment modalities.
As experimental studies of such a complex system require years of work, several groups have used computational tools to aid discovery and guide biological experimentation. Most of the existing models, however, focus on the late stages of neurogenesis and aim to understand the effects of new neuron incorporation into the dentate gyrus. These models have shown that new granule cells participate in pattern separation, avoiding interference between memories while older ones are not greatly disturbed [35–37]. However, existing models have not addressed all the processes that occur throughout the neurogenic cascade, specifically all three possible fates of cells that are part of the neurogenic cascade  proliferation, differentiation and cell death. Hence, here we propose a comprehensive computational model of all stages of neurogenic cascade, including transition, proliferation, differentiation and survival of newborn cells from the NPC stage to the stage of a matured neuron. To model early stages of neurogenesis, we use the theory of branching processes [38]. In a hierarchical system such as that formed by proliferating and differentiating ANPs, the theory allows to formulate explicit analytic solutions in the terms of multiple but finiteorder convolutions of distributions of transit times through different phases of cell cycle. This latter feature makes modeling particularly transparent, and allows avoiding purely numerical simulations. Branching process theory can be traced to the social scientists in the 19th century studying the extinction of family lines. From that time on, a large number of biological problems have been modeled by branching processes, particularly in the analysis of evolutionary cell population and population genetics. For example, during the evolution of a population of some reproducing particles, each particle lives for a lifetime, independently of the others, and produces a random number of new offspring. If each particle lives for a constant unit of time and produces progeny upon death, then the process is called a GaltonWatson branching process. If each particle has an exponentially distributed lifetime independent of the offspring distribution, then the process is called a Markov continuoustime branching process. If the lifetime of each particle is a random variable with an arbitrary distribution, independent of lifetimes of the offspring, then this process is named an agedependent (BellmanHarris) branching process [38]. Here, using the Multitype BellmanHarris branching model, we provide for the first time estimates of the early stages of the neurogenic cascade, focusing on the apoptosis and transit times of cells, from birth to incorporation into the hippocampal circuitry.
Methods
Experimental methods
Animals
Wildtype (C57BL/6) or transgenic NestinCFPnuc mice, in which CFP is expressed in the nuclei of both neural stem cells (NSCs) and ANPs [21], were used. All mouse studies were approved by the Baylor College of Medicine Institutional Animal Care and Use Committee and performed in accordance with institutional and federal guidelines. Unless otherwise stated, animals were 1 month old.
Cell labeling with Bromodeoxyuridine (BrdU)
When studying neurogenesis, the most accepted method to estimate the net effect of the neurogenic cascade is to use BrdU, which labels cells in S phase of the cell cycle, to trace proliferation and differentiation (Fig. 2). BrdU is injected and during the circulating 15 min time, it gets incorporated into a proliferating DNA. Over the course of the neurogenesis, the BrdU can be traced in cells that are the lineage of the initial proliferating cell. BrdU labeling can be done as a single or cumulative injection paradigm (Fig. 3). In single labeling experiments, animals were injected with 250 mg/kg BrdU at t=0. In cumulative labeling experiment, performed to obtain the highest yield of the apoptotic cells, animals were injected with 150 mg/kg BrdU every 3 h in the first 24 h, totally 9 injections including the one at t=0. Animals were then sacrificed at different time points, when the total number of BrdU cell as well as the percentage of cells in each stage of the cascade was quantified.
Histology
Mice were transcardially perfused with phosphate buffer saline (PBS) followed by 4% paraformaldehyde (PFA). The brains were dissected out submerged into 4% PFA for 4 h at room temperature (RT) and sectioned using a vibratome. For immunofluorescence, freefloating sections were immunolabeled according to conventional procedures. The brains were dissected out, and then transferred to a cryoprotectant solution (30% sucrose, 30% ethylene glycol in PBS) and stored at 20 °C. Once brains were well cryoprotected, six series of 50 μm lateral sections were collected using a vibratome. A full series of freefloating sections were immunostained following conventional procedures [39]. Briefly, all washes and incubations were done in PBS containing 3% bovine serum albumin (BSA; SigmaAldrich) and 0.3% Triton X100 (SigmaAldrich). An antigen retrieval step (2 N HCl, 15 min, 37 °C) for BrdU detection was performed, followed by extensive washes with borate buffer (0.1 M). Sections were preincubated in PBS containing 3% BSA, 5% normal goat serum (Vector Labs) and 0.3% Triton X100 for 2 h at RT, followed by overnight incubation with primary antibodies (see below) at 4 °C. After extensive washing, sections were incubated with the appropriate secondary antibody conjugated with Alexa 488 (Molecular Probes), Rhodamine RedX and Cyanine 5 (Jackson Immunoresearch) together with DAPI (5 μg/mL, SigmaAldrich) for 2 h at RT. They were then washed, and mounted on slides with Fluorescent Mounting Medium (Dako). The following primary antibodies were used: BrdU (1:400, Accurate) to detect proliferating cells in S phase; DCX (Cell Signaling, 1:200) to detect neuroblasts/immature neurons; GFAP (1:1,000, SigmaAldrich) to detect primary neural stem cells (NSCs) and distinguish them from the ANPs (GFAP); NeuN (1:1,000, Chemicon) to detect mature neurons of the dentate gyrus, granule cells.
Confocal Microscopy
Sections were imaged with a Zeiss LSM or a Leica SP5 confocal microscope. The number of apoptotic cells and/or BrdU+ cells per zstack was estimated via the optical dissector method [39]. Blind analysis was performed with AxioVision 4.5 (Zeiss) or LAS AF Lite (Leica). Twothree 20 μm zstacks (consisting of 30 optical slices of 0.8 μm thickness) were obtained from every section. The number of cells was evaluated as a function of the volume of the SGZ, defined as the SGZ length in the image multiplied by an optical thickness of 20 μm and a height of 20 μm (which we defined in these experiments as a layer of cells expanding 5 μm into the hilus and 15 μm into the granular layer), then extrapolated to the volume spanned by the SGZ in the hippocampus.
Computational methods
Cell transit times
Cell transit time is defined as the duration time of a cell spent in the phase or stage before it transits into the next phase or stage. Instead of the commonly used exponential distributions, we use shifted gamma distributions of the transit times [40] to model the cell transit times through phases G1, S and G2M, as well as the lifetime of cells. Advantage of these is the ability to independently specify the minimum transit times, mean transit times, and variances of transit times. Among other, this allows avoiding occurrence of cells that live for an indefinitely short period of time with certain probability. Such distribution has three parameters, (k,s,v), with its probability density function given as f(xk,s,v)=(x−v)^{k−1} e ^{−(x−v)/s}/(Γ(k)s ^{k}), where k is the shape parameter, s is the scale parameter and v is the shift value (minimum duration), and Γ(k) is Euler gamma function [40]. Cells in different cell cycle phases are assigned shifted gamma distributions with different sets of parameters.
Classification of cell stages
The hypothesized neurogenic cascade consists of i) multiple cell types as different states, and ii) probabilities of cells progressing from one type to another as transition rates between states. Our model includes the following cell types: primary NSCs, ANPs, neuroblast/immature neurons (NBs), granule cell neurons (GCs) and apoptotic cells (Apop).
NSC (neural stem cell)
NSCs provide the ultimate influx of newborn ANPs, which massively proliferate to drive the entire cell population to produce mature granule cells. The majority of NSCs are quiescent while activated ones divide asymmetrically to enrich the pool of newborn ANPs. The baseline NSCANP influx can be modeled as a homogeneous Poisson Process. Even though the influx rate may change as the animal ages, we assume it to be fixed at 1monthofage, as all our experimental data are acquired at this age. To quantify the number of proliferating NSCs, we apply BrdU pulseandchase experiment, where 150 mg/kg BrdU is given intraperitoneally to the mouse. Its half life is about 15 min; thus all activated NSCs that are in the S phase are labeled. We call these BrdU+ NSCs  labeled NSCs. Encinas et al., (2011) indicate that each newly activated NSCs proliferates three times to produce three ANP progeny asymetrically and eventually becomes an astrocyte.
ANP (amplifying neuroprogenitor)
Each newborn ANP proliferates several times (~2.45 on average, estimated by Encinas et al., 2011). After the cell divides for a minimum number of times (≥1), it can either keep proliferating, differentiate into a neuroblast, or die. We model ANP progression by specifying minimal/maximum number of divisions and renewal probability (p), which denotes the probability of an ANP continuing to proliferate after finishing its minimal number of divisions. If an ANP chooses to proliferate, it enters typical cell cycle phases of G_{1}, S and G_{2}M; otherwise, it becomes a nonproliferating ANP which may choose to differentiate to neuroblast (NB) by entering ANPNB transition state or commit programmed cell death (apoptosis or Apop) by entering the ANPApop transition state. Note: We assume that cells that are in either ANPNB or ANPApop stage are still ANPs, but nonproliferating.
NB (neuroblast)
NBs are nonproliferating cells that are differentiated from ANPs. The celll duration in NB stage is relatively long (2–10 days) and eventually, each NB may choose to differentiate to immature neuron or enter apoptosis.
(immature neuron)
Similarly as NB, any immature neuron lives for a period of duration time and at the end, differentiates into a mature neuron (granule cell) or dies.
GC (granule cell)
GCs are fully differentiated and mature neurons of the dentate gyrus that remain in dentate gyrus to form neuronal connections with existing neurons of the dentate and hippocampal circuitry. Once a GC is formed, it cannot die or differentiate anymore.
ANPNB
This is an intermediate state of transition time from a nonproliferating ANP to a NB, where the duration is modeled in the same way as duration of any other cell type, as a shiftedgamma distribution.
ANPApop
Is an intermediate state of transition time for a nonproliferating ANP to become an apoptotic cell.
Apop (apoptotic cell)
We assume that apoptosis may occur at the end of any cell stage (G_{1}, S, G_{2}M, NB, IN) along the neurogenic cascade, except for the granule cells. After BrdU pulse, all cells that are dividing and in S phase will be labeled. From the observation of apoptoticBrdU cell labeling curve, there are no labeled apoptotic cells observed at either 2 h or 12 h (Sierra et al., 2010) (Table 1), indicating that proliferating NSCs and their first progeny do not undergo apoptosis or that the removal of the apoptotic cells at these times is so fast that it escapes detection. The estimated duration of ANP_G_{2}M phase is about 2 h [41]. A proportion of newly labeled cells that are in their final allowed division and also transiting from S to G_{2}M would be captured by BrdU. These observations can indicate the existence of ANPApop stage, otherwise the apoptoticBrdU cells should be seen at 2 h or 12 h after BrdU injection. Also, they imply that the apoptotic rate of cells in either S or G_{2}M phase is close to 0; otherwise, cells that are labeled in late S phases can enter cell death immediately after BrdU injection while approaching the end of S or G_{2}M phases.
Model assumptions

The process of proliferation and maturation is driven by a steady influx of generation 1 ANPs (ANP1).

Arriving ANP1 cells enter the G1 phase of the cell cycle

Subsequently, the ANP1 cells proceed through G1, S and G2M phases before they split into two ANP2 cells.

Each ANP2 cell proceeds through the G1, S and G2M phases before it splits into two cells, each of which may become a NB or a ANP3 cell.

Each ANP3 cell proceeds through the G1, S and G2M phases before it splits into two cells, each of which becomes a NB.

NBs exist for the time needed for them to become neurons.

At each cell cycle phase, cells may enter apoptosis.

Apoptotic cells are quickly engulfed by MG and eliminated.

The transit times through the G1 phase is exponentially distributed with expected value T_{1}, whereas transit times through phases S and G2M are deterministic equal to T _{2} and T _{2} respectively.

The lifetime of the ANP is exponentially distributed with expected value T _{4}, whereas the lifetime of the NB is exponentially distributed with expected value T _{5}.
We denote a and b as the minimum and maximum number of divisions of each newborn ANP, where a is the required minimum number of divisions and b is the maximum allowed number of divisions. We further denote p as the renewal probability of each ANP (probability of proliferating after dividing a times) and denote X as the random variable of number of progeny produced by each new born ANP. Therefore, we obtain P(X=2^{a})=1−p, P(X=2^{a+i})=p ^{i}(1−p), for 1≤i≤b−a−1 and P(X=2^{b})=p ^{b−a}, and \(E(X)=2^{a}(1p)+\sum _{i=1}^{ba1}p^{i}(1p)\) 2^{a+i}+2^{b} p ^{b−a}. For a<b, the expected number of ANP divisions can be derived as log _{2} E(X). If a=b, the expected number of division is a.
Modeling of the neurogenic cascade using the BellmanHarris branching process
Our goal was to build a model of the early stage neurogenesis, from NSCs to newborn neurons, considering the factors influencing the ANP fate selection and cell death rates, such as transition to a cell cycle, number of divisions before differentiation into NB and probability of cell death at each step. We chose Multitype BellmanHarris branching process to model the neurogenic cell population and resulting BrdU labeling curves. BellmanHarris process is frequently employed to model proliferation of systems of biological cells [38], and in our model, it is necessary to distinguish cells in different cell cycle phases.
The structure of the model, constructed based on our experimental observations [32], is presented in Fig. 4. We model the hierarchical structure with transition probabilities of cells from a stage to the next possible stage. Thus, cell death rates for different cell phases were modeled by the corresponding transit probabilities into apoptosis, where symbol d _{ i } denotes the cell death rate of the cell type i.
If we denote the multivariate pgf of number of particles of all types present in the process initiated by an ancestor of type i with F _{ i }(s,t), we obtain the BellmanHarris integral equation for this scenario as
Differentiating both sides of the equation with respect to s _{ j } and setting s _{1}=s _{2}=...=s _{ I }=1, we may obtain the following equation for the matrix M(t)=[M _{ ij }(t)], where M _{ ij }(t) is the expected number of particles of type j at time t, in the process initiated by the ancestor particle of type i at time 0
where m _{ ij }(t) is the expected number of progeny of type j of a particle of type i, δ _{ ij }=1, if i=j, otherwise, δ _{ ij }=0.
The convolution of functions is defined by the following LebesgueStieltjes integral if A(t) and B(t) are two rightcontinuous functions with locally bounded variation on [0,∞)
Using the convolution notation, the above equation can be expressed as
In the matrix notation, we obtain
where I is the identity matrix and G=diag(G _{1},...,G _{ I }). This is a renewaltype equation, which has a unique solution of locally bounded variation if G(0)=0, expressed by the infinite series
and it yields the fundamental solution of the mathematical modeling of the neurogenic cascade as a matrix function of time for the number of cells of each type. Here, M _{ ij }(t) is the expected number of cells of type j at time t, in the process initiated by an ancestor cell of type i at time 0, T=diag(T _{1},...,T _{ I }) denotes the diagonal matrix with the distribution function of lifetime (or duration) T _{ i } of each cell i, m is the transition matrix and m _{ ij }(t) is the expected number of progeny of cell of type j produced by a cell of type i, obtained via the multivariate pgf of numbers of progeny produced by the type i cell, and I is the identity matrix.
Based on the experimental observation and model assumptions, we have the transition matrix m as (e.g. when minimum/maximum number of ANP divisions are 1 and 3, respectively)
where AN represents ANPNB stage, AA is for ANPApop stage, \(\overline {d_{i}}=1d_{i}\), \(x^{*}=2\overline {d}_{G_{2}M}(1p)\overline {d}_{ANP}\), y ^{∗∗}=\(2\overline {d}_{G_{2}M}(1p)d_{ANP}\), \(z^{*}=2\overline {d}_{G_{2}M}\overline {d}_{ANP}\) and \(w^{**}=2\overline {d}_{G_{2}M}\) d _{ ANP } (d _{ ANP } is the cell death rate of nonproliferating ANPs).
Furthermore, to model the NSC to ANP influx, we assume that any ‘arrival’ of a new ANP is independent of all previous ‘arrivals’ and the number of new ANPs arrived during a period of time is only dependent on the length of that period times the intensity of the influx, λ. Thus, such process is a Poisson process with intensity parameter λ, and the probability that the number of new ANPs arrived during a time unit (u) being equal to n is expressed as
where N(t) is the number of new ANPs arrived before time t and N(t+u) is the number of the new ANPs arrived until time t+u.
Modeling and simulation of cell labeling curves
The experimental data for our computational modeling have been obtained in three independent sets of timecourse labeling experiments of pulse BrdU injection, with derived curves of total and partial cell counts at different times of measurements: 1) Single BrdU pulseandchase was used to quantify total BrdU+ cell populations over a 32 day period and total number of apoptotic cells were used from the published manuscript [32]. Quantification was done at 12 different timepoints (t=2h r,12h r,1d,2d,3d,4d,8d,11d,15d,18d,22d,32d, assuming t=0 at the moment of BrdU injection) (Table 2). 2) Single BrdU pulseandchase was used to quantify BrdU+ NSCs and ANPs using GFAP and GFP immunostaining to differentiate between them. BrdU+ GFP+ NSCs were identified by their localization in the SGZ, radial GFAP+ process, and fine eGFP+ terminal arborizations in the molecular layer. BrdU+ GFP+ ANPs were identified by their localization in the SGZ, round morphology with no processes, and no GFAP staining. Quantification was done at t=2h r,1d,2d,4d, and 8d (Table 3). 3) Single BrdU pulseandchase was used to quantify NB, IN, and GC using DCX and NeuN immunostaining and morphology. Newborn NBs were BrdU+ DCX+ NeuN or NeuN+ round cells with small processes. Newborn GC were BrdU+ DCX Neu+ mature neurons within the granule cell layer. Quantification was done at t=1d,2d,4d,8d,15d (Table 3). In all experiments, mice were 1 month old at the time of BrdU injection (N = 25 mice per timepoint).
Given the estimated number of cells during the Sphase in each stage at the beginning of BrdU injection, we may calculate the number of labeled cells of each type at any moment by Eq. (1). However, solving it in analytical form is cumbersome. An approach alternativee to computationally producing the BrdU labeling curves is the eventbased simulation. Assuming that we have computed the numbers of cells in different stages at the moment of BrdU injection (t=0), we trace the fate of labeled cells at unit time points by recording their behaviors. Briefly, a series of random numbers are generated for the random times for which labeled cells stay in particular stages and the probability of the cells transiting to the next stage, until the cells enter apoptosis or become a matured neuron. Beginning to trace the entire process from t=0, we reproduce the labeling curves in silicoby accumulating the fates of all labeled cells up to particular moments (e.g. times of measurements in experiments).
A simulation program carrying out tasks outlined above has been written in the Python programming language. This program computes distribution of initial cell population and generates BrdU labeling curves. Figure 5 illustrates the eventbased simulation scheme of the dynamics of a neuroprogenitor cell.
Parameter search and goodnessoffit
To discover parameter combinations that can best fit the experimental data, we adapted a genetic algorithm as the searching heuristic since the parameter space is too complex to be searched by enumeration of all possible combinations. A genetic algorithm is a searching heuristic that mimics the process of natural evolution. It is used to generate optimized solutions to search problems in complex nonlinear systems. Each parameter range is encoded by a bit vector of length a, yielding 2^{a} possible values. An initial pseudo “population” was created by setting X randomly chosen parameter combinations as X “individuals”. The value of each modeling parameter in any “individual” has been converted to the binary format to become a 01 sequence. During the search, any parameter set is evaluated by computing variance weighted least square error \(\sum \frac {(ES)^{2}}{\sigma ^{2}}\) to test how well the simulated results fit the experimental labeling curves. E and σ ^{2} are mean and variance of experimental data at a given time point, whereas S is the corresponded simulation result. \(\sum \) sums over all available time points of measurements. The list of model parameters is shown in Table 1. Some of them can be experimentally estimable, whereas most are not observable.
Results
Cell proportions
At time points t= 1, 2, 4and 8 days, the sum of observed proportions of all cell types (NSCs, ANPs, NBs, GCs and apoptotic cells) is greater than 1 (Tables 2 and 3). This discrepancy is due to the nonspecific labeling or identification of cells during transitions from ANPs to NBs. Cells that are in the intermediate stage can be labeled by both ANP and NB markers. Thus, we assume that the estimated proportions on NSCs and GCs are realistic, while those of ANPs and NBs at intermediate time points are inflated due to nonspecific labeling. Those quantifications need to be adjusted. At early (2hr) or late time points (15 and 32 days), we retained the original data since at those time points the labeled cells are NSCs and GCs. At all other timepoints, the labeled cells are either ANPs or NBs, exclusively.
For any time point t (t = 1, 2, 4 or 8 days), we adjusted the sum of proportions of all types of BrdU+ cells (NSCs, ANPs, NBs, GCs and apoptotic cells) to be equal to 1 (impact of BrdU+ astrocytes is negligible since the proportion of BrdU+ astrocytes is very small when t≤8 days [41]). The proportion of the excessive amount of cells (denoted by d _{ t }) between ANPs and NBs is equal to the sum of the observed proportions minus 1. We denote α _{ t } to be the ratio defined so that α _{ t } d _{ t } is the proportion of double labeled cells that belong to ANPs, whereas (1−α _{ t })d _{ t } belong to NBs.
Let X=(X _{1},X _{2})^{T} and X ^{′}=(X1′,X2′)^{T} as two random vectors to represent the proportions of ANPs and NBs before and after transformation, respectively. Therefore, we obtain
and
where \(\mathbf {A}=\left (\begin {array}{cc} \alpha _{t} & \alpha _{t}1\\ \alpha _{t} & 1\alpha _{t} \end {array}\right)\), \(\mathbf {B}=\left (\begin {array}{c} (1\alpha _{t})d\\ \alpha _{t} \end {array}\right)\).
Thus, we can compute the mean and variance of X ^{′} as
where Σ[X ^{′}] and Σ[X] are covariance matrices of X ^{′} and X with Σ[X]=(V[X _{1}],V[X _{2}]) by assuming that C O V[X _{1},X _{2}]=0. Using the equations above, we calculated the adjusted estimates of means and sems of ANP and NB cell proportions (Table 4). Note that for any interval time point t (t=1d,2d,4d,8d), α _{ t } is assumed to be 1/2 since there is no prior knowledge about it.
Transforming cell proportions to cell counts
For the optimization algorithm, the total number of BrdU+ cells and the estimated number of BrdU+ cells of each specific type are required to evaluate the goodnessoffit. Nonzero data points expressed as proportions (random variables; Table 3) were transformed back to cell counts with recalculated means and variances.
For any time point t, we assumed that the count of total number of BrdU+ cells in an animal is a normally distributed random variable \(Y\sim N\left (\mu _{Y},\sigma _{Y}^{2}\right)\). Thus, for a size of n _{ Y } samples, we obtain that \(\hat {\mu _{Y}}=\overline {Y}\) and \(\hat {\sigma _{Y}^{2}}=S_{Y}^{2}\) (data from Table 1 were used to estimate \(\hat {\mu _{Y}}\) and \(\hat {\sigma _{Y}^{2}}\)).
For any specific cell type i (e.g. ANPs), if we denote X as the number of BrdU+ type i cells at time t, we have
where we assume that the observed proportion of type i cell, P, is Gaussian distributed, s.t. \(P\sim N\left (\mu _{P},\sigma _{P}^{2}\right)\). For a size of n _{ P } samples, \(\hat {\mu _{P}}=\overline {P}\) and \(\hat {\sigma _{P}^{2}}=S_{P}^{2}\) \(\left ({\vphantom {\hat {\sigma _{P}^{2}}}}\hat {\mu _{P}}\right.\) and \(\hat {\sigma _{P}^{2}}\) are estimated from data shown in Table 4).
Assuming that Y and P are independent random variables, we obtain that
and
Since
we obtain
Therefore, we estimate the mean and variance of the number of BrdU+ labeled cells of any type i (E[X] and V[X]) by equations below:
A summary of data in estimated cell numbers is shown in Table 5.
Initial distribution of cell population
Before deriving the model to simulate BrdU labeling curves, we focused on the quantification of the cell populations at the moment of the BrdU injection and aimed to obtain the distribution of cell populations at t=0.
In matrix M (Eq. (1)), the fundamental solution of the model M _{ ij } as a function of time, allows finding the number of cells at time t in compartment j, given that the population was seeded by a single cell in compartment i. However, under physiological conditions, the system is fed by a steady influx of freshly activated ANPs. Under such assumption, the number of of cells of each cell type at time t is given by
Since the system is fed only by newborn ANPs at their first division in G _{1}phase, only the top row of matrix \(\widetilde {M}(t)\), denoted by \(\widetilde {M}^{(1)}(t)\), is required. In labeling experiments, we treat 1 month old mice with BrdU injections. We assume that the snapshot of the neurogenic cell population in the animal’s brain is under a “steady state” at the moment of injection. Thus, if we tend with t to infinity in the Eq. (3), we can obtain a stationary distribution of cell numbers in the snapshot, which yields numbers of cells of different types, given by
where π is the vector of numbers of cells in all modeling compartments (see Matrix Eq. (2) for the description of modeling compartments), diagonal matrix E[T] has entries equal to expected duration times in all modeling compartments, (·)^{−1} denotes the inverse of a matrix, (·)^{(1)} denotes the top row of the matrix and λ,T,m,I,and∗ are as previously defined. Although in the long run the intensity parameter λ declines as the animal gets older, we assume that in a short period of time, within which the snapshot of the BrdU injection occurs, the influx rate of NSCs to newborn ANPs is constant.
BrdU labels all cells that are in Sphase, thus we know how many cells are labeled at the moment of injection, which is equal to π _{ s }, where subscript s stands for the compartment or a set of compartments which represents the cells in the Sphase. Under the assumption that in vivo descendants of labeled cells remain labeled (BrdU dilution is negligible), the BrdU pulse labeling curve of the number of labeled cells at a given time is equal to π _{ s } M(t). This expression is technically true only under the assumption that labeled cells concentrated at the beginning of the Sphase, which does not make much difference for times longer than the joint duration of S and G _{2} M. In a real situation, the time remaining for each cell in Sphase at any moment is a random variable. Additionally, the kfold convolution of matrix (T m) in the matrix M(t) (Eq. (1)) becomes analytically too complicated as k increases. Therefore, as explained before, instead of obtaining the timecourse labeling curve analytically, we decided to generate it by simulation in a more convenient and straightforward manner.
Simulation of BrdU labeling curves and parameter search
To obtain the set of model parameters that yields best fit to the experimental observation of BrdU pulseandchase labeling curves, we applied the genetic algorithm as the search heuristic in the simulation program. Table 6 lists model parameters and their assigned ranges of possible values. For example, the maximum number of ANP divisions ranges from 2 to 8, renewal probability of ANPs ranges from 0 to 1, expected duration of ANP Sphase ranges from 5 to 12 hr, shape parameter of distribution of ANP Sphase duration ranges from 5 to 40, minimum duration of ANP Sphase ranges from 1 to 4 and apoptotic rate for ANP Sphase can range from 0 to 0.99.
Simulation results of BrdU labeling curves that best fit our data are illustrated in Fig. 6. Parameter values that yield such fits were obtained by the genetic algorithm, which is implemented for the purpose of searching the parameter space and optimizing the goodnessoffit function. BrdU pulse labeling experiments have provided us with 40 nontrivial (nonzero) independently measured experimental data points. Nineteen model parameters are varied during the parameter search. Figure 7 demonstrates that the residuals were equally distributed along x axis and showed no systematic trend, which suggests that the model fit is good.
Among all model parameters, apoptotic rates are the most critical ones since they have not been estimated in prior earlystage hippocampal neurogenesis studies. Based on our model and simulation results at each cell state during early stages of hippocampal neurogenesis in 1monthold mice, we estimate that the apoptosis rates are low in proliferating ANPs whereas once ANPs become nonproliferating, about one third of them undergo apoptosis (Fig. 8). During the NB stage, apoptosis reaches maximum. A vast majority of the NBs die (97% undergo apoptosis) and only a few of them (estimated about 3%) will differentiate into the GCs. NSCs do not undergo apoptosis [32] and once a NSC is activated, it undergoes a number of asymmetric divisions after which it eventually becomes an astrocyte.
Prediction of dynamics of neurogenesis under reduced apoptosis
Taking values of model parameters from data fitting results, we carried out additional simulations to predict the overall changes in the BrdU labeling curves by inhibiting apoptosis (reducing apoptotic rates). While apoptotic rates at all cell stages were consistently reduced by a hypothetical amount (25%, 50%, 75% or 100%) all other model parameters remain unchanged. From predicted BrdU labeling curves depicted in Fig. 9, we observed that at the end of 32 days, the total number of BrdU+ cells and the number of BrdU labeled granule cells increased 3.4 and 11.5 fold, respectively, when apoptotic rates were reduced by 25% only. These numbers continue to increase sharply if apoptosis can be reduced even further. Under the extreme scenario, when apoptosis can be completely inhibited, the simulation results indicate that 14.3 times more of total BrdU+ cells and 61.0 times more of BrdU+ GCs are expected as the net outcome of neurogenesis, compared with the case when physiologic apoptotic rates are employed. Our study thus indicates that reducing apoptosis in any amount substantially increases adult hippocampal neurogenesis.
Discussion
In this study, we developed a computational model that as accurately as possible reflects the neurogenic cascade and specifically, apoptosis. Our goal was to estimate apoptotic rates at each stage of the neurogenic cascade, the distribution of cell type duration – including apoptosis, and the renewal probability of ANPs. We reasoned that these parameters were most important to design targeted experimentation to improve survival of newborn cells and net outcome of hippocampal neurogenesis. Since there is unavoidably a large amount of model parameters providing challenges and obstacles to unbiased estimation, we employed immunohistochemistry and statistical computation approaches to combine experimental and computational data. Furthermore, we computationally estimated experimentally unobservable parameters, such as the probability of ANP to proliferate and the rates of cells at different stages undergoing apoptosis. Our model indicates that apoptosis is low in the ANP stage and high in the NB stage. Regardless of origin, apoptotic cells have a short life, estimated to be around 1.4hrs. ANPs are predicted to divide 1–4 times; however, their renewal probability is low, at 0.1. Finally, the NB stage has the largest variance of the transit time. None of the estimates could be derived experimentally, and thus, our computational model represents a foundation upon which we can design novel biological experiments to increase neurogenesis based on targeted action on ANPs and newborn cell survival. Encinas et al. (2011) carried out labeling experiments (both single and cumulative labeling) to study adult hippocampal neurogenesis. They modeled neurogenic cascade similarly as we do, although they used 2 month old mice whereas we used 1 month old ones. While Encinas et al. (2011) determined division and duration related parameters, they did not infer any information on apoptosis. Their results were calculated from inferring the decay rate of each type of cell over a long period of time (800 days). In comparison, the parameter values that yield best fit in our study were comparable with respect to expected cell durations (Table 7). In addition, our model and simulation approach are able to provide estimates on apoptotic rates, minimum durations, shapes of duration distributions, and number of NSC and ANP divisions. More recent works [42, 43] investigates the regulatory mechanisms of neurogenesis, based on knockout experiments, which modify the dynamic behaviour of this process. Evaluating these knockout is a nontrivial task owing to the complicated nature of the hippocampal neurogenic niche. Unlike the model proposed herein, they model neurogenesis as a multicompartmental system of ordinary differential equations based on experimental data. To analyse the results of knockout experiments, they investigated how changes of cell properties, based on cells labelled by the cell division marker BrdU. Among other, they found that changing cell proliferation rates or the fraction of selfrenewal, may result in multiple time phases in the response of the system, such as an initial increase in cell counts followed by a decrease. Because of different experimental setup and modeling framwork used, these results are not directly comparable to ours. One of the obstacles is the difficulty in observing and recording the fates of the individual cells in vivo.
Conclusion
In sum, our computational model of the adult neurogenesis provides new information on the early stages of this phenomenon. It is our hope that the estimates of the properties of ANPs, NBs, and apoptotic cells will guide biological investigations and development of better experimental tools to utilize this unique process for the benefit of human health.
Abbreviations
 ANP:

Amplifying neuroprogenitor
 Apop:

Apoptotic cell
 BrdU:

Bromodeoxyuridine
 BSA:

Bovine serum albumin
 CFP:

Cyan fluorescent protein
 DCX:

Doublecortin
 GC:

Granule cell
 GFAP:

Glial fibrilary acidic protein
 GFP:

Green fluorescent protein
 IN:

Immature neuron
 NB:

Neuroblast
 NPCs:

Neural stem and progenitor cells
 NSC:

Neural stem cell
 PBS:

Phosphatebuffered saline
 PFA:

Paraformaldehyde
 SGZ:

Subgranular zone
References
 1
Lois C, AlvarezBuylla A. Proliferating subventricular zone cells in the adult mammalian forebrain can differentiate into neurons and glia. Proc Natl Acad Sci. 1993; 90(5):2074–7.
 2
Cameron H, Woolley C, McEwen B, Gould E. Differentiation of newly born neurons and glia in the dentate gyrus of the adult rat. Neuroscience. 1993; 56(2):337–44.
 3
Sierra A, Encinas JM, MaleticSavatic M. Adult human neurogenesis: from microscopy to magnetic resonance imaging. Front Neurosci. 2011; 5:47.
 4
Semerci F, MaleticSavatic M. Transgenic mouse models for studying adult neurogenesis. Front Biol. 2016; 11(3):151–67.
 5
Dupret D, Revest JM, Koehl M, Ichas F, De Giorgi F, Costet P, et al.Spatial relational memory requires hippocampal adult neurogenesis. PloS one. 2008; 3(4):e1959.
 6
Zhang CL, Zou Y, He W, Gage FH, Evans RM. A role for adult TLXpositive neural stem cells in learning and behaviour. Nature. 2008; 451(7181):1004–7.
 7
Deng W, Aimone JB, Gage FH. New neurons and new memories: how does adult hippocampal neurogenesis affect learning and memory?Nat Rev Neurosci. 2010; 11(5):339–50.
 8
Gould E, Beylin A, Tanapat P, Reeves A, Shors TJ. Learning enhances adult neurogenesis in the hippocampal formation. Nat Neurosci. 1999; 2(3):260–5.
 9
Van Praag H, Kempermann G, Gage FH. Running increases cell proliferation and neurogenesis in the adult mouse dentate gyrus. Nat Neurosci. 1999; 2(3):266–70.
 10
Kempermann G, Kuhn HG, Gage FH. More hippocampal neurons in adult mice living in an enriched environment. Nature. 1997; 386(6624):493.
 11
Revest J, Dupret D, Koehl M, FunkReiter C, Grosjean N, Piazza P, et al. Adult hippocampal neurogenesis is involved in anxietyrelated behaviors. Mol Psychiatry. 2009; 14(10):959–67.
 12
Parent JM, Timothy WY, Leibowitz RT, Geschwind DH, Sloviter RS, Lowenstein DH. Dentate granule cell neurogenesis is increased by seizures and contributes to aberrant network reorganization in the adult rat hippocampus. J Neurosci. 1997; 17(10):3727–38.
 13
Mohapel P, Ekdahl CT, Lindvall O. Status epilepticus severity influences the longterm outcome of neurogenesis in the adult dentate gyrus. Neurobiol Dis. 2004; 15(2):196–205.
 14
Sierra A, MartínSuárez S, ValcárcelMartín R, PascualBrazo J, Aelvoet SA, Abiega O, et al. Neuronal hyperactivity accelerates depletion of neural stem cells and impairs hippocampal neurogenesis. Cell Stem Cell. 2015; 16(5):488–503.
 15
Abiega O, Beccari S, DiazAparicio I, Nadjar A, Layé S, Leyrolle Q, et al. Neuronal hyperactivity disturbs ATP microgradients, impairs microglial motility, and reduces phagocytic receptor expression triggering apoptosis/microglial phagocytosis uncoupling. PLoS Biol. 2016; 14(5):e1002466.
 16
Erramuzpe A, Encinas J, Sierra A, MaleticSavatic M, Brewster A, Anderson AE, et al. Longitudinal variations of brain functional connectivity: A case report study based on a mouse model of epilepsy. F1000Research. 2015;4.
 17
Eisch AJ, Barrot M, Schad CA, Self DW, Nestler EJ. Opiates inhibit neurogenesis in the adult rat hippocampus. Proc Natl Acad Sci. 2000; 97(13):7579–84.
 18
Eisch AJ. Adult neurogenesis: implications for psychiatry. Prog Brain Res. 2002; 138:315–42.
 19
Noonan MA, Choi KH, Self DW, Eisch AJ. Withdrawal from cocaine selfadministration normalizes deficits in proliferation and enhances maturity of adultgenerated hippocampal neurons. J Neurosci. 2008; 28(10):2516–26.
 20
Duman RS, Malberg J, Nakagawa S. Regulation of adult neurogenesis by psychotropic drugs and stress. J Pharmacol Exp Ther. 2001; 299(2):401–7.
 21
Encinas JM, Vaahtokari A, Enikolopov G. Fluoxetine targets early progenitor cells in the adult brain. Proc Natl Acad Sci. 2006; 103(21):8233–8.
 22
David DJ, Samuels BA, Rainer Q, Wang JW, Marsteller D, Mendez I, et al. Neurogenesisdependent andindependent effects of fluoxetine in an animal model of anxiety/depression. Neuron. 2009; 62(4):479–93.
 23
Duan X, Chang JH, Ge S, Faulkner RL, Kim JY, Kitabatake Y, et al. DisruptedInSchizophrenia 1 regulates integration of newly generated neurons in the adult brain. Cell. 2007; 130(6):1146–58.
 24
Reif A, Fritzen S, Finger M, Strobel A, Lauer M, Schmitt A, et al. Neural stem cell proliferation is decreased in schizophrenia, but not in depression. Mol Psychiatry. 2006; 11(5):514–22.
 25
Weissman IL, Anderson DJ, Gage F. Stem and progenitor cells: origins, phenotypes, lineage commitments, and transdifferentiations. Annu Rev Cell Dev Biol. 2001; 17(1):387–403.
 26
Kempermann G, Jessberger S, Steiner B, Kronenberg G. Milestones of neuronal development in the adult hippocampus. Trends Neurosci. 2004; 27(8):447–52.
 27
Malinow R, Hayashi Y, MaleticSavatic M, Zaman S, Poncer J, Shi S, et al. Introduction of green fluorescent protein into hippocampal neurons through viral infection. 1999.
 28
Malinow R, Hayashi Y, MaleticSavatic M, Zaman SH, Poncer JC, Shi SH, et al. Introduction of green fluorescent protein (GFP) into hippocampal neurons through viral infection. Cold Spring Harb Protocol. 2010; 2010(4):pdb–prot5406.
 29
Kempermann G, Kuhn HG, Gage FH. Genetic influence on neurogenesis in the dentate gyrus of adult mice. Proc Natl Acad Sci. 1997; 94(19):10409–14.
 30
Ge S, Goh EL, Sailor KA, Kitabatake Y, Ming Gl, Song H. GABA regulates synaptic integration of newly generated neurons in the adult brain. Nature. 2006; 439(7076):589–93.
 31
Semerci F, Choi WT, Bajic A, Thakkar A, Encinas JM, Depreux F, Segil N, Groves AK, MaleticSavatic M. Lunatic fringemediated Notch signaling regulates adult hippocampal neural stem cell maintenance. Elife. 2017:6. doi:10.7554/eLife.24660.
 32
Sierra A, Encinas JM, Deudero JJ, Chancey JH, Enikolopov G, OverstreetWadiche LS, et al. Microglia shape adult hippocampal neurogenesis through apoptosiscoupled phagocytosis. Cell Stem Cell. 2010; 7(4):483–95.
 33
Garthe A, Behr J, Kempermann G. Adultgenerated hippocampal neurons allow the flexible use of spatially precise learning strategies. PloS one. 2009; 4(5):e5464.
 34
Dupret D, Fabre A, Döbrössy MD, Panatier A, Rodríguez JJ, Lamarque S, et al. Spatial learning depends on both the addition and removal of new hippocampal neurons. PLoS Biol. 2007; 5(8):e214.
 35
Aimone JB, Wiles J, Gage FH. Computational influence of adult neurogenesis on memory encoding. Neuron. 2009; 61(2):187–202.
 36
Aimone JB, Deng W, Gage FH. Adult neurogenesis: integrating theories and separating functions. Trends Cogn Sci. 2010; 14(7):325–37.
 37
Weisz VI, Argibay PF. A putative role for neurogenesis in neurocomputational terms: inferences from a hippocampal model. Cognition. 2009; 112(2):229–40.
 38
Kimmel M, Axelrod DE. Branching Processes in Biology, 2nd Edn. New York: Springer; 2015.
 39
Encinas JM, Enikolopov G. Identifying and quantitating neural stem and progenitor cells in the adult brain. Methods Cell Biol. 2008; 85:243–72.
 40
Casella G, Berger RL. Statistical inference. vol. 2. CA: Duxbury Pacific Grove; 2002.
 41
Encinas JM, Michurina TV, Peunova N, Park JH, Tordo J, Peterson DA, et al. Divisioncoupled astrocytic differentiation and agerelated depletion of neural stem cells in the adult hippocampus. Cell Stem Cell. 2011; 8(5):566–79.
 42
Ziebell F, MartinVillalba A, MarciniakCzochra A. Mathematical modelling of adult hippocampal neurogenesis: effects of altered stem cell dynamics on cell counts and bromodeoxyuridinelabelled cells. J R Soc Interface. 2014; 11(94):20140144.
 43
Ziebell F, Dehler S, MartinVillalba A, MarciniakCzochra A. Revealing AgeRelated Changes of Adult Hippocampal Neurogenesis. bioRxiv. 2017:112128.
Acknowledgements
We thank the members of MaleticSavatic lab for comments and critical reading of the paper.
Funding
The research was supported in part by the Baylor College of Medicine Microscopy Core (P30HD024064 Intellectual and Developmental Disabilities Research Grant from the Eunice Kennedy Shriver National Institute of Child Health and Human Development). A.L. was supported by the NLM Training Program in Biomedical Informatics (T15LM007093) and the Baylor College of Medicine Medical Scientist Training Program. Publication of this article was funded in part by the McKnight Endowment for Science, the CPRIT grant (RP130573CPRIT) and the NIH grant GM12003301 (M.M.S.), and the NIH grant R01HL128173 and the Polish National Science Center grant DEC 2012/04/A/ST7/00353 (M.K.). The funding bodies did not have any role in the design or conclusions of this study.
Availability of data and materials
All data and materials will be shared in accordance with the NIH Grants Policy on Sharing of Unique Research Resources.
About this supplement
This article has been published as part of BMC Systems Biology Volume 11 Supplement 5, 2017: Selected articles from the International Conference on Intelligent Biology and Medicine (ICIBM) 2016: systems biology. The full contents of the supplement are available online at https://bmcsystbiol.biomedcentral.com/articles/supplements/volume11supplement5.
Author information
Affiliations
Contributions
BL designed the computational model, analyzed the data and wrote the manuscript. AS, JJD, FS performed the biological experiments, collected data and participated in validation of the computational model. AL prepared the final manuscript. MK and MMS designed and supervised all studies, analyzed and interpreted the data, provided financial support, and wrote the manuscript. All authors have read and approved the manuscript.
Corresponding authors
Correspondence to Marek Kimmel or Mirjana MaleticSavatic.
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
All authors declare that they have no competing interests.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License(http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver(http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Li, B., Sierra, A., Deudero, J.J. et al. Multitype BellmanHarris branching model provides biological predictors of early stages of adult hippocampal neurogenesis. BMC Syst Biol 11, 90 (2017). https://doi.org/10.1186/s1291801704683
Published:
Keywords
 Hippocampus
 Adult neurogenesis
 Apoptosis
 Computational modeling
 Multitype BellmanHarris branching process