A tutorial review of mathematical techniques for quantifying tumor heterogeneity

1.   Introduction

The rapid fall in the survivability of different important species is a major issue in ecology. The unplanned activity and indiscriminate use of biological resources lead to the unprecedented depletion in the abundance of species ^[1]. The detail study regarding conservation of different species is important since the loss of a strongly interactive species can lead to potential changes in ecosystem composition, structure, and diversity ^[2,3,4,5]. For example, the killing of large amount of great whales by industrial whaling makes changes in the krill-consumer dynamics in the Southern Ocean, and it is observed that, whaling is the main cause of megafaunal collapse in the North Pacific Ocean ^[6]. Moreover, the functional dynamics of coastal marine ecosystems worldwide have been extremely altered due to over fishing of large herbivores and predators ^[7]. It is frequently observed that the functional extinction of species interactions takes place well before the species themselves become extinct. In ocean ecosystem, large number of interactive species persists as rare adults, or small or juvenile individuals and do not interact as large adults. In addition, many strongly interactive species such as mammals are becoming extinct from the areas which they occupy one or two century ago ^[8]. The problem continues to persist since most conservation laws, such as the US Endangered Species Act (ESA) of 1973, fail to capture the effects of extinction of strongly interacting species. It is to be noted that the ESA does not include the inter-species interactions altogether in the criteria for recovery of endangered mammal species ^[4]. It only emphasizes the short-term, single-species, demographic viability in only a few areas. In reality, most of the current recovery plans do not call for increase in numbers of individuals or numbers of populations and geographic range ^[9,10]. Several works are also performed ^{[11,12,13,14]} to discuss this important issue. The persistence time of the population for a dynamical system is obtained in ^[15,16] which can be used to find the time to extinction of a population. However, as far as my knowledge goes the mathematical form to find the probability of extinction of a species belonging to any trophic level is still missing.

Stochastic population models serve as fundamental mathematical framework in modern ecological theory and is applied in various research areas such as population ecology ^[17], epidemiology ^[18], conservation biology ^[19] etc. Among all these areas, epidemiology is a major platform where the knowledge of stochasticity is extensively used in ^[20,21,22] and ^[23]. It is observed that discrete and continuous time Markov chain models and stochastic differential equation models are used in many areas of population biology. The Kolmogorov differential equations (often termed as Master equation ^[24]) are used to compute the transition probability distributions when the rates of different possible transitions are provided (e.g. birth and death rates).

It is frequently observed that to model a dynamical problem under stochastic environment, a system of Itostochastic differential equations is formed and studied. In many cases discrete stochastic models are developed by studying changes in the system components over a small interval of time such as birth-death process.

The prey-predator models regulated by deterministic differential or difference equations are crucial in quantitative studies of the dynamics of natural populations. In general, the dynamical characteristics of the systems, such as stability of equilibrium points, stability of limit cycles etc. can be predicted uniquely in the deterministic set up. However, to apply these models in natural populations, it is more realistic to consider the stochastic dynamics of the system, where the randomness is generated either from variation among individual growth rates or due to environmental fluctuations. For such cases, the characteristics of the deterministic models are replaced by the probability statements. As for example, we may consider stochastic version of the logistic model, where the equilibrium point is replaced by the equilibrium distribution, which is represented by the cloud points around the equilibrium point rather than a fixed value ^[25]. The difference is that, in logistic growth, the carrying capacity is the only stable point, hence the population never reaches zero for positive value of intrinsic growth rate, whereas, in stochastic environment there is always a probability associated with the extinction of the population. Thus the stochastic model framework is more appealing in real life scenarios ^[26,27].

