## Abstract

The paper presents the classical age-dependent approach of the morphogenesis in the framework of the von Foerster equation, in which we introduce a new constraint and study a new feature: (i) the new constraint concerns cell proliferation along the contour lines of the cell density, depending on the local curvature such as it favours the amplification of the concavities (like in the gastrulation process) and (ii) the new feature consists of considering, on the cell density surface, a remarkable line (the null mean Gaussian curvature line), on which the normal diffusion vanishes, favouring local coexistence of diffusing morphogens, metabolites or cells, and hence the auto-assemblages of these entities. Two applications to biological multi-agents systems are studied, gastrulation and feather morphogenesis.

## 1. Introduction

Morphogenesis processes involve many mechanisms such as proliferation, differentiation, migration, etc. Each of these can be studied and modelled separately, but explaining the growth of an organ or a tissue requires a global approach, implying a common mathematical framework and numerical treatments made of both discrete iterations (in an event or multi-agent form) and continuous equations (in ordinary or partial differential form). We attempt, in this paper, to offer a common formalism for cell growth and tissue morphogenesis in a hybrid approach by showing the correspondence, at a global level (e.g. concerning global parameters such as growth or synchrony indices), between discrete and continuous approaches. The Malthusian growth parameter is the dominant eigenvalue of the proliferation operator, and we will define a new synchronizability index, the cell entropy, measuring the kinetic heterogeneity of the cell cycle in the context of an age-dependent proliferation. Complete information about demographic parameters and operators can be found in McCallum (2000) and Kot (2001).

The morphogenesis processes give birth to many different forms of tissues that are well adapted to their functions. We can combine these processes to obtain, from a large variety of initial conditions, the observed morphologies as asymptotic behaviours (attractors) of dynamical systems based on four types of endogeneous components, i.e. physical terms such as diffusion and convection, chemical terms such as mass action reaction and ionic attraction or repulsion, biological terms such as proliferation and differentiation and pure information agents such as enzymatic effectors (activator and/or inhibitor) or transcription/translation control factors. This internal dynamics of morphogenetic processes can be driven or entrained by exogeneous forces such as external fields of different nature, e.g. electromagnetic, gravitational, chemotactic, haptotactic, etc. By combining these elements and by using different time scales and spatial confinement constraints (represented by boundary conditions), we can explain most of the forms used by the nature for ensuring the main living functions, e.g. nutrition and thermoregulation, gastrulation and feather morphogenesis are some of the embryonic steps involved in these functions.

In the present paper, we will first deal in §2 with the modelling of the cell cycle dynamics (proliferation, growth and apoptosis) by proposing both a discrete and a continuous approach, which are compared in §3 by showing analogies between explicit formulations of the main thermodynamical parameters, i.e. the Malthusian growth parameter *r* and the cell entropy *H*. Section 4 presents some properties of these parameters, especially the ability of *H*^{−1} to quantify the synchronizability of a cell population, and §5 gives some indications about the estimation of *H*. Section 6 introduces, in cell cycle dynamics, other morphogenetic processes, such as differentiation and migration. The whole morphogenesis system is then a hybrid system, which we describe as far as possible with a common formalism. Section 7 describes the problem of synchronizability in the general context of the biological age. Sections 8 and 9 are respectively devoted to gastrulation and feather morphogenesis applications, for which we need to introduce the curvature-dependent proliferation and the null-diffusion front waves successively.

## 2. Cell cycle dynamics: proliferation, growth and apoptosis

### (a) General approach

The cell cycle describes the successive physiological states a cell crosses from the time of its birth to the instant at which it divides (Baserga 1980). Five successive phases can be distinguished: a resting transient phase G_{0}, the growth phase G_{1}, during which the different cell components grow, the nucleic acids synthesis phase S, the pre-mitotic control phase G_{2} and the mitotic phase M, during which the cell undergoes a change in structure with condensation of the chromosomes, resulting in a discontinuous process in which two daughter cells are produced. The five phases of the cell cycle can be further decomposed into different states, by taking into account critical control points sensitive to the regulation in the frame of genetic networks or to external or internal metabolite, morphogen or antibody fluctuations (e.g. Daszynkiewiez *et al.* 1980; Bertuzzi *et al*. 1981; Cohen *et al.* 2001; Cinquin & Demongeot 2002*a*,*b*, 2005; Aracena *et al.* 2003, 2004*a*–*c*; Basse *et al.* 2003; Demongeot *et al.*2003*a*–*c*; Aracena & Demongeot 2004; Baum *et al.*2004, 2006; Laforge *et al.* 2005; Guardavaccaro & Pagano 2006) accelerating the progress to the mitosis or provoking either cell death (apoptosis) or cell differentiation (a dramatic change of cell lineage).

In the mathematical models proposed for the cell cycle, the fundamental parameters considered are the rates at which the cell crosses the different physiological states of its cycle (figure 1) and the main variable is the state density function. We assume, in discrete models, that the system is described in terms of *n* states, corresponding to well-defined biochemical and metabolic criteria, then the state density of the cell population at time *t* is described by the state vector *u*(*t*)=[*u*_{1}(*t*),…,*u*_{i}(*t*),…,*u*_{n}(*t*)], where *u*_{i}(*t*) is the size of the fraction of the cell population in state *i* at time *t*. In continuous models, the state density is described by a non-normalized distribution function *u*(*x*,*t*), which is interpreted as size or concentration. We will describe successively two models: the first continuous, introduced by McKendrick and von Foerster (McKendrick 1926; von Foerster 1959) in terms of a partial differential equation (PDE) and the other discrete, firstly defined by Hahn (1970) in terms of a matricial equation.

### (b) Continuous approach

The physiological states of a cell can be described as a continuous multi-dimensional state variable, based on the internal concentrations of nucleic acids, expressed genes and proteins resulting from this expression, i.e. of the major elements that characterize the cell functions and change continuously in time (Kimmel 1980; Goldbeter 1991; Tyson & Novak 2001), but for the sake of simplicity, we will use in the following a scalar (one-dimensional) internal cell state *x* (stage of progress in the cell cycle or simply physiological cell age) and the external observation time *t*. We will use the formalism of the state density, described by the distribution function *u*(*x*,*t*), defined from [1,*M*] to *R*_{+}, where *M* is the last mitotic state and *u*(*x*,*t*) is the probability for a cell to be in state *x* at time *t* multiplied by the total population size or concentration *N*(*t*),
2.1Different parameters characterize the transition through the different cell cycle states

the maturity coefficient

*β*(*x*), which describes the maturation velocity, i.e. the local relationship between the physiological age*x*and the chronological observation time scale*t*andthe probability that a cell disappears at age

*x*from the population per unit time. Disappearance has two causes, apoptosis with a frequency*μ*(*x*), and for the cells surviving, post-mitosis differentiation with a frequency (1−*Q*(*x*))*β*(*x*) (including the mitotic abortion rate), where 2*Q*(*x*) (0<*Q*(*x*)≤1) is the mitotic rate at age*x*(i.e. the mean number of daughter cells of the same lineage as the cell mother of age*x*).

