Link with the continuous time Markov chains

Get Complete Project Material File(s) Now! »

Cell average models (used in Chapter 2 and 3)

The simplest way to describe the changes in cell dynamics is to follow the total cell num-ber along time. The natural associated mathematical formalism is the Ordinary Differential Equations system (ODEs), and the cell number dynamics is then described as a continuous function of time. We introduce the main models using the finite dimension differential sys-tem formalism to describe biological structures. We also present the modeling principles used during this thesis to represent the cell dynamics involved in the formation of an ovarian follicle.
The best-known (and simplest) model of population dynamics is likely to be the Malthus growth model, published in 1798 [4]. In this model, the change in speed of a growing pop-ulation is proportional to the population size (see Table I.1). A generalized version of this model is the Birth–Death model, where the population growth speed depends on both the population, and time-dependent birth and death rates. The Malthus growth and Birth–Death models belong to the linear model category in which each individual cell is independent from the other cells. Since the ODE formalism only represents the total cell number, it may be unnatural to think about independency between cells. This property is a direct conse-quence of the additivity property verified by the linear ODE solution: suppose that we have a population of N0 cells at time t = 0 following the Malthusian growth of parameter α. We consider the trajectory of each individual in the population and thus solve N0 systems following the Malthus law of parameter α and starting from 1 cell. Each of these systems verifies N(t) = eαt, for all t ≥ 0. Summing those N0 systems (by additivity property of the linear ODE solution), we obtain that the total cell number is, for all time t ≥ 0, the expected function N(t) = N0eαt, which is the trajectory of population following a Malthus growth and composed initially of N0 cells.
The cell interactions are actually important mechanisms that allow the development of a tissue. These various mechanisms are generally the result of local events: for example, chemical factors exchanged to help develop a tissue or growth factors that are collected from the local cell environment to facilitate the tissue development. The integration of interaction phenomena leads to non-linear models that are often difficult to analyze or even simulate. With regard to the ODE formalism, the cell interactions can be represented by a non-local term consisting of nonlinear functions depending on the population size. For example, in the Gompertz or Verhulst models, the population grows modulated by an intrinsic growth rate until it reaches a steady state K corresponding to the maximal capacity of the population (see Table I.1). Another case of cell interactions is treated in Chapter 2 where we propose a model dedicated to the ovarian follicle activation.
One may want to represent different cell types such as proliferative versus quiescent cells. This can be done with the compartmental models, which consist in decomposing the cell population into subcategories representing the selected cell types that explain the whole pop-ulation dynamics. For example, the famous Lotka–Volterra model [5, 6], that represents birth and competition events occuring between two species, is a good example of compartmental model with interaction phenomena. As regards to cell dynamics, the biological ODE system models presented in [7, 8] aim to analyze the differentiation and proliferation of cells involved in the primary immune response. In [9], a model reproducing the sequence of divisions that a pool of progenitor cells undergoes during the neurogenesis phase is presented. In [10], a time-dependent compartmental model is proposed to describe the cellular changes in the somatic cells during the terminal phase of the ovarian follicle development. The model com-partements represent the proliferative, the differentiated (exited from the cell-cycle) and the apoptotic cells.
Once the model is built, some dynamic characteristics can be inferred from its outputs such as the doubling time (for example, equals to ln(2)/α for the Malthus growth model ). Another way to characterize the trajectories of the built-model is to the search for some remarkable solutions such as the Exponentially Stable States.
Exponentially stable state analysis: general theory
One can look for solutions of the form t 7→N eλt. In particular, when λ is positive and max-imal (Malthus parameter), such solutions are called the Exponentially Stable States (ESS) and N ∈ RN is then called the stable state vector. The ESS indicates that the overall cell number evolves exponentially following a Malthusian growth.
Once a remarkable solution of a dynamical system has been found, it may be interesting to investigate to what extent this trajectory is representative of the system: do all the trajectories of the dynamical system converge towards this particular trajectory (attractive solution)? If so, after how long? A technique based on a function called entropy has been used to measure the difference between two solutions of a dynamical system over time. If this entropy function has appropriate properties, then the convergence of the two solutions towards each other can be shown. This method is called the Generalized Relative Entropy (GRE) [11, 12]. We introduce below some elements to understand its application to the specific case of a finite dimensional system inspired by the work done in Section 6.3 of [12].
Suppose that the dynamics of the whole population n(t) := (nk(t))k∈J1,KK follows the ODE system:
d n(t) = An(t), n(t = 0; θ) = n0 (I.1) where A ∈ MK (R)4 is a positive matrix, i.e. for all 1 ≤ i, j ≤ K, the element Aij of matrix A verifies Aij > 0. We can thus apply the Perron-Frobenius theorem (see [13] for the positive matrix case) and deduce that A has a first (maximum) eigenvalue λ > 0 associated with a positive right eigenvector N ∈ RK , and a positive left eigenvector φ ∈ RK
Proposition I.1 (Proposition 6.6, extract from [12]). Let λ be the first (positive) eigenvalue associated with matrix A and H be a convex function on R, then the solution to Eq. (I.1) satisfies dtdH[n(t)e−λt] ≤ 0.
The use of the right and left eigenvectors N and φ, the convexity of function H and the positivity of matrix A are crucial to prove that the GRE decreases with respect to time. On the contrary, the positivity of the Malthus parameter λ is not required here to prove the GRE decays. In Chapter 3, we will deal with a case of non-positive matrix A.
From the decay of the entropy H[n(t)e−λt] with respect to time, one can deduce that the trajectory t 7→n(t)e−λt approaches the right eigenvector N in a given norm. For example, taking H(x) = x2, a convergence with square norm can then be deduced.
One can consider an entropy approach to show the convergence of all trajectories n to the (intermediate) state eλtN. Since the left and right eigenvectors are not positive, the entropy function proposed in the usual theoretical framework (Eq. (I.2)) cannot be used in this case and must be adapted. We tackle this issue in Chapter 3 (in the more general context of infinite differential systems).
In this subsection, we have just seen that a first way to model and analyze a biological system consists in selecting suitable structuring states that are markers of a biological function (cell types) and that have a significant role in population dynamics. The result is a differential system with a finite size of K. When K becomes very large, the analysis of the model often becomes very complicated and finally calls for another formalism where the structuring states are described by continuous variables.