The stability related issues of a predator-prey system is well discussed from the theoretical work of Lotka ^[28], Volterra ^[29] and the experimental work of Gause^[30]. In many cases the dynamics involving only two interacting species may not capture the real world scenario. The behavior of a complex system may only be understood when the interactions among a large number of species are incorporated ^[31,32]. Thus, it should be noted that what amount of risk of extinction the interacting species face for a given initial population size of prey and predator, is an important aspect of population dynamics. Here the risk of extinction is quantified as the probability of reaching the extinction equilibrium, before reaching to a sufficiently large population size so that extinction can be avoided. Theoretically, this large population tends to $\infty$^[33,34]. Recently ^[34] used the stochastic differential equation model including demographic stochasticity and predicted the time to extinction and extinction probability for the Atlantic Herring populations. They calculated the mean time to extinction from the distribution of the sojourn time, the amount of time spent by the species at each population size given an initial abundance ^[35].

In this paper, we explore the idea of obtaining the chance of extinction of species at any trophic level of the $n$ trophic food chain model. Discrete stochastic model is used to study the changes in the system components over a small time interval to develop the probability of extinction. The theoretical findings are validated by numerical simulation.

In next section (section $2$), the background of the knowledge to be applied is discussed. In section $3$, the main result of the paper is established with detail analysis followed by some case studies to understand the applicability of the findings. Section $4$ is the numerical part followed by section $5$ which deals with conclusion and future direction.

2.   Background of the work

In this section a two dimensional general predator prey system is considered as below,

where $x$ and $y$ denote the sizes of prey and predator population at time $t$, respectively; $b_1(.)$ is the birth rate of the prey species; $d_i(.) (i = 1, 2)$ are the natural death rates of the prey and predator populations respectively. The initial conditions are, $x(0)\ge 0, y(0)\ge 0$. The birth and death rates may be density dependent or density independent, described by the biological background of the two interacting species. The function $g(x, y)$ denotes the inter-species interaction between the prey and predators, known as predator's functional response. $\alpha$ is the conversion rate of ingested prey into new predators and assumed to lie between 0 and 1.

The stochastic formulation of the single species dynamics in absence of predators are studied by many authors ^[33,34]. In absence of predators, model (2.1) reduces to a single species model with birth and death rates $b_1(x)$ and $d_1(x)$ respectively. The logistic model under stochastic environment are studied in different ways in the literature ^{[36,37,38,39,40]}. In general, there are infinitely many choices of functions available for birth and death rates for different choices of parameters. The choices of the functions are mainly driven by the collected data on the species under investigation ^[41]. In addition, when we study real populations, generally some other information is also available such as, life expectancy. For example, the estimate of average life expectancy was used to obtain the estimates of intrinsic rate parameters^[37] to model the dynamics of muskrat population in Netherlands and it help them to uniquely chose the value of the rate parameters.

The generalized form of the model (2.1) in $n$-trophic food chain serves as the deterministic skeleton to formulate analogous stochastic models which account for the variation in births, deaths, transmission and recovery. In the stochastic formulation of the model (2.1) we assume $X(t)$ and $Y(t)$ are two discrete random variables representing the number of prey and predator populations respectively at time $t$ and these two random variables take values in the state space $\{(x, y) : x = 0, 1, \ldots M; y = 0, 1, \ldots N\}$ where $M$ and $N$ are the maximum sustainable population sizes of prey and predator population respectively. Here we assume that the time step $\Delta t$ is sufficiently small so that at most one change in the population size is possible in this small time gap, i. e. $\Delta X(t) = X(t+\Delta t)-X(t) \in \{-1, 0, 1\}$ and $\Delta Y(t) = Y(t+\Delta t)-Y(t) \in \{-1, 0, 1\}$. Thus the change in population sizes $\Delta X(t)$ and $\Delta Y(t)$ preclude the possibility of more than one birth, death, or any other transformations in small time $\Delta t$. In fact, these events have probabilities of order $(\Delta t)^2$. The infinitesimal probabilities associated with the birth and death of the populations for the model (2.1) is depicted in table 1. For example, $(\Delta X(t), \Delta Y(t)) = (1, 0)$ denotes the event of the birth of a prey and no change in predator abundance during time $\Delta t$. The probability of this event is obtained by multiplying the probability of one birth of prey ($b_1(x)\Delta t$) with the probability of no death of predator ($(1 - d_2(y)\Delta t)$). Similarly, the probabilities of the other events can be defined.

