A Model for Cell Proliferation in a Developing Organism.

In mathematical biology, there is a great deal of interest in producing continuum models by scaling discrete agent-based models governed by local stochastic rules. We discuss a particular example of this approach: a model for the proliferation of neural crest cells that can help us understand the development of Hirschprung's disease, a potentially-fatal condition in which the enteric nervous system of a new-born child does not extend all the way through the intestine and colon. Our starting point is a discrete-state, continuous-time Markov chain model proposed by Hywood et al. (2013a) for the location of the neural crest cells that make up the enteric nervous system. Hywood et al. (2013a) scaled their model to derive an approximate second order partial differential equation describing how the limiting expected number of neural crest cells evolve in space and time. In contrast, we exploit the relationship between the above-mentioned Markov chain model and the well-known Yule-Furry process to derive the exact form of the scaled version of the process. Furthermore, we provide expressions for other features of the domain agent occupancy process, such as the variance of the marginal occupancy at a particular site, the distribution of the number of agents that are yet to reach a given site and a stochastic description of the process itself.

Introduction

Cell motility and cell proliferation play an indispensable role in collective cell spreading, which is essential to many key biological processes, medicine and pathological mechanisms including foetal tissue growth, tumour growth and wound healing (see e.g., Maini et al. 2004a; Deng et al. 2006; Cai et al. 2007; Friedl and Gilmour 2009), cell mobility in orchestrating morphogenesis during embryonic development (Gilbert 2003; Keller 2005; Binder et al. 2008; Vulin et al. 2014), immune responses Madri and Graesser (2000), tissue domain expansion (Crampin, Hackborn and Maini 2002), cancer Weinberg and Hanahan (2000) and vascular disease Raines (2000).

Cell proliferation and cell division process lead to domain expansion. This can change the cellular spatial distribution and consequently lead to mechanisms of transport for proliferating cells. Understanding this mechanism in controlling collective cell migration can, potentially, enable us to prevent any abnormalities such as developmental anomalies or cancer metastasis. Agent-based cellular models have widely been used to model proliferative tissue growth (Binder and Landman 2009; Hywood et al. 2013b; Ross et al. 2016). The methods that are typically used to develop these models are based on using a master equation or by using the Fokker-Planck partial differential equation approach, Pillay et al. (2017). Applications of using these developed mathematical models have extensively been used for bone and tumour growth (Czarnecki et al. 2014; Monteagudo and Santos 2015), embryonic tissue growth by assuming time evolution of the length of the tissue domain was known (Binder et al. 2008; De Oliveira and Binder 2019) and fungal and yeast colonies (Matsuura 2000; Tronnolone et al. 2017). Therefore, to improve our understanding of proliferation cell processes, mathematical models have been formulated for the behaviour of the cell population as a whole, based on the behaviour of the individual cells within the population. For example, in the model for the growth of the foetal gut proposed by Hywood et al. (2013a), neural crest cells populate the gut during its growth and remain in the gut tissue to form the neurons that construct the enteric nervous system (ENS). Failure of the neural crest cells to invade the gut tissue completely can cause imperfect formation of the ENS, and lead to the existence of Hirschsprung’s disease Newgreen and Young (2002).

According to Grosfeld (2008), the first report of this disease may have occurred back in 1691. However, the disease is named after the Danish pediatrician Hirschsprung (1981), who described it. This disease can result in serious colon complications, such as enterocolitis and toxic megacolon, which can be life threatening, and it is vital that it be diagnosed early so that medical care can be commenced in a timely manner.

In the literature, there are two main approaches to explaining collective cell spreading: (i) discrete models, and (ii) continuum models Murray (2012). The first class, discrete, individual-based models, usually involve discretizing time and/or space. These models capture the stochasticity and nonuniformity observed in experiments Simpson et al. (2007). The benefit of these models is that they can predict features based on information Murray et al. (2011) and fluctuation Simpson et al. (2014) that can be measured in experimental data Simpson et al. (2010).