The dynamical equation describing the evolution of the distribution *u*(*x*,*t*) is a McKendrick–von Foerster-like PDE
2.2with the boundary conditions
2.3where *f*(*x*) denotes the initial age distribution, in general defined by *f*(1)=1, *f*(*x*)=0, for *x*>1.

### Proposition 2.1.

*The asymptotic rate of increase of the cell population, the Malthusian parameter, is the unique real root r for the equation expressing u*(1,0) = 1
2.4*where l(x) is the survival function (i.e. the probability that a cell born at age* 0 *survives at age x*) *given by* *; u(x,t)* = e^{r(t−h(x))}*l(x) is a solution of equation (*2.2*) and a stationary density is given by w(x)* = e^{−rh(x)}*l(x), where* *is the mean observed time between age* 1 *and age x*.

### Proof.

The uniqueness of *r* comes from the fact that the integral in equation (2.4) is a decreasing function of *r*. It is also a consequence of the continuous Perron–Fröbenius theory (Greiner 1984; Sgibnev 2001; Piazzera & Tonetto 2005; Doumic 2007; Clairambault *et al.* 2008) and *u*(*x*,*t*)=e^{r(t−h(x))}*l*(*x*) is a solution of equation (2.2) because ∂*u*/∂*t*=*ru*(*x*,*t*), ∂*u*/∂*t*=−*r* ∂*h*/∂*x* *u*(*x*,*t*)+ e^{r(t−h(x))}∂*l*/∂*x*, ∂*h*/∂*x*=1/*β*(*x*), ∂*l*/∂*x*=(−*μ*(*x*)/*β*(*x*))*l*(*x*); hence, *β*(*x*)∂*u*/∂*x*+ ∂*u*/∂*t*=−*ru*(*x*,*t*)−*μ*(*x*)*u*(*x*,*t*)+*ru*(*x*,*t*)=−*μ*(*x*)*u*.

Solutions of similar cell division operators can be found in Michel *et al.* (2004) and Michel (2006); d*a*/*β*(*a*) is the mean expectation time of reaching the state *a*+d*a* in a negative binomial Levy process (Kozubowski & Podgorski 2007) of parameter *β*(*a*), hence *h*(*x*) is the mean observed time between age 1 and age *x*. ▪

There exists a fundamental decomposition of the Malthusian parameter *r*: it suffices to develop the formulae below to prove it or to see Demetrius (1983), Demetrius & Demongeot (1984), Demongeot & Demetrius (1989)
2.5where *H* and *Φ* are called, respectively, the cell cycle entropy and potential, and are defined by
2.6and
2.7where *V* (*x*)=*l*(*x*)*m*(*x*) is the number of new cells produced by a cell surviving until the age *x*, i.e. the probability *l*(*x*) for a cell to be still alive at age *x* multiplied by the number *m*(*x*) of new cells produced at age *x*, which is called the fecundity function, and *p*(*x*)=e^{−rh(x)}*V* (*x*) is the probability that the mother of a newborn cell be at age *x*. In the above formulae, the quantity represents the mean cell lifespan of the cell population and is denoted by *T*.

We call *H* the entropy of the cell cycle or cell entropy, and we suggest *H*^{−1} as a measure of synchrony because *H*=0 when *p* is a Dirac distribution (maximal synchrony) and *H* is maximum when *p* is uniform (maximal asynchrony). The significance of decomposition (2.5) and the mathematical basis for proposing *H*^{−1} as a synchronizability measure is discussed in §4.

### c) Discrete approach

In the class of models proposed by Takahashi (1966, 1968), Hahn (1970) and Bertuzzi *et al.* (1981), the physiological state or age of a cell is a discrete variable *i* and the cellular dynamics is parametrized by discrete entities depending on the discrete time *t*. The state density of the cell population at time *t* is described by *u*(*t*)=[*u*_{1}(*t*),…,*u*_{i}(*t*),…,*u*_{n}(*t*)], where *u*_{i}(*t*) is the size of the cell population fraction in state *i* at time *t*. The parameters are the birth and death rates that represent transition between the different states of the cycle. Macroscopic variables, such as growth rate and entropy, can be calculated, and experimental data can be expressed in terms of the transition rates. Consequently, discrete models and their macroscopic variables constitute the basic formalism that underlies empirical evaluation of cell cycle, in which age *i* and time *t* belong, respectively, to {1,…,*n*} and {0,…,*T*}. We assume that there are *M* states in the cell cycle between birth and mitosis, and that each phase *P*_{i} of the cell cycle contains *nP*_{i}−*nP*_{i−1} internal states (with *nP*_{−1}=0, *nP*_{0}=1, *nP*_{1}=*nG*_{1}, *nP*_{2}=*nS*, *nP*_{3}=*nG*_{2}, *nP*_{4}=*M*). The total population size at time *t*, expressed in terms of the age distribution, is
2.8The dynamical system describing the change in age distribution of the cell is given by
2.9The elements in the matrix *A* are assumed to be constant in time (Demongeot in press); they represent the rates at which cells cross the ages of the cell cycle, which can be estimated from experimental data (De Boer & Perelson 2005). *α*_{i} is the probability to stay at state *i*, *β*_{i} is the probability to go to the next state (*i*+1) and *γ*_{i} to the next further state (*i*+2), in any lapse of observation time [*t*, *t*+1]. The general structure of *A* is
2.10where 0<*Q*≤1, *α*_{i}+*β*_{i}+*γ*_{i}=1−*μ*_{i}≤1, ∀*i*=1,…,*M*.

The special case, where *α*_{i}=*α*, *β*_{i}=*β* and *γ*_{i}=*γ* for any *i*, corresponds to a model in which cells remain in the same state or cross the different phases of the cycle at the same rates (Hahn 1970). *A* is then irreducible and primitive, hence the asymptotic properties of the dynamical system (2.9) can be completely characterized by invoking the discrete Perron–Fröbenius theorem: the asymptotic growth rate of the cell population is , where *λ* is the dominant eigenvalue *A*. The stable age distribution is the right-normalized (i.e. whose sum of components equals 1) eigenvector *w* associated with *λ*, and there exists a decomposition for the population growth rate analogous to equation (2.6).

Let us consider now the probability matrix *P* defined by *p*_{ij}=*a*_{ij}*w*_{j}/*λw*_{i}, which gives the probability to observe backward the state *j* after the state *i*. We have, because of the definition of the eigenvector *w*,
2.11

Let us denote by *π*=(*π*_{i})_{i=1,…,M} the stationary distribution of *P* (i.e. the normalized left eigenvector associated with the eigenvalue 1). Then, the Malthusian parameter *r* admits the decomposition between discrete cell entropy and potential
2.12where
2.13and
2.14*H* and *Φ* are the analogues of the corresponding functions in equation (2.6) and *H*^{−1} can be considered as a measure of the cell cycle synchronizability for the discrete model, as for the continuous one. The more precise significance of this decomposition and its mathematical basis will be discussed in §4. Before we proceed with this discussion, we give a simple model for which the function *H* can be explicitly calculated. If we consider the Hahn (1970) model with the condition *α*_{i}=*α*, *β*_{i}=*β* and *γ*_{i}=*γ* for all *i*, then the matrix *A* is quasi-circulant (Demongeot *et al.* 1995), and we have
2.15where *S*=(1−(2*Q*)^{−1/M})/(1−(2*Q*)^{−1}).

Then, the matrix *P* is given by
2.16and the invariant probability distribution is given by *π*=(*π*_{i})_{i=1,…,M}=(1/*M*,…,1/*M*).

Then, we can deduce the expressions of the cell entropy and potential
2.17and
2.18We will suppose in the following that *γ* = 0. If not, but if *α*(*x*)=*α*, *β*(*x*)=*β* and *γ*(*x*)=*γ* for any age *x*, we consider the second-order PDE
2.19awhere *γ* represents the maturation acceleration (Bertuzzi *et al.* 1981; Demongeot *et al.* 1995). The main interest of the continuous formulation (2.19a) is the possibility to add a diffusion term, when *u*(*x*,*t*,*s*) depends on age *x*, time *t* and space *s*
which reduces when *β*′(*x*)+2*γ*′(*x*)−*γ*′′(*x*)=0 to
2.19b

By introducing in equation (2.19b) the notion of the bio-demographic Dalembertian , which takes into account the diffusion both in age and in space, we obtain the most general operator including accelerated ageing and population motion.

## 3. Comparison between discrete and continuous approaches

Table 1 summarizes the correspondence between discrete and continuous approaches, by comparing continuous and discrete formulae when the Leslie matrix coefficients *α* and *β* do not depend on the age. The survival function *l*(*x*) is expressed, if *δ* is small enough (*δ* equals 1/*m*, where *m* is a sufficiently large integer) and if *μ*=1−*α*−*β*, as
3.1The last term in equation (3.1) corresponds to the discrete survival function. If *μ* and *β* depend on the state *x*, we can show, in the same way, that the survival function is
3.2if we denote *μ*(*k*)+*β*(*k*)=(1−*α*_{j})^{1/δ} and .

If we replace the Malthusian parameter *r* by its discrete equivalent , then the discrete equivalents of *h*(*x*) and of the expression e^{−rh(x)}*l*(*x*) become, as above,
3.3and
3.4because , if *λ*=*α*+*β*(2*Q*)^{1/n} is sufficiently near 1.

Let us consider now the fecundity function *m*(*x*). Its discrete equivalent is given by
and the discrete equivalent of *w*(*x*) is the probability distribution
3.5where *S*=(1−(2*Q*)^{−1/M})/(1−(2*Q*)^{−1}). It is easy to check that *w* is the discrete stationary density for the dynamics (2.9) found in equation (2.15). We have *p*(*x*)=e^{−rh(x)}*l*(*x*)*m*(*x*), hence
3.6where *k*=(2*Q*)^{1/M}/*λ*.

Finally, the correspondence between discrete and continuous formulae of the cell entropy and potential is ensured by
3.7Then, if *μ*=1−*α*−*β*≈0 and *α*≈0,
3.8where the discrete equivalent of the temperature is given by the mean cell lifespan
3.9We remark that *T*_{d} is approximately equal, if *α* is small, to *M*/*β*, which is the expectation of the waiting time before reaching *M* successes in a Bernoulli game of probability *β* to have a success (i.e. to not remain in the same state of the cell cycle). We eventually have the correspondence
3.10When *μ* and *α* are small, the above formulae show that continuous and discrete approaches are completely similar and the equation of the discrete dynamics (2.9) for the Hahn model can be considered as a difference equation associated with the von Foerster PDE (2.3).

## 4. Properties of the thermodynamical parameters *r* and *H*

We will perform now a complete identification between the cell cycle modelling and the classical statistical mechanics approach, by defining first the random process of successive mother ages within random sequences called genealogies. This process is a Bernoulli backward shift (i.e. generated by a translation on a set of random sequences), dual of the forward deterministic dynamical systems ruling the cell ages. Then we will show that *r* characterizes the forward stability of the stationary age distribution (for the Euclidean distance) and *H* both the backward stability of the Bernoulli shift invariant measure *μ* (for the Kullback distance) and its variability, i.e. its faculty to give asynchronized generations.

### (a) The space of genealogies and the reproduction age shift

A precise characterization of these two stability properties involves the definition of a (backward) genealogy *x*. Genealogy is the term used to describe one temporal sequence of states, such as those crossed by a cell during its life cycle, enabling us to explain the nature of the decomposition given by equations (2.5) and (2.12). We first introduce the concept and indicate how *H* and *r* can be derived. We then indicate how the notions of complexity and convergence can be precisely described in terms of this concept: a backward genealogy represents the sequence of states that a cell crossed before reaching its present state. A forward genealogy describes the sequence of states its descendants will enter as the population evolves. Formally, a genealogy is a sequence
4.1where the *x*_{i}’s correspond to the cell cycle states and *i* is the time.

Let Ω denote the set of all genealogies and let *μ* denote a probability measure on Ω that is invariant under the backward translation operator

Let *H*_{μ} denote the Kolmogorov–Sinai entropy of the dynamical system (). The decomposition (2.12) is derived from the fact that, if *λ* denotes the dominant eigenvalue of the matrix *A* given in equation (2.10), log *λ* satisfies a variational principle described on the space of genealogies (Demetrius 1983). We have
4.2where the potential *Φ*(*x*) is equal to , with *a*_{ij} being the general coefficient of *A*. The variational principle defined by equation (4.2) is formally identical to the principle of minimization of the free energy in statistical mechanics. There is a unique measure *μ** for which this maximum is attained; *μ** is an equilibrium measure and we have
4.3and
4.4where *P*=(*p*_{ij})_{i,j=1,…,M} is the transition matrix of the Markov chain associated with , *π* denotes its invariant probability and *w* is the eigenvector associated with the dominant eigenvalue *λ* of *A*.

Let us now consider the Markov subchain whose states are those observed just before a state 1 in a genealogy; this subchain, called reproductive, is indexed backward by the times 0,−1,…,−*k*,… . The associated Markov transition matrix *R* gives the probability *r*_{ij} to pass backward from an age of reproduction *j* at a certain generation to an age of reproduction *i* at the previous generation, and we have the following proposition.

### Proposition 4.1.

*The reproductive Markov subchain associated with the Hahn model is a Bernoulli shift whose probability transition matrix is made of the succession of identical lines defined by r*_{ij} = a_{1j}*w*_{j}/λ*w*_{1}*. The reproductive entropy of this subchain is equal to H*_{R}*=H*_{μ*}*.*

### Proof.

If *X* denotes the number of states observed backward between states *j*, the first, and *i*, the next after state 1 (*j* and *i* are, respectively, the age at which a mother and a grandmother have got their offspring), we have

Because the lines of *R* are identical, the corresponding Markov chain is a Bernoulli shift whose entropy is the same as those of the Markov chain generating the genealogies, and the corresponding invariant probability distribution is given by *ρ*=(*ρ*_{i})_{i=1,…,M}=(*α*/*λ*,0,…,0,(2*Q*)^{2/M}*γ*/*λ*,(2*Q*)^{1/M}*β*/*λ*), and we have . ▪

### Remark.

If α = γ = 0, the Hahn model becomes a Leslie–Lotka model, and we get the same result; the Leslie–Lotka model also has a Bernoulli shift for the process of its reproduction ages, but some differences are apparent. The analogue of the Hahn matrix *A* is the Leslie–Lotka matrix *L* defined by
where at age *i*, *m _{i}* denotes the fecundity coefficient and

*β*≤ 1 the survival coefficient. Then, by using the same reasoning as above, we have

_{i}In this case, *H*_{μ*} = *π*_{1}*H*_{R} = *H*_{R}/*B*, where *B* = S_{i}*ip*_{i} is the mean biological age of the cell population (useful in modelling the ageing (Demongeot in press)), then *H*_{μ*}<*H*_{R}, which corresponds to the fact that, in the sequences of generations, only the part corresponding to the reproduction states is random (the other part is made of the strict succession of ages).

### (b) Analytical properties of *r*

The mode of convergence in the space of genealogies reflects the actual dynamics that occurs in the cell cycle. It is quite distinct from the mode of convergence based on the Euclidean distance. Thus, if *y*(*t*) denotes an arbitrary age distribution and *w* the stable age distribution (observed in the exponential growth phase), we have
4.5where *λ*′ is the first subdominant eigenvalue of *A* and the norm is Euclidean (Chiorino *et al.* 2001).

### (c) Analytical properties of *H*

The analytical properties of *H*, namely its description in terms of variability of the cell cycle and its characterization as a measure of the rate of decay of synchrony, can be made explicit by considering the genealogies. A case of finite-range phase transition for *P* (e.g. when *α*=*γ*=0, *β*=*Q*=1) permits to exhibit an asymptotic *M*-cycle of densities *π*^{(k)}, and to draw the landscape of confiners and stochastic isochrons corresponding to the successive measures *π*^{(k)} of this cycle, in order to find the stimuli capable of causing a lasting synchrony, if the entropies of the measures are small (Demongeot *et al.* 1995). But, it is rarely the case because, at cell division, daughter cells receive identical genetic material. However, the cytoplasmic components are, in general, unequally distributed: daughter cells receive different amounts of mitochondria, proteins and RNA. As documented by several authors (Prescott 1976; Baserga 1980), the variability in cytoplasmic components generates a heterogeneity in the rate at which the daughter cells cross the different phases of the cell cycle. Individual variations in mitotic rates of a genetically homogeneous population of cells thus constitute an intrinsic property of cellular cultures and underly the rapid decay in synchrony characterizing bacterial populations in which synchrony is initially induced by physicochemical means. The problem of characterizing this property of asynchrony has been addressed by several authors (Zeuthen 1958, 1964; Rueckert & Mueller 1960; Blumenthal & Zahler 1962; Burns 1962, 1964; Engelberg 1964; Barnett-Hall & Waugh 1967; Hahn 1970) and several statistical indices have been proposed. These indices have been expressed in terms of the stationary distribution of cells and are not explicitly derived from the dynamical process that underlies kinetic heterogeneity. A good asynchrony index must have an important feature: it must precisely describe the variations in the rates at which the cells cross the different physiological states of the cell cycle and also characterize the rate of decay of synchrony when the stationary distribution is perturbed. The index we propose, *H*^{−1}, we call it the entropy index because it is derived from the thermodynamic formalism above, provides a parameter that can accommodate and classify the experimental data characterizing the heterogeneity of cellular cultures. The significance of *H*^{−1} as an synchrony index is derived from two properties.

#### (i) H determines the rate of convergence of the cell population to its stable stationary state

The parameter *H* describes the stability of the cell cycle in that it determines the rate of convergence to its stable stationary state. Let *μ** denote the equilibrium invariant probability of the backward Markov transition on the space of cell cycle states [1,…,*M*], and *μ*_{t} the probability measure on this space generated after a time *t* from an initial probability *μ*_{0} by the cell dynamics. It follows from results due to Goldstein (1981) that
4.6for *t* sufficiently large, where the norm refers to the Kullback distance between probability measures. Hence, when *μ*_{0} is a Dirac and *H*=0, the rate of convergence to zero of the majorant *K*e^{−Ht} of the distance of *μ*_{t} to the equilibrium is zero, ensuring the persistence of the initial synchrony. When *μ*_{0} is a Dirac and *H*≫0, the rate of convergence to the equilibrium is of the order of magnitude *O*(e^{−Ht}), provoking a rapid loss of synchrony.

In the case where *β* is small, *H* is related to *λ*′, the first subdominant eigenvalue of the Hahn matrix *A* (used in equation (4.6)). Let us consider the case where *M*=3, we have

Then, *λ*′=*α*+(−1/2+*i*√3/2)(2*Q*)^{1/3}*β* and because *x*=(2*Q*)^{1/3}*β*/*α* is small, we get

In general, *H* is implicitly linked to the whole subspectrum of *A* (Tuljapurkar 1982, 2004).

#### (ii) H describes the complexity of the cell cycle

The parameter *H* also characterizes the complexity of the cell cycle in that it represents the variations in rates at which the cells cross the cycle. First, let us notice that, if *γ*=0, the variance *V* (*W*) of the waiting time *W* before reaching *n* successes in a Bernoulli game of probability (2*Q*)^{1/M}*β*/*λ* to have a success (i.e. to not remain in the same last mitotic state) is equal to *nαλ*/[(2*Q*)^{1/M}*β*]^{2} and quantifies the variability in the succession of states. Another quantity depends on this variability, the mean doubling time *D*, defined as the time required for the cell population *N*(*t*) to double, when it grows with its exponential dynamics,
which implies .

By using the variable *x*=(2*Q*)^{1/M}*β*/*α*, we have *α*/*λ*=1/(1+*x*) and (2*Q*)^{1/M}*β*/*λ*=*x*/(1+*x*). Hence, we have
If *α* is small, *β*≈1 and *Q*=1,*x* is large and we have *H*≈1/*x*, *D*≈*M* and *V* (*W*)≈*M*/*x*≈*DH*. For a given mean doubling time, *DH* varies in this case in the same sense as *V* (*W*). If *α*≈1, *β* is small and if *Q*=1, *x* is small and *H*≈*x*/(1+*x*), *D*≈log 2(1+*x*)/*x* and *V* (*W*)≈*M*/*x*^{2}≈*D*/*H* log 2. For a given mean doubling time *DH* varies in this case, in the opposite sense to *V* (*W*).

The cases studied above show that the relationship between *H* and *V* (*W*) is not simple. Hence, to make the idea of cell cycle complexity more precise in the context of the genealogies, let us consider the set of partial genealogies generated until time *t* by a cell in the last mitotic state *M* at time 0. This set of genealogies can be ordered in terms of their probabilities, and genealogies can be selected from this set starting from the one with the highest probability until the sum of the probabilities of the genealogies selected exceeds a given positive number (1−*ε*). Let denote the total number of partial genealogies selected: we call it the effective genealogy size until time *t* of observation. The asymptotic behaviour of characterizes, in a precise sense, a measure of complexity of the cell cycle.

First, we consider the case where all cells in the cycle pass through successive physiological states at the same rate, without remaining in a given state (figure 1). This cell dynamics describes populations in which synchrony persists when all cells start at the same physiological state: the effective genealogies set is made of a single element and .

Consider the case where the cells pass through the different states at variable rates. This cell dynamics describes populations in which there is a decay of synchrony. In this case, the number for each

*t*, will depend on the variations in the doubling time of the different cells in the population. The asymptotic behaviour of the number ) can be described in terms of entropy. We have (Demetrius 1983) 4.7

Hence, the parameter *H* measures, in a precise sense, the complexity of the cell cycle. The case *H*=0 corresponds to synchronous populations if cells have the same initial physiological state.

The measures of synchrony have, in general, a minimal value in the domain of exponential growth considered to be the situation of maximal asynchrony (Barnett-Hall & Waugh 1967). The main indices of synchrony are as follows.

Scherbaum’s index (Scherbaum 1957; Williams & Scherbaum 1959; Scherbaum 1960; Chung

*et al.*1973), equal to the coefficient of variability of the empirical growth rate*r*defined by , in which*R*is the fraction of cells showing mitosis,*T*_{M}=*M*−*nG*_{2}is the time required for mitosis and*G*is the time from formation of a new cell by division to its termination by division (*G*can be called the interdivision time, mean generation time, mean cell life time or mean cell lifespan).Rueckert–Mueller–Burns’s index, equal to (

*G*/2 −*D*_{1/2})/(*T*_{M}/2) (Rueckert & Mueller 1960; Burns 1962), where*D*_{1/2}is the median of the distribution of the generation time.Blumenthal–Zahler’s index, equal to

*N*(*G*,*t*)/*N*(0,*t*)−2^{G/D}(Blumenthal & Zahler 1962), where*D*is the cell doubling time.Engelberg’s index, defined as the ratio

*A*_{1}/*A*_{2}, where*A*_{1}is the area under the curve of dlog(*N*(*t*))/d*t*when it is strictly over*r*and*A*_{2}is the total area under this curve during the cell doubling time*D*(Engelberg 1964; Hahn 1970).Zeuthen’s index of synchrony (Zeuthen 1958, 1964), equal to where

*G*is the mean cell life time and*N*(*t*), the size of the cell population at time*t*, tends to zero, when the cell culture reaches its asymptotic exponential growth as*t*tends to infinity. Indeed when the cell population is in its exponential growth regime, the mean relative velocity dlog*N*(*t*)/d*t*tends to*r*, and*Z*(*t*) tends to 0, when*t*increases. It is the case in the following example (Zeuthen 1964), with a relative velocity in the cell cycle defined byThis case occurs when and

*N*(*t*)=*N*(0),*r*=0,*Z*(*t*) tends to 0 when*t*increases and the population is not synchronizable.

Another example of asynchronous population, but not synchronizable, is given by Hahn’s model in which all coefficients *p*_{ij} of the matrix *P* have the same value 1/*M*; in this case, , *r*=0 and *Z*(*t*) tends to 0, when *t* increases. Contrarily, if we have in Hahn’s model ∀*i*=1,…,*M*, *β*_{i}=1, then *H*=0, , *G*=*M* and *Z*(*t*) is indefinite, the asynchronous character is only due to uniform initial conditions ∀*i*=1,…,*M*, *u*_{i}(0)=*k*≥1, and the population is synchronizable.

To conclude, indices of synchrony, such as Zeuthen’s index above are, for a given culture, a relative indicator of the proximity to its own exponential growth. *H*^{−1} (used, for example, in Park *et al.* 2006) is an absolute and non-relative index of synchronizability in the sense that we can compare this index for cultures having the same growth rate *r*. More precisely, if *H*=0, for any growth rate *r*, the cells will remain synchronized if they have all the same initial physiological state because they will have the same genealogies, i.e. the same succession of observed states at regular observation times (it is the case in the last example above). If the cells do not have the same initial state, the distribution on the possible initial ‘phases’ causes a desynchronization in the culture, but perturbation techniques (such as temperature shock) or external forces (such as alimentation or day/night alternance) can restore the synchronization, if *H* is small. In this sense, we can say that *H*^{−1} is an absolute degree of synchronizability.

Finally, let us consider the mitotic index *I*(*t*), defined as the ratio between the number of mitotic cells (whose states are between *nG*_{2}+1 and *M*) and the total number *N*(*t*) of cells at time *t*. If *Q*=1 and if we neglect the rate of apoptosis, we can estimate *I*(*t*) in the exponential growth phase (Jagers 1975) by
where *υ* is the fraction of the mean cell life time dedicated to mitotic states. This formula shows that, in identical conditions for *Φ*, *υ* and *T*_{d}, the mitotic index varies as *H*. It is the case when we consider a process of cell populations mixing without interactions, where the entropy increases (see also Michel *et al.* (2004) for the entropy dissipation). More precisely, if we select in the *i*th subpopulation, the genealogies whose probability exceeds a given positive number (1−*ε*), let denote the total number of such genealogies (we have called the *i*th effective genealogies size) observed until time *t* in the subpopulation *i*. Then, we have the following proposition.

### Proposition 4.2.

*Let us consider m cell populations, each of them ruled by a Hahn matrix A _{i}, i* = 1

*,…,n, which are supposed to be defined by the parameters α*= 1 −

_{i},β_{i}*α*

_{i}and Q_{i}. Then, if the mean entropy H(t) is defined by , then*increases:*

*, and if the A*

_{i}

*values have the same dominant eigenvalue λ, the limit for the growth of H is given by*4.8

### Proof.

We can deduce from the result (4.7) that and hence that
and then , where the variance is taken for the probability distribution . The inequation (4.8) is then obtained by applying Jensen’s inequality to the concave function *f* on [0,1]
where , with *S*_{i}=(1−(2*Q*_{i})^{−1/M})/(1−(2*Q*_{i})^{−1}). ▪

The last good property of *H* (and hence of *H*^{−1}) is its robustness; the following result shows that slight perturbations on the Hahn matrix *A* induce small changes in the value of *H*.

### Proposition 4.3.

*If we perturb the Hahn matrix by changing in its coefficients k times the parameters α and β* = 1 − *α, in choosing at the states j*_{1}*,…,j _{k}* ≠ 1,

*M*: ∀

_{s}= 1,…,

*k, αj*=

_{s}*α*+

*εj*= 1 −

_{s}and βj_{s}*αj*

_{s}where εj_{s}is small, then*where*

*is the invariant probability measure of P**

*perturbed from π*=(

*π*

_{i})_{i}_{= 1,…,n}= (1/

*M*,…,1/

*M) and defined by*

*and where H*

_{0}

*(respectively, H*

_{i}

*) is the entropy before (respectively, after) the perturbation α*

_{i}

*given by*

### Proof.

Let us consider *A**, the perturbed Hahn matrix, defined by , except for *i*=*j*=*j*_{s}, *s*=1,…,*k*, where and . Then, by calculating the determinant of *A**−*λ***I*, we can show that and *λ**=*λ*+*η*/[(2*Q*)^{(M−1})^{/M}*β*^{M−1}*R*], with and (*ξ*_{i}−1), where *ξ*_{i} is the *i*th *M*-root of the unity difference of 1. We finish by calculating the matrix *P** and its invariant measure *π**. ▪

## 5. Estimation of the index of synchronizability *H*^{−1}

For estimating the index of synchronizability *H*^{−1}, we choose the following procedure: (i) to estimate the parameters of the Hahn matrix *A* by remarking that *α*, *β* and *γ* can be replaced by empirical frequencies and *Q* can be estimated by the empirical fraction of non-abortive mitoses (cf. De Boer & Perelson 2005) and then (ii) to calculate *H* from its explicit formula, and finally to calculate *H*^{−1}. If we observe only some synchronized genealogies during a long time, we assume that *A* is constant over this time and we identify its coefficients as the transition frequencies of the cells from the same lineage between the successive states of their cell cycle, e.g. by estimating the parameter in a one-dimensional exponential statistical structure with boundary conditions (Demongeot 1981). If we observe many desynchronized cells during a short lapse of time by using specific imaging techniques of segmentation (Brunie *et al.* 1995; Cottet *et al.* 1997; Ronot *et al*. 2000; Aracena *et al.* 2004*c*), we suppose that all cells have the same transition frequencies and we estimate the parameter *a*_{ij} from the empirical conditional frequency of going from *j* to *i*. In any case, a longitudinal or vertical sampling allows the estimation of the coefficients of *A* and then the calculation of *H*^{−1}.

## 6. Relationships between migration, differentiation and proliferation

We summarize in equation (6.1) the relationships between the different morphogenetic processes, i.e. proliferation, migration and differentiation. The influence of external fields (thermic, electromagnetic, chemotactic and gravitational) allows either a positive collaboration between these morphogenetic processes facilitating to reach the asymptotic behaviour of the growth process (converging to the morphogenetic ‘attractor’) or a global inhibitory effect preventing the appearance of the normal phenotype of the form. For example, the thermal diffusion in a morphogenetic dynamics can cause, if it is too important, the complete cancelling of the dissipative structures created by highly nonlinear endogeneous kinetic terms (Demongeot & Laurent 1983) or by the interaction between internal elements and external fields (Castets *et al.* 1990). The action of environmental fluctuations can cause changes in the morphogenetic process leading to the occurrence of local singularities: a convex parabolic singularity on the frontier of a form can become concave, provoking a local increase of proliferation if this concavity is turned toward the nutritive compartment, in order to respect the invariance of the ratio cell perimeter/cell volume (Thom 1972; Forest *et al.*2004, 2006*a*,*b*; Forest & Demongeot 2006, 2008), as observed in pathologic tree growth and in gastrulation or neurulation morphogenesis (Ermentrout 1991; Spirov 1993; Forest *et al.*2004, 2006*a*–*c*; Jiang *et al.* 2004). More generally, if *u*(*s*,*a*,*t*) denotes the density of cells at age *a*, time *t* and location *s*, *C*_{u}(*s*,*a*,*t*) denotes the mean Gaussian curvature of the surface *u* at (*s*,*a*,*t*) and ∇*ρ* (respectively, ∇*c*) denotes the gradient of cell adhesion sites (respectively, chemoattractant), the global morphogenetic equation is
6.1The mutual influence of proliferation, migration and differentiation processes (table 2) can be seen through changes in parameter values of marginal dynamics in the global equation (6.1) above under the influence of slow (with respect to the morphogenesis velocity) parameter dynamics. Presented here as thermodynamic descriptors of the cell cycle dynamics, the growth rate *r* and the entropic parameter *H* can also be used for describing the dynamics of self-organizing biological subsystems or of simplified biochemical systems *in vitro*, and also for modelling the dynamics of mesoscopic scale phenomena, such as tissue morphogenesis or self-ordering dynamics of social systems (Blaschko 1901; Borges 1984; Bolognia *et al.* 1994). Then, the state vector *u* is easily assimilated to the concentration vector of chemicals or of communicating agents in populations of molecules in a solution (Ronot *et al.* 2000; Namy *et al*. 2004), cells in a tissue (Murray 1981, 1989; Tranquillo & Murray 1992; Thüler 2002; Zeng *et al.* 2004; Tracqui 2006), tissues in an organism (Kraissl 1951; Ridge & Wright 1996; Laforge *et al*. 2005) or organisms in a society (Bieniawski *et al*. 2004; Kroo 2004), as proposed in §7.

## 7. Synchrony and synchronizability in biological multi-agent systems

The measure of the synchronizability index *H*^{−1} we described above for the cell cycle can be extended to any biological system at different scales (and more generally to whatever complex system where processes or elements are shown to be coupled and are able to synchronize their activities). It informs us about the intensity of the relationships that unify different biological elements in space or time by giving the level of their synchronizability. Moreover, it gives information on the ability of those elements (or processes) to sense the synchronizing signals (coupling signals) from the other elements (their synchronizability), as the different elements of a cell do during the cell cycle.

During embryogenesis or the growth of an organism (e.g. plant growth), populations of cells synchronize their differentiation and their proliferation within specific tissues. Moreover, synchronization also occurs between differentiated tissues for their setting up in the organism, i.e. their respective positioning and the emergence of specific communication pathways between them. This setting up is under the control of one or several chemical or mechanical signals in the direct neighbourhood of the tissues emitting these signals, e.g. chordal induction causes the differentiation of several tissues, such as the neural tube that develops close to the chorda, usually in the form of gradients, e.g. gradients of hormones such as auxin causing plant cells differentiation (Thellier *et al.* 2004). The chorda also causes lower differentiation of other tissues farther from the neural tube. Such systems can be modelled by sets of differential equations, for example, by following variables such as concentrations of cells through the spatial and/or temporal dimensions. In such models, cells ‘communicate’ (in other terms they adapt their current states in the function of external chemical signals) by way of diffusing chemicals, or they migrate due to chemotaxis behaviour. In those systems, the principal factor that determines the manner they will order is the ratio between the term of reaction (*R*) and the term of active and/or diffusive transport (*D*). In classical reaction–diffusion systems, this is expressed as the ratio *R*/*D* (reaction/diffusion) (Kolmogorov *et al.* 1937; Rashevsky 1938; Turing 1952; Castets *et al.* 1990; Harrison 1993). Those systems are good examples where the synchronizability measure *H*^{−1} would be very useful. Organisms, organs, cells and finally cell-compartmented areas, organites or sub-systems such as regulatory and metabolic networks, are able to synchronize with other systems of the same type (e.g. cells with other cells within a tissue or an organ) but also with their subsystems (global control of the subsystems by centralization and integration of their activities) and reciprocally with the system that contains them. For example, all the activities (metabolic, genetic, mechanical, etc.) are synchronized within a cell, but the cells are also synchronized within their tissue, and so on over different levels. Moreover, those systems can be considered as individually functioning agents that behave as sources and receptors of information and that can process it. Within a population (of organites, cells, tissues, etc.) evolving in a certain environment, those agents ‘communicate’ and in that way synchronize each other. Many biological systems begin to be studied as populations of interlinked individuals, i.e. as complex or collective systems. This is the case, of course, for societies of biological agents, such as social insects or human societies (Theraulaz *et al.* 2002), but also for populations of molecules or molecular complexes (Vincent & Thellier 1983; Thellier *et al.* 2003) having individual reacting activities and coupled by way of their chemical environment.

A particularly interesting class of collective systems is that which corresponds to the ‘collective trail systems’. Well-known examples are human trail systems (paths in the grass formed by pedestrians (Helbing *et al.* 1997; Goldstone & Roberts 2006)), ant trails (Theraulaz *et al.* 2002; Johnson & Rossi 2006), hydrodynamic trails (e.g. birds in a flight formation (Bieniawski *et al.* 2004; Kroo 2004)), but other systems could correspond to that definition, such as some chemotactic bacterial or cellular systems (e.g. immune cells that follow the gradients of cytokines), coupled transmembrane proteins such as those of the inner mitochondrial membrane (Demongeot *et al.* 2007*a*) or cytoskeleton-based systems (listeria-actin comets, individual or bundles of microtubules or actin filaments, etc. (Glade *et al.* 2002; Glade 2008)). Those systems are composed of numerous agents that are both sources and sensors of information. The agents move in their medium, randomly (e.g. bacteria, ants) or on the contrary as active walkers (e.g. humans), and release along their trajectory a track of various nature (chemical, visual, hydrodynamic). Moreover, they are permanently sensing their surrounding medium and its variations (composition, structure, concentration, etc.). Thereby, the other agents (and the emitting ones) will sense such trails and that will modify their activity. This is a kind of communication pathway, mediated by variations in the environment, between individuals of various nature. The trails are however submitted to an evolution as a function of time. Once produced, if nothing reinforces them, they diffuse (e.g. pheromones) or are progressively degraded and modified by the environment (e.g. when grass germinates, the pedestrian paths progressively disappear). Spatial self-organizing behaviours can emerge from the trail system when self-reinforcing activities occur. This is the case in ant colonies, or in human trail systems where the agents (individual ants or humans) themselves use those pathways and reinforce them. Such trails are useful for the individuals concerned because they help synchronizing activities, yielding optimized behaviours at the level of the collection.

The efficiency of such trail systems depends on several other factors compared with classical reaction–diffusion systems. The shape of the trail depends (i) on the rate of the agent, (ii) on the intensity of its released signal (the source of the trail), (iii) on the persistence ability of the emitted trail (related to degradation or diffusion), and (iv) on its spreading ability (due to transport phenomena). Moreover, it depends on whether the trail can be sensed by the agents or not, as it depends (v) on the sensitivity of the individuals to the ‘trail signal’. This can be reduced to only two parameters; the synchronizability *H*^{−1} is intrinsic to the agents and corresponds to point (v), whereas the points (i)–(iv) are components of a function *E* expressing the efficiency of the trail signal. In these systems, emergence of self-organization behaviours depends on both *H*^{−1} and *E*. The relationships between these two parameters can be viewed as discrete reformulations of the relationships between *R*, *D* and *E* in continuous systems because, here, the agents are discrete elements in the environment and behave as punctual sources and processors of information.

Human or ant-trail systems are very efficient in space because the trails are intense, sharp-shaped along the trajectories of the agents and because the individuals are very sensitive and reactive to those signals. In these systems, self-ordering phenomena are frequent. On the contrary, molecular trails such as bacterium or worth microtubule trails are much more less efficient, as they are less intense (they are composed of small amounts of molecules) because diffusion tends to rapidly homogenize the medium and is high when compared with the size of the trail, and because the individuals are less sensitive to such signals. To summarize, in molecular trails, *R* and *E* are terms weaker and *D* is higher than in human trails. The following example describes this problem well. Because microtubules or actin fibres are dynamic supramolecular assemblies of protein bricks, they are intrinsically trail systems; when disassembling, they release their subunits along their trajectory, and these subunits can be potentially used by other fibres. Nevertheless, although they are structurally made to behave as trail systems, the question of their efficiency to produce order by this mechanism remains uncertain. In fact, it depends on the level of reactivity of the fibres. In parallel to studies of microtubules in realistic conditions (*in vivo* or in cellular extracts), solutions of purified microtubules were studied for their ability to produce self-ordered temporal or spatial behaviours. Those solutions are composed only of tubulin subunits (the bricks of the microtubules) and of guanosine triphosphate (GTP) (an energetic nucleotide necessary for the microtubular reactions) diluted in a buffer containing ions, and particularly magnesium, which has been shown to control the reactivity of microtubules (the more magnesium, the more reactive they are). Three remarkable behaviours were observed in similar conditions of reaction, but by varying the levels of magnesium: (i) at 20 mM magnesium, propagating waves (1 mm large) of concentration of microtubules form, moving at about 1 mm min^{−1} through the sample (Mandelkow *et al*. 1989), (ii) at 10 mM, damped temporal oscillations of the concentration of assembled tubulin (microtubules) were observed over about 10 cycles, but this time over the whole sample (Carlier *et al.* 1987; Pirollet *et al.* 1987), and (iii) at 1 mM magnesium, no periodic behaviour occurred, but spatially self-organized striped morphologies developed progressively during hours (Tabony 1994). Although the very dense and viscous solutions in (iii) were shown to result from coupled growth of stable microtubules and mechanics (Liu *et al.* 2006), and not from the trailing properties of microtubules (Glade 2008; Glade *et al.* 2009), the phenomena observed in (i) and (ii) are due to synchronizing effects between microtubules through the solutions. As shown recently (Glade 2008), the trails produced by individual microtubules are very weak (variations of very few molecules) and made of tubulin-guanosine diphosphate (GDP) (the inactive form of tubulin that can trigger the disassembly of a microtubule, if assembled). Moreover, their spreading around the emitting end of a microtubule is very large (up to 1 *μ*m) compared with the size of the source (20 nm in diameter). Thereby, these ‘trails’ cannot be considered as thin and long tracks, as they are in human or ant systems, but more as classical diffusion fronts spreading from the disassembling ends of microtubules (or as round and symmetrical trails). However, when the reactivity of the microtubules is sufficient and when mechanics are insignificant, the ‘microtubular trails’ play a role in the synchronization of microtubules throughout the sample. Here, the synchronizability of the agents (microtubules) and the intensity of the sources are directly linked to the reactivity of the system. When it is high, disassembling microtubules rapidly produce a local increase of tubulin-GDP that cannot diffuse much during this period. If a node of microtubules starts to disassemble somewhere, it produces a large local increase of tubulin-GDP, thus inducing a local inhibition of the assembly. As in classical reaction–diffusion-based excitable media, the inhibition propagates (here, it induces the disassembly of neighbouring microtubules that, in their turn, cause the formation of an inbibition area), while it decreases in the previously inhibited regions (tubulin-GDP is regenerated into tubulin-GTP and microtubules reassemble (Sept 1999)).

When the ratio of *R*/*D* decreases, *H*^{−1} and/or *E* decreases; for example, by reducing the concentration of magnesium slowing down the reactivity of microtubules, the diffusion has time to homogenize the medium before the microtubules react, so the oscillations only occur at the level of the whole sample and cannot be observed locally. This is the case for the temporal oscillating solutions of microtubules (Carlier *et al.* 1987; Pirollet *et al.* 1987). Finally, one can expect that, in more reactive conditions such as those used in Mandelkow *et al*. (1989), microtubules could form nice trails and self-organize by local communications as ants. This cannot be done experimentally, but simulations can virtually verify such assumptions by allowing dynamic fibrillar supramolecular assemblies to react at non-realistic rates, 100–1000 times faster than real ones. In this case, the resolution of the trails increases and the synchronization too, at lower scales (Glade 2008). This example perfectly illustrates the notion of synchronizability in multi-agent systems, where the measure of *H*^{−1} makes sense.

## 8. Gastrulation

The gastrulation process is a good example of a multi-agent system in which a cell (called a bottle cell) differentiates due to transcription and/or translation factor fluctuations, becoming concave (figure 2), in contrary to the other exothelial cells that remain convex (i.e. having a nutritive interface with the amniotic liquid more important than the basal interface). By using the global equation (6.1) for the cell densities (with big *r*/small *H* for the concave bottle cells and small *r*/big *H* for the convex exothelial cells, realized by imposing in equation (6.1) a proliferation that is proportional to the local concave curvature *C*_{u} and null in the convex case, as in Thom (1972), Forest *et al.* (2004) and Forest & Demongeot (2006)) and a multi-agent system for the cell wall positioning (Forest & Demongeot 2008), we can simulate the first phase (figure 3) of the gastrulation process, which ends by the collapse of the external parts of the invagination.

## 9. Feather morphogenesis

The morphogenesis of feather primordia is an embryonic process, which allows the adult feathers to be correctly located. For example, the feather wheels of peacocks are well-arranged structures that are allowed by this self-organizing embryonic process.

The reaction–diffusion system corresponding to feather morphogenesis (Michon *et al.* 2008) rules three variables, the density *n* of migrant primordial cells and the concentration *u* (respectively, *v*) of an activator (respectively, inhibitor), the bone morphogenetic protein 7 (BMP-7) (respectively, BMP-2), following the equations
9.1with and , and Neumann boundary conditions.

For the sake of simplicity, we will use, in the following, a simplified equation for *n*,
9.2Then, it is possible to derive Turing’s (1952) instability necessary conditions explicitly as follows, where *u*_{s} (respectively, *v*_{s}) denotes the stationary concentration of *u* (respectively, *v*) and *n*_{s}=1:
9.3If *v*≪1 and *n* is near its stationary value, e.g. if *D*_{n}, *χ* and *b* are large, such as the system reaches rapidly its slow (*u*,*v*) manifold, we can decompose the two last equations of equation (9.1) in order to obtain a potential Hamiltonian system: ∂*u*/∂*t*=−∂*P*/∂*u*+∂*H*/∂*v*, ∂*v*/∂*t*=−∂*P*/∂*v*−∂*H*/∂*u*, with *P*=(*k*_{u}*u*^{2}+*k*_{v}*v*^{2})/2 and . Then, *c*_{1} and *c*_{2} (respectively, *k*_{u} and *k*_{v}) can be considered as more frequency (respectively, amplitude) modulating parameters (Demongeot *et al.*2007*b*,*c*; Forest *et al.* 2007; Glade *et al.* 2007), and the synchronizability can be estimated by considering the isochrons landscape of the simplified system (Demongeot & Françoise 2006).

Another very important parameter is the ratio between the diffusion coefficients *D*_{u}/*D*_{v}, which is less than 1, as is usual in lateral inhibition (Demongeot *et al.* 2009): if the ratio is equal to the critical value 0.06, we observe both in experiments (figure 4) and in simulations (figure 5) a spatio-temporal synchrony between the effectors *u* and *v*; in particular, experiments and simulations show a coincidence of their remarkable Gaussian lines, i.e. the projections of the null-curvature lines on the *u* and *v* concentration surfaces, verifying
and

The two remarkable lines for the effectors *u* and *v* coincide for the critical value of *D*_{u}/*D*_{v}=0.06 (figure 6). The two-dimensional projection of these lines forms a front wave moving in the same direction as the fronts of the concentration surface contour lines; hence, if the vertical curvature of the concentration surface at a point of a remarkable line has the same value as the mean Gaussian curvature at this point, then the diffusion vanishes and *u* and *v* are susceptible to form at this location an assemblage similar to the phospho-lipoproteic plasmic membrane or in the inner mitochondrial membrane (Demongeot *et al.* 2007*a*). The coexistence at this same least-diffusion location of migrant cells *n*, as well as morphogens *u* and *v*, permits an anatomic boundary to be built for future feathers, avoiding chemical reactions between these components, which change their physical nature and have general thermic fluctuations (hence no null diffusion). These phenomena are summarized in figure 7, which shows the coincidence (or the spatial synchrony) between the remarkable lines in one and two dimensions, suggesting that this mechanism can be met in several circumstances of formation of an anatomic boundary; for example, in Demongeot *et al.* (2007*a*), a lateral inhibition mechanism is also used to show a spatial synchrony between transmembrane proteins (the ATPase and the translocase), allowing the realization of a variational principle, which maximizes the mitochondrial production of adenosine triphosphate (ATP) and minimizes the mean free path of the adenylates inside the mitochondrion, by favouring the spatial vicinity between ATPase and translocase sites.

## 10. Conclusion

The kinetic heterogeneity that characterizes populations of replicating cells resides in the unequal distribution of metabolic components which occur at cell division. This kinetic heterogeneity is expressed by the differences in the rates at which the cell passes through the different phases of the cell cycle. Here, we have analysed the relations between two ways of modelling these phenomena, one model based on a representation of the physiological states in terms of continuous variables, the other in terms of discrete variables. The description in terms of discrete variables represents an approximation of the cellular dynamics that evidently unfolds in continuous time, but we have shown that the macroscopic variables (namely growth rate *r* and cell entropy *H*) describing the underlying phenomena have the same analytical expression, independent of the discrete or continuous character; hence, the present work provides a mathematical basis for using the discrete description of the changes in states of the cell cycle for explaining and predicting with the same precision as the continuous approach. The principal advance provided by the synchronizability index *H*^{−1} is to give information over time and space on how the different parts that constitute an organized complex system are synchronized between each other and on their ability to be synchronized. In other terms, it measures the sensitivity of the system with respect to the variations of states of its elements, within space and/or time, coupled by diffusion, convection, chemotaxis, haptotaxis or whatever other physicochemical coupling. This information can be simply extracted from such a coupling when described by the general formalism given above for morphogenetic modelling, thus furnishing a better understanding of how individuals, such as cells, work in their globality compared with classical analyses that focus on single elements or single relations between them, or that study the systems only at their stationary states. Populations of molecules or molecular complexes, as well as populations of biological agents such as social insects or human societies, that have individual reacting activities and are coupled by way of their chemical environment constitute a domain of reference for applying the mathematical tools introduced in this article.

## Acknowledgements

We appreciate the support of the EC NoE VPH project. This work has taken great benefit from numerous exciting discussions with L. Forest (accidentally deceased in the Alps in February 2008), whose intellectual creativity, integrity, passion and humanity remain unforgettable for ever.

## Footnotes

One contribution of 17 to a Theme Issue ‘From biological and clinical experiments to mathematical models’.

- © 2009 The Royal Society