We extend this idea for a tri-trophic food chain model given by,

where the transitions and the corresponding probabilities are given in the table 2.

We extend this further for an $n$ trophic food chain model also as given below.

The $n$ trophic food chain model considered here is given by:

Here the functions $b_i(x_i)$'s, $d_i(x_i)$'s and $g_i(x_i, x_{i+1})$ denote the birth rates, death rates and the interaction terms respectively.

It is assumed that $X_i(t) (i = 1, 2, \ldots, n)$ are the discrete random variables representing the number of prey and predator populations in the food chain. $X(t)$ is the vector random variable which represents the population size at time $t$, that is $X(t) = (X_1(t), X_2(t), \ldots, X_n(t))$ and it takes values in the state space $\{x = (x_1, x_2, \ldots, x_n): x_i = 0, 1, \ldots M_i\}$. $M_i$ is the maximum sustainable population size of the population belonging to $i^{th}$ trophic level. As before the time step $\Delta t$ is taken sufficiently small so that only one change is possible in this small time interval, i. e. $\Delta X_i(t) = X_i(t+\Delta t)-X_i(t) \in \{-1, 0, 1\}$. The infinitesimal probabilities of the birth death formulation of model (2.3) is depicted in table $3$. For example, $\Delta X(t) = (1, 0, \ldots, 0)$ denotes the event of the birth of a prey and no change in the abundance of remaining species in time $\Delta t$. The probability of this event is ($b_1(x)\Delta t$). Similarly probabilities of the other events are defined in table 3.

It should be noted that the total number of transitions for an $n$ trophic food chain model is $2n+1$. The expected values of the changes in small time interval $\Delta t$ may be obtained as follows.

Similarly, we can find

Finally we have,

Now we will find the expressions for variance and co-variance terms for the changes.

We have,

Now it can be shown that,

It should be noted that, here the randomness in population growth rate is assumed to be demographic stochasticity which is the chance variation in the number of individual births and deaths and usually modeled by birth and death process. So the above model does not depict the random variations in the environmental conditions although the environmental fluctuation has significant effect on the population dynamics, for example it may have effect on the survival and reproduction rates ^[42].

3.   Chance of extinction

In this section the expression for the probability of extinction of a particular species population in the given food chain model is derived. Let $T_{x}(s)$ be the first passage time of the random variable $X$ to the state $s$. This is also known as persistent time or the first exit time in engineering literature. Thus the event ${T_{x}(0) < T_{x}(b)}$ represents the event that starting with population size $x$, the species goes to extinction before reaching to the population size $b$.

An $n$ dimensional food chain model as given below is considered here.

Suppose we are interested to find the chance of extinction of $i^{th}$ species. We start with the $1^{st}$ species (here the prey population) and subsequently do the same for $i^{th}$ species.

3.1.   Probability of extinction for prey population

Define,

In this case it should be noted that, according to the definition $u_1$ is not a random variable. We now derive the pde for $u_1$. It can be shown that,

see (2.4), (2.5), (2.6), (2.7), (2.11)

Now using the expressions of $\mathbb{E}(\Delta X_i)$, $\mathbb{E}(\Delta X_{i}^2)$ and $\mathbb{E}(\Delta X_i \Delta X_{i-1})$ and noting that

we get,

Theoretically, probability of extinction of the prey species can be obtained by letting $M_1 \longrightarrow \infty$. When the system is in steady state $\frac{\partial u_1}{\partial t} = 0.$

So, in steady state the above pde takes the form,

The boundary conditions are given by,

If the system starts with $0$ population for $1^{st}$ predator, the whole system becomes a single species dynamical system consisting of the prey population only. Because in this case all the predator population will die with probability 1. Hence, the said system will be governed by single species model and there will be no interaction term. In that case $u_1(x, \ldots, 0) = f_1(x)$. The expression for $f_1$ is given by ^[35],