In parallel to discrete models, much attention has also been given to developing continuum models, illustrating the behaviour of populations of cells using reaction-diffusion equations. Continuum models can be used to clarify the relationships among the model parameters so that they can be investigated numerically Maini et al. (2004b). The interplay between pattern formation and domain growth, studied in Crampin, Gaffney and Maini (2002), and Woolley et al. (2011), makes use of both deterministic and stochastic aspects.

A multiscale conception of the complex processes driving cell migration can be attained by linking these two modelling approaches together in an equivalent framework; this gives intuition into the interplay between the individual level and population level models that enable us to use either modelling scenario to study appropriate properties of the system. Hywood et al. (2013a) developed such a linkage by studying a discrete, stochastic, domain growth model where discrete ‘domain agents’ representing the gut cells proliferate. This model is a one dimensional semi-infinite lattice model for tissue growth with lattice spacing . Let N(t) be the number of domain agents at time t and be the length of the gut at this time. It follows that is the position of the ith lattice site. In this model a proliferation event can occur, which results in the inclusion of an additional domain agent into the lattice. Consequently, the number of agents rises by one with each proliferation. In this process, if the domain agent at site i is selected to proliferate, it moves to site and pushes all domain agents in sites one place higher.

This process is depicted in Fig. 1 which demonstrates a realization of two proliferation events with a subset of blue marked and green unmarked agents. Whenever a marked or an unmarked agent is selected to proliferate, it moves one site to the right and an unmarked domain agent is inserted into the lattice. If a marked domain agent is selected to proliferate, then the addition of the unmarked domain agent will not only extend the lattice but possibly also separate marked domain agents.

Fig. 1

(Color online) A diagrammatic representation of a growing lattice made up of domain agents, some of which are marked (blue [dark gray]). Initially there are eleven domain agents in total. The domain agents occupying sites 4-8 are marked. The remaining are unmarked (colored green). The new unmarked domain agent is added to the lattice, occupying site 8. The arrows show the movements of respective domain agents. In the next step, the unmarked domain agent at site 3 is selected to proliferate. The new unmarked domain agent is added to the lattice, occupying site 3

By letting the spacing shrink to zero, Hywood et al. (2013a) derived a PDE (their Eq. (14)) to describe how the continuous expected occupancy evolves in time. The coefficient of the diffusion term in this PDE is a linear function of , and hence it should be taken to be zero as the lattice spacing shrinks. However, Hywood et al. (2013a) observed a very good fit between simulations and the trajectories of this PDE, when is taken to be small and positive.

Our aim in this work is to analyse the model of Hywood et al. (2013a) in a different way, by observing that the proliferation process of Hywood et al. (2013a) can be interpreted as a superposition of Yule-Furry processes. We are able to derive expressions for the expectation and variance of the occupancy at a given point, the distribution of the number of agents yet to reach a given point and a stochastic description of the process itself.

The Proliferation Process

In this paper we work on the proliferation process introduced by Hywood et al. (2013a). Let be a random vector that models the state of a proliferation process at time t, with N(t) the number of sites at time t and if site i is occupied by a marked domain agent at time t and equal to zero otherwise. Let denote the time at which the k-th proliferation event occurs with and the time between the kth and st events. We assume that the random variable is exponentially distributed with rate , so that the expected value of is , and that the location of the st proliferation event has a discrete uniform distribution on . If the proliferation event occurs at site j then .

The above assumptions imply that the probability that a proliferation event occurs in , when the state is X(t), is and the location of such a proliferation is to the left of site i with probability . Noting that if a proliferation event occurs exactly at site i, it follows that the expected occupancy of site i at time satisfies the equationTaking expectations of both sides, we can derive an equation for ,We can model an initial condition where there is a contiguous segment of lattice sites that are occupied by marked agents by taking .