Structuring physiological variables at the population level (used in Chapter 3 to 5)

Structured population models first appeared almost two hundred years after the first average population models. The structuring state generally chosen corresponds either to a physio-logical characteristic that changes continuously over the life of the individual (size, age, size increment, etc.) or fixed at birth (genetic mutations). Some of those models were constructed from the discrete differential system by applying a method named continuum limit (see for instance [16, 17]). However, the models were built most of the time from scratch. During this thesis, we focus on age-dependent models and therefore present a non-exhaustive state-of-the-art of the different approaches considered in that case.
The McKendrick–VonFoerster and renewal equations
The age-structured models have the particularity of assuming two biological properties: the dynamic behavior of each individual is independent of that of his or her parents and an individual’s ability to reproduce or die does not depend on past experience [18]. The initial genesis of age-structured models is too complex to be described here, but the reader can refer to [19] for a complete presentation.
We solely introduce the partial differential equation (PDE) model named as McKendrick-VonFoerster (MK-VF) equation. First introduced by McKendrick in 1926 for demographic studies [20], it was reintroduced almost 30 years later by Von-Foerster in a cellular biology context without any link to the first model [21]. It is now known as the McKendrick– VonFoerster equation [19], and consists of a non-conservative transport equation with a non-local term at the boundary condition:
∂tρ(t, a) + ∂aρ(t, a) = −δ(a)ρ(t, a)
ρ(t, 0) = R0+∞ µ(a)ρ(t, a)da,
ρ(0, a) = ρ0(a)
where δ and µ are continuous positive functions defined on R+ corresponding to the death and birth rates, respectively. The age-structured population is represented by the density function ρ(t, •) ∈ L1(R+). The initial condition ρ0 is usually assumed to belong to the L1(R+) space. Another classical assumption is often made on the death rate:
Z a (I.4)
limδ(s)ds = +∞.
a→+∞ 0
Hence, the probability that an individual is still alive at age a verifies exp(− R0a δ(s)ds) and assumption (I.4) guarantee that this individual eventually dies.
This model is a special case where the structuring variable evolves at the same speed as time. There exist other models representing the evolution of temporal variables such as mat-uration that do not evolve at the same speed as time (see for instance the cell maturation models [22, 23, 24]). Interestingly, the Leslie model [25] which corresponds to an age and time-discrete version of the MK–VF model was also proposed after the first introduction of the MK–VF. A retrospective link between these two models is proposed in [16] where the author uses the continuum limit to move from Leslie’s model to the MK–VF model.
The MK–VF model is mainly a demographic model. In this thesis, we focus on its “cell model” counterpart known as the renewal equation. The renewal equation is a special case of the MK–VF where the death rate is linked to the birth rate:
∂tρ(t, a) + ∂aρ(t, a) = −b(a)ρ(t, a),
ρ(t, 0) = 2 R0+∞ b(a)ρ(t, a)da, (I.5)
ρ(0, a) = ρ0(a),
where both ρ0 and ρ(t, •) belongs to the L1(R+) space. The function b is supposed to be positive and represents the division rate. It usually belongs to the C(R+) space (the set of continuous functions) but we also consider the case of non-negative functions belonging to the space Cc(R+) (the set of continuous compactly supported functions) in Chapter 5. In the same way as the MK–VF, the probability that an individual has not divided at age a verifies P [τ ≥ a] = exp − b(s)ds .
If one wants to model a population in which all cells eventually divide (non quiescent cells), the following hypothesis is necessary Z a (I.6) limb(s)ds = +∞. a→+∞ 0
Both the MK–VF and the renewal equation have been extensively studied and modified since their introduction. Regarding the theoretical studies, general results based on operator theory can be found in [26, 18]; [19, 27] provide analytical solutions using the method of characteristics. In the case of the renewal equation, this representation formula has the form
( ρ(t a, 0) exp ( R −a b(s)ds) when a t.
ρ(t, a) = ρ0 (a − t) exp − aa t b(s)ds when a ≥ t (I.7)
− − R 0 ≤
(see the details of the proof in Chapter 5). Note that this formula still depends on the ρ func-tion. However, this particular shape helps to design a numerical scheme [28]. Other methods are also available: escalator boxcar train (see for instance [29, 30]), finite-difference method (see for instance [31, 32]) and finite-volume method [33], presented in Chapter 3. A review of the main numerical schemes available for the age-structured equations can be found in [32].
More recently, studies have been done on measure solutions. They consider measure ini-tial conditions such as Dirac mass instead of the classical L1 function. Such conditions allow to consider synchronized initial conditions and will be more widely discussed in Chapter 5. In [34, 35, 36, 37], the authors study the existence and uniqueness of such solutions while numerical approaches are proposed in [38, 39].

The Lotka integral equation

Even before the first appearance of the age-structured PDE model, Sharpe and Lotka intro-duced in 1911 a model dedicated to the age distribution of a population [40]. Their aim was to represent the age fluctuation based on demographic observations: “certain age-distributions will practically never occur” [40]. This model links the number of births per unit time B to the fraction l of individuals surviving to age a : B(t) = B(t − a)l(a)b(a)da.
The function b is the birth rate per capita for mother of age a. Based on the other observation that “there must be a limiting stable type about which the actual distribution varies, and towards which it tends to return if through any agency disturbed therefrom”, they assume an exponential solution B(t) = Qert and deduce the Fredholm integral equation, called the characteristic equation [41], t e−ral(a)b(a)da = 1. 0
Sharpe and Lotka laid the foundations for the well-known result of the convergence of all the trajectories to an exponential stable state, but did not provide any clues for the proof. The proof of this result was ultimately demonstrated in 1941 by Feller [42] and calls for tools be-longing to the renewal theory approach. Another proof based on group theory was proposed later in [41].
The true “renewal equation” term refers to the integral equation of the form Z t (I.8) u(t) = g(t) + u(t − x)f(x)ds.
Eq. (I.8) is a linear Volterra equation of the second kind where f is sometimes called the kernel and g is a source term. This type of equation appears under different forms mainly in population theory, the theory of industrial replacement and self-renewing aggregates, and are usually analyzed and solved by Laplace transform techniques. The term renewal equation used above and referring to a MK–FV like equation is thus a little abusive although the two models are the same. A way to pass from the MK–FV equation to the Lotka integral can be found p.148 of [18].

Stable state analysis

Again, one way to analyze the differential system is to look for remarkable solutions and then see if these particular solutions are attractive (i.e. any trajectory converges towards these solutions in long time). The range of remarkable solutions usually studied for differential systems of infinite dimension are the solutions with separable variables, and the renewal equation (I.5) does not deviate from the rule. We are therefore interested in the stationary age profile, anticipated by Sharpe and Lotka, (t, a) 7→ρˆ(a)eλt: the stationary age profile is multiplied by exponential growth. One can obtain the well-known eigenvalue problem by injecting such a solution in Eq.

Table of contents :

Acknowledgements 
I Introduction 
I.1 Modeling cell dynamics
I.1.1 Cell average models (used in Chapter 2 and 3)
a) Exponentially stable state analysis: general theory
b) Stable state analysis: illustration of a failure case of the general theory
I.1.2 Structuring physiological variables at the population level (used in Chapter 3 to 5)
a) The McKendrick–VonFoerster and renewal equations
b) The Lotka integral equation
c) Stable state analysis
I.1.3 Cell-based level (used in Chapter 2)
a) Link with the continuous time Markov chains
b) Simulation of the Poisson process
c) The branching property
I.1.4 Structuring physiological variables at the individual level (used in Chapter 3)
a) The moment-generating functions
b) The Poisson point measure
I.1.5 Data and model matching (used in Chapter 2, 3 and 5)
I.2 Introduction to the early development of the ovarian follicle
I.2.1 The dynamics of follicle growth
a) The ovarian reserve of primordial follicles
b) Folliculogenesis stages
I.2.2 Dataset presentation
a) Follicle types in early development
b) Kinetics information: dataset building of Chapter 3 and Chapter 4
c) Dataset building of Chapter 2
I.3 Contributions and outline of the dissertation
II Modeling the ovarian follicle activation 
II.1 Model definition
II.2 Model analysis
II.2.1 Analysis of the deterministic model
II.2.2 Analysis of the extinction of the precursor cell population
a) Analytical expressions in the linear case
b) Upper bound of the stochastic model
c) Numerical scheme for the mean extinction time and mean number of proliferative cells at the extinction time
II.3 Parameter calibration
II.3.1 Dataset description
II.3.2 Likelihood method
II.3.3 Fitting results
a) Two-event submodels
b) Three-event submodels and complete model
c) Comparison of models
II.3.4 Model prediction
a) Distribution of the initial condition
b) Proliferative cell proportion: reconstruction of time
c) Mean extinction time, mean number of cells at the extinction time and mean number of division events before extinction
d) Biological interpretation
II.4 Conclusion
II.5 Appendix
II.5.1 MLE parameter sets
IIIModeling the compact growth phase 
III.1 Starting point
III.1.1 The [Clément-Michel-Monniaux-Stiehl] (CMMS) model [1]
III.1.2 From the [CMMS] model to our model
III.2 Introduction
III.3 Model description and main results
III.3.1 Model description
III.3.2 Hypotheses
III.3.3 Notation
III.3.4 Main results
a) Eigenproblem approach
b) Renewal equation approach
c) Calibration
III.4 Theoretical proof and illustrations
III.4.1 Eigenproblem
III.4.2 Asymptotic study for the deterministic formalism
III.4.3 Asymptotic study of the martingale problem
III.4.4 Asymptotic study of the renewal equations
III.4.5 Numerical illustration
III.5 Parameter calibration
III.5.1 Structural identifiability
III.5.2 Biological application
a) Biological background
b) Dataset description
c) Parameter estimation
III.6 Conclusion
III.7 Complements on the dataset treatment (unpublished)
IVModeling the compact growth phase: complementary works 
IV.1 Back on the [CMMS] model
IV.2 Mechanical spatial-structured model
a) Free boundary problem: Hele-Shaw model
b) Link with the Keller-Segel model: from incompressible to compressible model
V Inverse problem for a structured cell population dynamics model 
V.1 Model and discretized solutions
V.1.1 Formal solutions of the direct problem
V.1.2 Numerical scheme to simulate the direct problem
V.2 Analysis of the inverse problem (IP) for continuous initial conditions
V.2.1 Single layer case
a) Well-posedness of the inverse problem (IP)
b) Numerical procedure for non-synchronized populations
V.2.2 Multi-layer case
V.3 Analysis of the inverse problem (IP) for Dirac measure initial conditions
V.4 Conclusion
VI Discussion and perspectives 
A Appendix
A.1 Additionnal materials of Chapter 3
A.1.1 Supplemental proofs: deterministic model
A.1.2 Supplemental proofs: Stochastic model
A.1.3 Supplemental proofs: Moment study
a) Stochastic simulation procedures
b) Deterministic simulation protocol
A.1.4 Construction of Figure III.5
A.1.5 Parameter estimation procedure
A.2 Additionnal materials of Chapter 5
A.2.1 Numerical scheme
A.2.2 Complement of proofs
A.2.3 Complementary illustrations of Algorithm
Bibliographie

GET THE COMPLETE PROJECT