It can be shown that $\phi_i$ satisfies (3.4) by setting $\Delta x^2_i = 0$ and $\Delta x_k \Delta x_i = 0$ $(k = 1, 2, \ldots, n, k\neq i)$ because we assume that when the $i^{th}$ species reaches the state $M_i$, the small change in the population size at $i^{th}$ trophic in time $\Delta t$ will be negligible compared to $M_i$.

3.2.   Probability of extinction of population at $i^{th}$ trophic level

To find the chance of extinction of the $i^{th}$ species we define,

where $X = (X_1, X_2, \ldots, X_n)$ and $x = (x_1, x_2, \ldots, x_n)$.

It should be noted that, in order to find the chance of extinction for $i^{th} (i = 2, 3, \ldots, n)$population, the condition that the $(i-1)^{th}$ species do not die out $(T_{x_{i-1}}(0) > T_{x_{i-1}}(M_{i-1}))$ is imposed. This condition is necessary because the loss of a species leads to the extinction of other dependant species ^[1]. $u_i(x)$ represents probability of extinction of $i^{th}$ species when $M_i \longrightarrow \infty.$ Now

where,

Following the same idea as applied in case of deriving the pde for $u_1$, it can be shown that $w_{i-1}$ satisfies the same pde (3.7) as $u_1$ does with the following boundary conditions.

This $\psi^k_{i-1}$ satisfies the pde obtained by setting $\Delta x^2_k = 0$ and $\Delta x_k \Delta x_j = 0$ in (3.4) for $j = 2, \ldots, n$ and $j\neq k$.

3.3.   Case study for two and three dimensional food chain model

So far the discussion revolves around the general $n$-trophic food chain model. In this section two particular cases viz. two dimensional and three dimensional food chain model for better understanding are discussed.

3.3.1.   Two dimensional food chain model

First, a two dimensional food chain model (2.1) is considered.

The probability of extinction for theses types of model is also provided in ^[43].

We define:

where $T_{x}(.)$ is defined earlier. $u(x, y)$ represents the probability of extinction when the dynamical system start from the state $(x, y)$ allowing $M\longrightarrow \infty$ ^[33]. Now,

From ^[35], it can be shown that

Now,

Taking expectation on both sides and using (3.22) we get,

Using the equation (2.1), we obtain the following equation

Now replacing the coefficients of $\frac{\partial u}{\partial x}$ and $\frac{\partial u}{\partial y}$ by $\frac{\mathrm{d}x}{\mathrm{d}t}$ and $\frac{\mathrm{d}y}{\mathrm{d}t}$ respectively, by virtue of the equation (2.1), and noting that, $\frac{\partial u}{\partial t} = \frac{\partial u}{\partial x}\frac{\mathrm{d}x}{\mathrm{d}t}+\frac{\partial u}{\partial y}\frac{\mathrm{d}y}{\mathrm{d}t}$, we obtain the following partial differential equation

Therefore the above expression can be written as

In steady state the pde (3.24) takes the form:

Lemma 3.1. The pde (3.25) is elliptic in nature.

Proof. The standard form of a second order partial differential equation is

In this case $A = \frac{1}{2}[b_1(x)+ g(x, y)+d_1(x)]$, $B = -\alpha g(x, y)$ and $C = \frac{1}{2}[\alpha g(x, y)+d_2(y)]$. Therefore, $B^2 - 4AC = -\alpha(1-\alpha) g^2(x, y)-d_2(y)(b_1(x)+d_1(x))-g(x, y)(\alpha b_1(x)+\alpha d_1(x)+d_2(y)) < 0$ since $g(x, y), b_1(x), d_1(x), d_2(y)$ are all positive and $\alpha < 1$.

Hence the pde (3.25) is elliptic.

To compute the first passage probabilities given the size of prey and predator, the above equation is solved numerically. We now proceed to develop appropriate boundary conditions for the above differential equation. We have,

Therefore