To derive a PDE, we make the transformation of variables and to C(x, t), divide both sides of Eq. (2) by , and take limits and so thatIf we want to take as an initial condition that on some interval (a, b], and zero otherwise, we need to let both r and s in the above model approach infinity in such a way that and .

Equation (3) is valid for , with the length of the domain at time t in the continuous model. With , observing that, for any t and in the discrete model,it follows that satisfies the differential equationand soWe see that the expected length of the domain in the discrete model grows exponentially. Via the limiting procedure used above, this is also true of the length of the deterministic domain of the PDE (3).

Hywood et al. (2013a) proposed that the expected occupancy C(x, t) is well-modelled by the second-order partial differential equation.Because of the form (5), the first order term on the right hand side of (6) is the same as that in (3). Hywood et al. (2013a) noted that the trajectories of Eq. (6) fitted simulated results well. However it is clear that the second order term on the right hand side should be zero, as . Indeed, the two limits and required to derive (6) cannot both be finite and non-zero. Thus, this equation has to be regarded as an approximation. In Fig. 2, we compare Eq. (6) derived by Hywood et al. (2013a) with simulation results averaged over 1000 of realizations for small values of . We can see that it does not fit well.

Fig. 2

Expected occupancy averaged over 1000 realizations (blue) and Solutions to Eq. (6) (green curve) with , and at times with initial particle mass on the interval [12, 18] and

In Sect. 3, we shall present an alternative view of the proliferation process that exploits the fact that individual particles move according to a Yule-Furry process. In the following Sect. 4, we shall scale this representation to derive the exact form of the continuum limit of the model. We provide the proofs of the main results in Sect. 5.

Yule-Furry Processes and Proliferation Processes

The pure birth process with a Markovian structure introduced by McKendrick (1914), is considered as one of the simplest branching processes. In the case that the birth rate is linear, it is called the pure linear birth, classical Yule or Yule-Furry process see, for example, Yule (1925), Bharucha-Reid (1997), De La Fortelle (2006) and references therein. The Yule-Furry process has been used to model stochastic population growth with a preferential attachment process as a model of speciation in biology. It was developed with the aim to model various stochastic dynamical systems processes in different fields ranging from epidemics, demography and population biology to server queuing modelling Elia and Taricco (1992), mathematical biology Baake and Wakolbinger (2015) and branching processes Kendall (1966). For instance, Romero-Arias et al. (2018) analysed a growth model of an avascular tumour that considers the basic biological principles of proliferation and genetic mutations of the cell. They modelled the cell mutations dynamics by using a Yule-Furry process and by assuming that cell mutation depends only on the previous genetic state.

Mathematically, the Yule or Yule-Furry process is a pure-birth continuous-time Markov chain on the state space with transition rates for some , which is known as the splitting rate, see for example (Taylor 1975, Page 122). It is well-known that if a Yule-Furry process starts in state j at time zero then its position at time t has a negative binomial distribution with parameters j and , that is the probability that it is in state at time t iswith if .

A domain agent in the proliferation process of Hywood et al. (2013a) moves according to a Yule-Furry process. Observing that the number of domain agents at a given position at time t is necessarily 0 or 1, it easily follows thatNow we are interested in the expected total number of particles at position k at time t, therefore, we consider of these processes, , independently evolving on the positive integers, with common splitting rate and with different starting points . Let the number of particles at position k at time t in the jth process be denoted by . Then the expected total number of particles at position k at time t is given byAn important observation, made by Hywood et al. (2013a), is that this expected total number of particles is given by (9) even if the processes are dependent, and that imposing a particular type of dependence results in the discrete space proliferation model described in Sect. 1.

To see this, observe that the position Y(t) at time t of a particle evolving according to a Yule-Furry Process with splitting rate and with satisfies the stochastic equationwhere is the number of points in an inhomogeneous Poisson process with an intensity measure depending on the evolution of Y on the interval [0, t) given byTo simulate a Yule-Furry process according to this viewpoint, we could use the fact that the integrand in (11) is piecewise constant and so, when , we draw an exponential random variable with parameter and increase the state to when this time has elapsed.(Fig. 3)