because, if the initial prey population is $0$, the event $T_x(0) < T_x(M)$ is a certain event.

Also,

The function $f(x)$ is given by (^[35])

Now

The differential equation for $g(x) = u(x, N)$ can be derived by setting the coefficients of $\frac{\partial^2 u}{\partial y^2}$ and $\frac{\partial^2 u}{\partial x \partial y}$ in (3.23) to be $0$ as $N$ is assumed to be the maximum possible size of predator population. In addition to this $\frac{\partial u}{\partial t} = 0$ for stationary state. Thus,

The boundary conditions are $g(0) = u(0, N) = 1, g(M) = u(M, N) = 0.$ Using this, the solution of the above differential equation is,

Remark 1. We have,

This indicates that if the predator population is at its maximum level, the extinction of the prey population is inevitable if the prey population is sufficiently small compared to the maximum size ($M$) of the prey population.

Now we will develop the method of finding the probability of extinction of predator population. We define

The condition $T_x(0) > T_x(M)$ is imposed since otherwise the extinction of predator species is certain. We can write,

where

It can be shown that (in a similar way as in the case of $u$) $w$ satisfies the same pde as $u$ with the boundary conditions given next.

We have,

$w(M, y)$ will be obtained as in case of $u(x, N)$ and is found to be,

Remark 2. We have

Now it should be noted that,

The implication of this result is that if the prey population size is at its maximum, the predator population will reach state $0$ before it reaches its maximum (i. e. $N$) if $N$ is very large compared to the initial population size of predator.

The remaining boundary conditions are,

3.3.2.   Three dimensional food chain model

In many cases tri-trophic food chain model is appropriate to capture more information for a complex system ^[44]. In order to find the chance of extinction of a population at any trophic level, we can extend the above idea for tri-trophic food chain model also. Here a tri-trophic food chain model of the form given below is considered,

Similar procedure may be used to define,

We can follow the same argument as in the case of general $n$-trophic food chain model to show that,

The six boundary conditions on the six surfaces are given by,

The function $f_1(x, z)$ can be found in the following way:

Since $M_2$ is assumed to be the maximum possible population size of $1^{st}$ predator, the decrease of it by one individual is negligible compared to $M_2$. Hence we may set $\Delta x_2 = 0$ and $(\Delta x_2)^2 = 0$ in the general equation (3.4) for $1^{st}$ predator.

This will give the following differential equation in stationary state,

$f_1$ satisfies the boundary conditions,

In case of $u_1(x, y, 0)$, the system is free from 3rd predator and hence it is equivalent to 2dimensional prey-predator system. This suggest that $u_1(x, y, 0) = u(x, y) = f_2(x, y)$.

Following the argument in order to derive the differential equation for $f_3$, it can be shown that,

The boundary conditions for (3.51) are,

In case of 1st predator population we define,

As in the case of two dimensional model, the expression for $u_2$ can be written in the form,

where,

$w_1$ satisfies the same pde as $u_1$ with suitable boundary conditions given below.

$w_1(0, y, z) = 0$. From the definition of $w_1$ we have,

Hence we have $w_1(x, 0, z) = 1-u_1(x, 0, z) = 1-f(x)$. It can be shown that $w_1(x, M_2, z) = 0.$ In order to find $w_1(x, y, 0)$ we follow the following procedure.

Let $w_1(x, y, 0) = \phi(x, y)$(say). It can be shown that $\phi(x, y)$ satisfies the pde (3.3.2) with the coefficients of $\frac{\partial^2\phi}{\partial z^2}$, $\frac{\partial^2\phi}{\partial y\partial z}$ to be $0$. Hence $\phi(x, y)$ satisfies at steady state,

$\phi(x, y)$ satisfies the following boundary conditions. $\phi(x, 0) = w_1(x, 0, 0) = 1-f(x)$ using the expression for $w_1(x, 0, z)$ given above. $\phi(x, M_2) = 0$, $\phi(0, y) = 0$ and $\phi(M_1, y) = 1-\frac{y}{M_2}$.