Fig. 3

Left-hand panel: the integral (11) evaluated with with respect to t. The jumps occur at exponential intervals with increasing rate. Right-hand panel: the integral (11) evaluated with respect to Y. The jumps are triggered by events in the Poisson processes , , and respectively

Another way to simulate the process would be to generate independent homogeneous Poisson processes with rate at each integer point i and move the particle at k one place to the right at any instant when there is an event in one of the Poisson processes with . This is equivalent to thinking of the random intensity measure in (11) as generated by the independent homogeneous Poisson processes . Specifically, by expressing the area represented by the integral in (11) as a sum of the areas of horizontal rectangles, we see thatwith for and for . These two viewpoints are illustrated in Fig. 3 for a case where . In the left hand panel the jumps are events in a Poisson process of increasing intensity. In the right hand panel, they are events in the superposition of an increasing number of Poisson processes.

Now let’s return to the processes driving Eq. (9). In the spirit of (12), we can think of each of these as being generated by a sequence of Poisson processes, the ith process corresponding to the integer point i. If these processes are independent for different values of j, then is the asymmetric random walk process of Hywood et al. (2013a).

However, if the processes are identical so that for all j, then the model is the proliferation process of Hywood et al. (2013a): points in the ith Poisson process correspond to proliferation events at the ith site and move particles from all processes that are to the right of i simultaneously. In this way, the proliferation process can be regarded as a coupled version of the asymmetric random walk. This leads us to the following theorem.

Theorem 3.1

Let be a random vector that models the state of a proliferation process as defined in Sect. 2 at time t, with N(t) the number of sites at time t and . Then where is the hypergeometric function (Abramowitz and Stegun 1965, Definition 15.1.1).

The expected number of points at site k in the proliferation process is given by in (9).

The variance in the number of points at site k in the proliferation process is given by .

For , let , with , so that measures the distance between the jth and st marked domain agents at time t. Then the random variables are independent with common probability mass function given by

For , let be the number of marked domain agents in positions at time t. Then

Proof

See the Appendix.

Asymptotic Properties

With a view to establishing a continuum limit of the model in Sect. 3, we suppose that there are particles initially situated at positionswithin the interval (a, b] with equal spacing . The idea is that there is an initial ’particle mass’ distributed evenly among these N(0) points. The point at position carries a mass , which we can think of as a rectangle of unit height over the interval .

The expected height of the rectangle at position y at time t,where with respect to (9), we have used the correspondence , and . Note that is a piecewise constant function, with jumps at points of the form . Furthermore, k, r and s increase at the same rate as goes to zero, in particular, and .

For and an initial probability mass on the interval [12, 18], the blue lines in the left-hand panel of Fig. 4 depict as a function of y for values of with n varying from 1 to 10. Decreasing corresponds to increasing the height of the curves and decreasing the tail mass. The right-hand panel of Fig. 4 shows as a function of for and (approximately where the peaks occur in the left hand panel of Fig. 5) for values of n up to 100. The blue lines in Fig. 5 provide another illustration of y plotted against with 107 marked agents initially located at sites 12 up to 118, and with n ranging from 1 to 5. We observe that, as decreases, the curves become flatter with decreased tail mass.

Fig. 4

Left-hand panel: values of at with initial particle mass on the interval [12, 18]. Right-hand panel: plotted against n,

Fig. 5

Left-hand panel: Evaluation of at with initial particle mass on the interval [12, 118] . Right-hand panel: plotted against different values of

The limiting profile C(y, t)

Figures 4 and the left hand panel of 5 suggest that, taken as a function of y for a fixed value of t, approaches a limit as . We investigate this behaviour in this section, commencing with an intuitive characterisation using the normal approximation to the binomial distribution. Equation (14) for the expected number of particles at position y at time t when the spacing is has the formwhere is a random variable with a binomial distribution with parameters and . The mean and variance of are, respectively, and .

Employing the normal approximation to approximate when is small, we getwhere Z is a standard normal random variable. The performance of this approximation is illustrated in red in the left hand panel of Fig. 4 and in the left-hand panel of Fig. 5. In our numerical studies, without loss of generality, we assume the value of to be 0.69. We plotted for different values of in the right hand panel of Fig. 5. In Fig. 6 we compare the normal approximation result in Eq. 16 with Eq. 6, derived by Hywood et al. (2013a). We see that Eq. 16 fits well with the simulation results averaged over 1000 of realizations and for , and . Figure 7 illustrates how expression (16) varies with time.

Fig. 6

Expected occupancy averaged over 1000 realizations (blue), Solutions to Eq. (6) (green curve) and Occupation density (normal approximation (red)) with , and at times with initial particle mass on the interval [12, 18] and

Fig. 7

Occupation density (normal approximation with at times and 5.0, with initial particle mass on the interval [12, 18]

The expressionappearing in (16) is equal to zero if and approaches or as depending on whether is negative or positive, that is according as or . Similarly,is equal to zero if and, as , approaches if or , and otherwise. Soapproaches 1/2 if or , 1 if and 0 otherwise. Thus we would expect the limiting profile for to be

Remark 4.1

Note that C(y, t), given by (17)

has the correct mass for all t, and

is a weak solution of the first-order partial differential Eq. (3).

It remains to establish the form (17) rigorously as . This is given in the following theorem.

Theorem 4.1

With the spacing set at , let be the expected number of particles at site at time t.

For all positive y and t, is given by (17).

For all positive y and t, where C(y, t) is defined in (17).

For and all positive t, let be the separation at time t of the agents that start at and . Then, for all , .

Proof

See the Appendix.

We started this section by observing that we can think of our initial condition as spreading a mass of over points, so that each point carries a mass of . In the limit as the mass at each point goes to zero but the number of points approaches infinity in such a way that the total mass approaches .

For , this motivates us to define in the model where the spacing is set to and is defined in part (iv) of Theorem 3.1. The limiting distribution of is given in the following theorem.

Theorem 4.2

For positive t, and ,

For positive t, and , where is the standard normal distribution function.

See the Appendix.

Remark 4.2

If , then Theorem 4.2(i) implies that for all . In particular, this implies that , which tells us that weakly converges to a point mass at zero. On the other hand, if then for all and weakly converges to a point mass at . Otherwise, weakly converges to a point mass at .

This agrees with the observation of Theorem 4.1 that the limiting expected number of agents C(y, t) has the form of a square wave of height over the interval .

Conclusion

In this paper we have analysed the cellular proliferation process model of Hywood et al. (2013a) by taking advantage of its relationship to the superposition of a set of coupled Yule-Furry processes, each with the same splitting rate. In Sect. 2, we presented the model of Hywood et al. (2013a) together with the argument of Hywood et al. (2013a) to derive an approximating second order PDE for the expected occupancy of the process at position y at time t.

In Sect. 3, we analysed the discrete model with positive , and used the model’s representation as a set of coupled Yule-Furry processes to derive an exact expression for the expected occupancy ) at site k at time t. We also derived expressions for the variance of the occupancy and the distribution of the number of sites between adjacent domain agents. Finally, we showed that at time t, domain agents occupy positions that are distributed as a geometric renewal process. In Sect. 4, we considered limiting versions as of the process that we analysed in Sect. 3. Theorem 4.1 contains expressions for the limiting expected occpuancy at position y at time t and its variance, while Theorem 4.2 presents Weak Law of Large Numbers and Central Limit Theorem results for the total mass of domain agents occupying positions to the left of y at time t.