Similarly, $w_1(x, y, M_3) = \psi(x, y)$ can be shown to satisfy

The boundary conditions for $\psi$ are:

which can be obtained in similar way.

In continuation of this process of formulating the probability of extinction we can find the same ($u_3$) for the 2nd predator also.

4.   Numerical simulations

In this section the theoretical findings are validated by using the following simple two dimensional Lotka-Volterra food chain model.

Here $r$ is the intrinsic growth rate of the prey population, $K$ is the carrying capacity, $\alpha$ is the conversion coefficient and $d$ is the natural death rate of the predator population. If we compare the model (4.1) with the model (2.1), we may consider the following birth-death and interaction terms.

The probability of extinction of prey population ($u$) and predator population ($v$) are generated numerically taking help of finite difference method using the procedure as given in section (3.3.1). Sufficiently large values of $M$ and $N$ are considered to have better approximation for probability of extinction. The surface of $u(x, y)$ and $v(x, y)$ are given in the figures 1 and 2.

One can notice in the figure 1 that the probability of extinction is high in the region where the initial population sizes of both prey and predator population are low and vice-versa. Note that it is natural for any species to be prone to extinction if its initial population size is very small. This natural phenomena is captured in this proposed method. The solution set of $u(x, y)$ and $v(x, y)$ provide the probability of extinction of prey and predator population respectively.

To extend the applicability of this method or to understand the theoretical findings in a better way, we also consider the effect of initial population size of either prey or predator separately on the extinction probability. We observe that if the prey population size increases, the probability of extinction of prey population diminishes for any fixed initial population size of predator(Figure 3). The interesting result we observe that, if the initial size of predator increases, the chance of extinction of prey decreases until the initial size of predator population reaches a certain value. The further increase in the initial size of predator increases the chance of extinction of prey population (see Figure 3 and 4). Thus the presence of a moderate amount of predator in the system is important to avoid extinction of prey which in turn helps the predator to survive also.

In figure 5, we observe that if the initial prey population increases, the chance of extinction ($v$ here) of predator decreases rapidly up to a certain value of initial prey. The change in the extinction probability of predator is not significant for further increase in the prey population and when initial prey population size exceeds a certain value, the chance of extinction of predator increases. Here also we observe that if a food chain model start with very high abundance of initial prey population predator population may become extinct with high probability. We also observe that for any fixed value of initial prey population, the higher the initial predator population, lower is the chance of extinction of it.

5.   Conclusion and future direction

Indiscriminate use of biological resources along with unplanned activity of human ^[45] has adverse effect on ecosystem. As a result, a good number of rare and endangered species are becoming extinct very rapidly. If this continues to happen it can be predicted that the species in the ecosystem will be extinct in an unprecedented rate ^[46,47]. In this context, the assessment of extinction threats of different rare and endangered species are needed so that corresponding safeguard can be provided.

In this paper we have formulated the technique to find the probability of extinction for an $n$-trophic food chain model. This process can be used to find the extinction probability of species belonging to any trophic level of the food chain and proper conservation policy can be implemented in advance. This strategy may help significantly to account for the extinction threat of species in widely recognized bio-diverse places like Atlantic Forest ^[48], Sumatra and Madagascar etc. ^[49]. It is observed that, the variables determining the extinction probabilities follow few partial differential equations some of which are very complex in nature and can not be solved explicitly. The complexity increases with the increase in the dimension of food chain. So numerical method may be used to find the solution of those differential equations to get the probability of extinction.

In this paper we consider only demographic variability. The impact of environmental noise on any food chain is also significant. It will be interesting and more appropriate if this work is extended incorporating environmental noise also.

Acknowledgments

I am thankful to Govt. College of Engg. & Textile Technology, Berhampore for this research. I also express my thank to fellow researchers with whom I discussed and take suggestions for this work and to the learned reviewers whose suggestions help me immensely to improve the manuscript.

Conflict of interest

The authors declared they have no competing interests.