# Linear Processes in Stochastic Population Dynamics: Theory and Application to Insect Development

Authors: Hernán G. Solari; Mario A. Natiello
Journal: The Scientific World Journal
 (2014)

---

## Abstract

We consider stochastic population processes (Markov jump processes) that develop as a consequence of the occurrence of random events at random time intervals. The population is divided into subpopulations or compartments.
The events occur at rates that depend linearly on the number of individuals in the different described compartments. The dynamics is presented in terms of  Kolmogorov Forward Equation in the space of events and projected onto the space of populations when needed. The general properties of the problem are discussed. Solutions are obtained using a revised version of the Method of Characteristics. After a few examples of exact solutions we systematically develop short-time approximations to the problem. While the lowest order approximation matches the Poisson and multinomial heuristics previously proposed, higher-order approximations are completely new. Further, we analyse a model for insect development as a sequence ofE developmental stages regulated by rates that are linear in the implied subpopulations. Transition to the next stage competes with death at all times. The process ends at a predetermined stage, for example, pupation or adult emergence. In its simpler version all the stages are distributed with the same characteristic time.


---

## Body

## 1. Introduction

Stochastic population dynamics deals with the (stochastic) time evolution of groups of similar individuals (called populations or subpopulations) immersed in a habitat. The goal is to describe the evolution of population numbers and the number of events associated with changes in the populations.Most of the attention has been paied in the past to the populations, considering events only by their extrinsic value as jumps in populations. However, this perspective must be revised. In any study of vector-transmitted diseases, such as in [1, 2], the most relevant statistics corresponds to the interactions between hosts and vectors. Such interactions are associated with events in the life cycle of the vector. In the case of arthropods transmitting diseases to mammals, blood feeding to complete oogenesis is the opportunity to transmit pathogens. In other problems, the development of insects is experimentally described by the sequences of transformations in their life cycle (such as moulting) [3–5]: each transformation must be considered an event. To our knowledge, explicit use of the event structure to study population changes in Markov Jump processes has been introduced in [6], along with a short-time approximation based on self-consistent Poisson processes.The present research originates in our efforts to model vector-transmitted diseases, such as Dengue, starting from a description of the life cycle ofAedes aegypti, the vector [7–9]. In modeling insect development in an aquatic phase (such as the case of Aedes aegypti and Drosophila melanogaster) a crucial point is to represent accurately the statistics of the duration of the aquatic phase (also known as the time of emergence statistics) as it determines the response of the population to climatic events such as rains [10]. A deterministic model for development was proposed in [11] as a sequential process with associated linear rates. The final sections of this paper describe a stochastic development model.While [6] deals with the stochastic approximation of general problems (achieving a method with error o
(
t
3
/
2
)), this paper specifically addresses the problem of formulating and solving population dynamics problems in event space with linear rates, that is, rates that depend linearly on the populations (in the chemical literature they are known as first-order reactions [12]). We provide both exact solutions and general approximation techniques with error o
(
t
n
). Results from [6] for the linear case are a modified version of one of the approximation schemes in this paper (Poisson approximation, Section 5). The linear problem is mathematically simpler than the general case and a few general results can be established before resorting to approximations.Linear processes are an important part of most biological dynamical descriptions. Using an analogy from classical population dynamics, linear (Malthusian) population growth [13] describes the evolution of a population at low densities (e.g., where the local environment is essentially unaffected by the population growth and the individuals do not interfere with each other). More generally, processes that consider only individuals in a somewhat isolated basis can be satisfactorily described with linear rates. However, whether to use linear descriptions or not is ultimately determined by the nature of the process (or subprocess). For example, finite environments will sooner or later force individuals to interact in growing populations and nonlinear rates become necessary such as in Verhulst’s description of population growth [14].Perhaps the oldest application of linear processes to population dynamics regardsbirth and death processes. Historically, according to Kendall [15], the first to propose birth and death equations in population dynamics was Furry [16], for a system of physical origin. Currently, such an approach is being used in several sciences: linear stochastic processes have been proposed for the development of tumors [17, 18]; drug delivery and cell replication are also described in these terms [19]. Concerning the tumors, the model proposed in [17, 18] is based on linear individual processes, a detailed treatment of which would require a paper of its own. The accumulation of individual processes is finally approximated with a “tunnelling process” corresponding to a maturation cascade that can be easily described with the tools of this paper (see (71)). Concerning the question of drug delivery, we will come back in detail to [19] in Section 7.A substantial use of population dynamics driven by jump processes has been done also in chemistry, particularly after the rediscovery by Gillespie [20] of Kendall’s simulation method [21, 22]. Poisson [23, 24], binomial [25, 26], and multinomial [27] heuristics for the short-time integrals have been proposed. There is, however, an important difference between chemistry and biology concerning the issues of this paper. On the one hand, linear rates are seldom the case in chemical reactions except perhaps in the case of irreversible decompositions or isomeric recombination. Attention to linear rates in chemistry is relatively recent [12] when compared to over 200 years of Malthusian population growth. On the other hand, detailed event statistics does not appear to be fundamental either: highly diluted chemical reactions on a test tube may easily generate 1
0
14 events of each type. Consequently, attention to event based descriptions appears to be absent from the chemistry literature. This might, however, change in the future, since the technical possibilities given by for example, “in-cell chemistry” and “nanochemistry” allow dealing with a lower number of events where tiny differences in event count could be relevant for the general outcome.In the mathematical literature attention to jump processes starts with Kolmogorov’s foundational work [28] and its further elaboration by Feller [29]. A substantial effort to relate the stochastic description to deterministic equations was performed by Kurtz [30–34]. However, attention has always focused on the dynamics in population space, rather than event space.This work focuses on the event statistics associated withlinear processes. Even for cases where the detailed event statistics is not of central interest, such an approach offers an alternative to the more traditional formulation in population space.Further, we discuss in detail the stochastic version of a linear developmental model. Specifically, we deal with a subprocess in the life cycle of insects that admits both a linear description and an exact solution, namely, the development of immature stages, which is a highly individual subprocess. From egg form to adult form, the development of insects goes through several transformations at different levels. For example, at the most visible level, insects undergoing complete metamorphosis evolve in the sequence egg, larva, pupa, and adult. Further, different substages can be recognised by inspection in the development of larvae (instars) and even more stages when observed by other methods [35]. The immature individual may die along the process or otherwise complete it and exit as adult.In Section2 we formulate the dynamical population problem in event space. In Section 2.1 we present the Kolmogorov Forward Equations (for Markov jump processes) for the most general linear systems of population dynamics in the form appropriated to probabilities and generating functions. The general solution of these equations will be written using the Theorem of the Characteristics, reformulated with respect to the standard geometrical presentation [36] into a dynamical presentation. We will analyse the general structure of the equations and in particular their first integrals of motion and the projection onto population space (Sections 2 and 3). Section 4 contains two classical examples. In Section 5, Poisson and multinomial approximations (both with error o
(
t
)) are derived and a higher-order approximation with error o
(
t
n
) (with n as an arbitrary positive integer) is elaborated. The larval development model is developed in Sections 6 and 7, treating the two general cases. We also revisit in Section 7 the question of drug delivery introduced in [19], in terms of the present approach. Concluding remarks are left for the final Section.


## 2. Mathematical Formulation

The main tools of stochastic population dynamics are a list of (nonnegative, integer) populationsX
j, where j
=
1
⋯
N a list of events α
=
1
⋯
E an (integer) incidence matrix
δ
j
α describing how event α modifies population j, and a list of transition rates W
α
(
X
1
,
…
,
X
n
,
t
) describing the probability rate per unit time of occurrence of each event.An important feature of the present approach is that the dynamical description will be handled in event-space. In this sense, the “state” of the system is given by the arrayn
(
s
,
t
)
=
(
n
1
,
…
,
n
E
)
t indicating how may events of each class have occurred from time s up to time t. Hence, event-space is, in this setup, a nonnegative integer lattice of dimension E. The connection with the actual population values at time t is (s denotes the initial condition)
(1)
X
j
(
t
)
=
X
j
(
s
)
+
∑
α
‍
δ
j
α
n
α
(
s
,
t
)
.Definition 1 (linear independency).
One says that the populations arelinearly dependent if there is a vector u
≠
0 such that for all t
≥
s one has ∑
j
‍
u
j
X
j
(
t
)
=
0. Otherwise, one says that the populations are linearly independent and abbreviate it as LI.
In similar form, we say that thepopulations have linearly dependent increments if there is u
≠
0 so that for all t
′
>
t
≥
s it holds that ∑
j
‍
u
j
(
X
j
(
t
′
)
-
X
j
(
t
)
)
=
0.
The populations are said to havelinearly independent increments (abbreviated LII) if they do not have linearly dependent increments.Along this work we will assume the following.Assumption 2.
The populations havelinearly independent increments; that is, all involved populations are necessary to provide a complete description of the problem; none of them can be expressed as a linear combination of the remaining ones for all times.Lemma 3.
(a) If the populations have linearly independent increments, then the populations are linearly independent.
(b) The matrixδ has no nonzero vectors u such that ∑
j
‍
δ
j
α
u
j
=
0 for all α if and only if the populations have linearly independent increments.Proof.
(a) The assertion is equivalent to the following: if the populations are linearly dependent, then they have linearly dependent increments, which is immediate after the definition of linearly dependent increments.
(b) Assume that for allα there exists u
≠
0 such that ∑
j
‍
u
j
δ
j
α
=
0; then
(2)
∑
j
‍
u
j
(
X
j
(
t
′
)
-
X
j
(
t
)
)
=
∑
j
‍
u
j
(
∑
α
‍
δ
j
α
n
α
(
t
′
,
t
)
)
=
0
,
and the populations have linearly dependent increments.
Conversely, assume that the populations have linearly dependent increments and select time intervals such thatn
β
(
t
β
′
,
t
β
)
=
1 and n
α
(
t
β
′
,
t
β
)
=
0 for α
≠
β; then, we have that
(3)
∑
j
‍
u
j
δ
j
α
=
0
,
∀
α
.
Statement (b) is the negative statement of these facts.Strictly speaking the last step in the above proof takes for granted that there exist time intervals where only one event occurs. This is assured by the next assumption (third item). Also the initial conditions and the problem description should be such that all defined events are relevant for the dynamics (e.g., the event “contagion” is irrelevant for an epidemic problem with zero infectives as initial condition).


### 2.1. Dynamical Equation—Basic Assumptions

Population numbers evolve in “jumps” given by the occurrence of each event (birth, death, etc.). In the sequel, we will use Latin indicesi, j, k, and m for populations and Greek indices α, β, and γ for events.The ultimate goal of stochastic population dynamics is to calculate the probabilitiesP
(
n
1
,
…
,
n
E
,
s
,
t
) of having n
1
,
…
,
n
E events in each class up to time t, given an initial condition {
X
j
(
s
)
}.The stochastic dynamics we are interested in correspond to aDensity Dependent Markov Jump Process [37, 38]. For a sufficiently small time interval h, we consider the following.Assumption 4.
For each event,(i)
event occurrences in disjoint time intervals are independent;
(ii)
the Chapman-Kolmogorov equation [28, 29] holds;
(iii)
P
(
n
α
+
1
,
t
+
h
,
t
)
=
W
α
(
X
,
t
)
h
+
o
(
h
);
(iv)
P
(
n
α
,
t
+
h
,
t
)
=
1
-
W
α
(
X
,
t
)
h
+
o
(
h
);
(v)
P
(
n
α
(
h
)
+
k
,
t
+
h
,
t
)
=
o
(
h
), k
≥
2. Here h is the elapsed time since the time t and n
α are the events of type α accumulated up to t.It goes without saying that probabilitiesP
(
⋯
) are assumed to be differentiable functions of time. W
α
(
X
,
t
) denotes the probability rate for event α, which is assumed in the following lemmas to be a smooth function of the populations, not necessarily linear (yet). Two well-known general results describe the dynamical evolution under these assumptions.Lemma 5 (see Theorem1 in [29]).
Under Assumptions2 and 4, the waiting time to the next event is exponentially distributed.Theorem 6 (see [28] and Theorem  1 in [29]).
Under Assumptions2 and 4, the dynamics of the process in the space of events obeys the Kolmogorov Forward Equation, namely;
(4)
P
˙
(
{
n
1
,
…
,
n
E
}
,
s
,
t
)
=
∑
α
‍
P
(
{
…
,
n
α
-
1
,
…
}
,
s
,
t
)
W
α
(
X
(
{
…
,
n
α
-
1
,
…
,
t
}
)
)
-
P
(
{
n
1
,
…
,
n
E
}
,
s
,
t
)
(
∑
α
‍
W
α
(
X
(
{
n
1
,
…
,
n
E
,
t
}
)
)
)
.Proofs of these results are given in AppendixA.


### 2.2.  Linear Rates

We will assume that the transition rates depend—in a linear fashion—on the populations only. To distinguish from general rates (denoted above by the letterW), we will use the letter R to denote linear transition rates. Hence,
(5)
R
α
=
∑
k
‍
r
α
k
X
k
(
t
)
.
Since rates are nonnegative, we have that the proportionality coefficient between event α and subpopulation k satisfies r
α
k
≥
0. These proportionality coefficients convey the environmental information and as such they may be time-dependent since the environment varies with time. Whenever possible, we do not write this time dependency explicitly to lighten the notation.Lemma 7.
For the case of linear rates, population space{
X
k
≥
0
} is positively invariant under the time evolution given by (1) and the previous assumptions if and only if δ
i
α
≥
-
1 and δ
i
α
=
-
1 only if R
α
=
r
α
i
X
i
(
t
) (no sum).Proof.
The probability of occurrence of eventα in a small time interval h is
(6)
h
R
α
+
o
(
h
)
=
h
∑
k
‍
r
α
k
X
k
(
t
)
+
o
(
h
)
.
After an occurrence, population X
k is modified to X
k
+
δ
k
α. Population space is positively invariant under time evolution if and only if the probability of occurrence of event α is zero for those population values X
¯
k such that X
¯
k
+
δ
k
α
<
0 (we are just saying that events pushing the populations to negative values do not exist). Since this should hold for arbitrary h, we have that R
α
=
0.
Assume thatδ
k
α
≤
-
1. Then, R
α
=
0 for 0
≤
X
¯
k
<
-
δ
k
α. Letting X
¯
k
=
0 we realise that ∑
j
≠
k
‍
r
α
j
X
j
(
t
)
=
0 regardless of the values of X
j, and hence it holds that r
α
j
=
0 for j
≠
k and therefore R
α
=
r
α
k
X
k (no sum). This proves the second statement. If δ
k
α
<
-
1 we have that R
α
=
0 also for X
¯
k
=
1. Hence r
α
k
=
0 as well; R
α
≡
0 and the event α can be ignored since it does not influence the dynamics.


## 2.1. Dynamical Equation—Basic Assumptions

Population numbers evolve in “jumps” given by the occurrence of each event (birth, death, etc.). In the sequel, we will use Latin indicesi, j, k, and m for populations and Greek indices α, β, and γ for events.The ultimate goal of stochastic population dynamics is to calculate the probabilitiesP
(
n
1
,
…
,
n
E
,
s
,
t
) of having n
1
,
…
,
n
E events in each class up to time t, given an initial condition {
X
j
(
s
)
}.The stochastic dynamics we are interested in correspond to aDensity Dependent Markov Jump Process [37, 38]. For a sufficiently small time interval h, we consider the following.Assumption 4.
For each event,(i)
event occurrences in disjoint time intervals are independent;
(ii)
the Chapman-Kolmogorov equation [28, 29] holds;
(iii)
P
(
n
α
+
1
,
t
+
h
,
t
)
=
W
α
(
X
,
t
)
h
+
o
(
h
);
(iv)
P
(
n
α
,
t
+
h
,
t
)
=
1
-
W
α
(
X
,
t
)
h
+
o
(
h
);
(v)
P
(
n
α
(
h
)
+
k
,
t
+
h
,
t
)
=
o
(
h
), k
≥
2. Here h is the elapsed time since the time t and n
α are the events of type α accumulated up to t.It goes without saying that probabilitiesP
(
⋯
) are assumed to be differentiable functions of time. W
α
(
X
,
t
) denotes the probability rate for event α, which is assumed in the following lemmas to be a smooth function of the populations, not necessarily linear (yet). Two well-known general results describe the dynamical evolution under these assumptions.Lemma 5 (see Theorem1 in [29]).
Under Assumptions2 and 4, the waiting time to the next event is exponentially distributed.Theorem 6 (see [28] and Theorem  1 in [29]).
Under Assumptions2 and 4, the dynamics of the process in the space of events obeys the Kolmogorov Forward Equation, namely;
(4)
P
˙
(
{
n
1
,
…
,
n
E
}
,
s
,
t
)
=
∑
α
‍
P
(
{
…
,
n
α
-
1
,
…
}
,
s
,
t
)
W
α
(
X
(
{
…
,
n
α
-
1
,
…
,
t
}
)
)
-
P
(
{
n
1
,
…
,
n
E
}
,
s
,
t
)
(
∑
α
‍
W
α
(
X
(
{
n
1
,
…
,
n
E
,
t
}
)
)
)
.Proofs of these results are given in AppendixA.


## 2.2.  Linear Rates

We will assume that the transition rates depend—in a linear fashion—on the populations only. To distinguish from general rates (denoted above by the letterW), we will use the letter R to denote linear transition rates. Hence,
(5)
R
α
=
∑
k
‍
r
α
k
X
k
(
t
)
.
Since rates are nonnegative, we have that the proportionality coefficient between event α and subpopulation k satisfies r
α
k
≥
0. These proportionality coefficients convey the environmental information and as such they may be time-dependent since the environment varies with time. Whenever possible, we do not write this time dependency explicitly to lighten the notation.Lemma 7.
For the case of linear rates, population space{
X
k
≥
0
} is positively invariant under the time evolution given by (1) and the previous assumptions if and only if δ
i
α
≥
-
1 and δ
i
α
=
-
1 only if R
α
=
r
α
i
X
i
(
t
) (no sum).Proof.
The probability of occurrence of eventα in a small time interval h is
(6)
h
R
α
+
o
(
h
)
=
h
∑
k
‍
r
α
k
X
k
(
t
)
+
o
(
h
)
.
After an occurrence, population X
k is modified to X
k
+
δ
k
α. Population space is positively invariant under time evolution if and only if the probability of occurrence of event α is zero for those population values X
¯
k such that X
¯
k
+
δ
k
α
<
0 (we are just saying that events pushing the populations to negative values do not exist). Since this should hold for arbitrary h, we have that R
α
=
0.
Assume thatδ
k
α
≤
-
1. Then, R
α
=
0 for 0
≤
X
¯
k
<
-
δ
k
α. Letting X
¯
k
=
0 we realise that ∑
j
≠
k
‍
r
α
j
X
j
(
t
)
=
0 regardless of the values of X
j, and hence it holds that r
α
j
=
0 for j
≠
k and therefore R
α
=
r
α
k
X
k (no sum). This proves the second statement. If δ
k
α
<
-
1 we have that R
α
=
0 also for X
¯
k
=
1. Hence r
α
k
=
0 as well; R
α
≡
0 and the event α can be ignored since it does not influence the dynamics.


## 3. Kolmogorov Forward Equation

Let us recast the dynamical equation given by Theorem6 in terms of the generating function [39] in event-space, defined as
(7)
Ψ
(
z
1
,
…
,
z
E
,
s
,
t
)
=
∑
{
n
α
}
‍
(
∏
α
‍
z
α
n
α
)
P
(
n
1
,
…
,
n
E
,
s
,
t
)
.
For the case of linear rates, we define
(8)
R
α
0
=
∑
k
‍
r
α
k
X
k
(
s
)
,
F
α
β
=
∑
k
‍
r
α
k
δ
k
β
.
As a consequence,
(9)
R
α
=
R
α
0
+
∑
β
‍
F
α
β
n
β
.
Using Theorem 6 to rewrite ∂
Ψ
/
∂
t and using further n
β
Ψ
=
z
β
(
∂
Ψ
/
∂
z
β
) (a matter of replacing in the definition of Ψ), we obtain
(10)
∂
Ψ
∂
t
|
s
=
∑
α
‍
(
z
α
-
1
)
(
R
α
0
Ψ
+
∑
β
‍
F
α
β
z
β
∂
Ψ
/
∂
z
β
)
.
Since the generating function reflects a probability distribution, we have the probability condition Ψ
(
1
,
…
1
,
s
,
t
)
=
1. Further, the solutions of (10) can be uniquely translated into those given by Theorem 6 for corresponding initial conditions: P
(
m
1
,
…
,
m
E
,
s
,
s
)
=
1 corresponding to Ψ
(
z
1
,
…
,
z
E
,
s
,
s
)
=
z
1
m
1
z
2
m
2
⋯
z
E
m
E, where the integers m
1
,
…
,
m
E denote how many events of each class had already occurred at t
=
s. The natural initial condition becomes Ψ
(
z
1
,
…
,
z
E
,
s
,
s
)
=
1, corresponding to P
(
0
,
…
,
0
,
s
,
s
)
=
1.


### 3.1. The Method of Characteristics Revisited

Lemma 8.
Equation (10) with R
α
0
=
0 and with the initial condition Ψ
(
z
1
,
…
,
z
E
,
s
,
s
)
=
Φ
0
(
z
1
,
…
,
z
E
) admits the solution
(11)
Ψ
(
z
1
,
…
,
z
E
,
s
,
t
)
=
Φ
0
(
ϕ
(
z
,
t
,
t
-
s
)
)
,
where ϕ
(
z
,
t
,
τ
) is the flow of the following ODE with initial condition x
(
0
)
=
z, s
(
0
)
=
t:
(12)
d
x
β
d
τ
=
∑
α
‍
F
α
β
(
x
α
-
1
)
x
β
≡
b
β
(
x
,
s
)
,
d
s
d
τ
=
-
1
.Proof.
It is convenient to renames as s
=
x
E
+
1, introduce b
E
+
1
=
-
1, and recall that b
β depends only on the extended variable array x. However, when necessary, we will write explicitly x
(
0
) as (
z
,
t
). We have that
(13)
∂
Ψ
∂
t
|
s
=
∑
β
E
+
1
‍
∂
Φ
0
∂
x
β
|
ϕ
(
x
,
t
-
s
)
d
ϕ
β
(
z
,
t
,
t
-
s
)
d
t
=
∑
β
E
+
1
‍
∂
Φ
0
∂
x
β
|
ϕ
(
x
,
t
-
s
)
×
(
b
β
(
ϕ
(
z
,
t
,
t
-
s
)
)
+
∂
ϕ
β
(
z
,
u
,
t
-
s
)
∂
u
|
u
=
t
)
.
Further, consider the monodromy matrix [40]
(14)
M
γ
β
=
∂
ϕ
β
(
x
,
τ
)
∂
x
γ
.
M transforms the initial velocity into the final velocity; hence
(15)
∑
γ
‍
M
γ
β
b
γ
(
z
,
t
)
=
b
β
(
ϕ
(
z
,
t
,
τ
)
)
.
Thus, we get
(16)
∂
Ψ
∂
t
|
s
=
∑
β
E
+
1
‍
∂
Φ
0
∂
x
β
|
ϕ
(
x
,
t
-
s
)
(
∑
γ
E
+
1
‍
M
γ
β
b
γ
(
z
,
t
)
+
∂
ϕ
β
(
z
,
u
,
t
-
s
)
∂
u
|
u
=
t
)
=
∑
β
E
+
1
‍
∂
Φ
0
∂
x
β
|
ϕ
(
x
,
t
-
s
)
×
(
∑
γ
E
‍
M
γ
β
b
γ
(
z
,
t
)
-
M
E
+
1
β
+
∂
ϕ
β
(
z
,
u
,
t
-
s
)
∂
u
|
u
=
t
)
=
∑
γ
E
‍
∂
Ψ
∂
z
γ
b
γ
(
z
,
t
)
.
From the definition of monodromy matrix above it follows that
(17)
-
M
E
+
1
β
+
∂
ϕ
β
(
z
,
u
,
t
-
s
)
∂
u
|
u
=
t
=
0
,
and by (11) and the chain rule,
(18)
∑
β
‍
∂
Φ
0
∂
x
β
|
x
=
ϕ
(
x
,
t
-
s
)
M
γ
β
=
∑
β
‍
∂
Φ
0
∂
x
β
|
x
=
ϕ
(
x
,
t
-
s
)
∂
ϕ
β
(
z
,
t
,
t
-
s
)
∂
z
γ
=
∂
Ψ
∂
z
γ
.Lemma 9.
Equation (10) with R
α
0
≠
0 and with the initial condition Ψ
(
z
1
,
…
,
z
E
,
s
,
s
)
=
Φ
0
(
z
1
,
…
,
z
E
) admits the solution
(19)
Ψ
(
z
1
,
…
,
z
E
,
s
,
t
)
=
Φ
0
(
ϕ
(
z
,
t
,
t
-
s
)
)
Φ
1
(
z
1
,
…
,
z
E
,
s
,
t
)
,
where Φ
0
(
ϕ
(
z
,
t
,
t
-
s
)
) is a solution of the homogeneous problem with R
α
0
=
0 and
(20)
Φ
1
(
z
1
,
…
,
z
E
,
s
,
t
)
=
exp
⁡
∫
s
t
‍
(
∑
α
‍
R
α
0
(
τ
)
(
ϕ
α
(
z
,
t
,
t
-
τ
)
-
1
)
d
τ
)
.Proof.
The proof is nothing more than an application of the method of the variation of the constants. The details are as follows.
Introducing the proposed solutionΨ into (10) and using Lemma 8 to eliminate Φ
0 we obtain that Φ
1 satisfies (10) with initial condition Φ
1
(
z
,
s
,
s
)
=
1. Computing directly from (20) we have that
(21)
∂
Φ
1
∂
t
|
s
=
∑
α
‍
R
α
0
(
ϕ
α
(
z
,
t
,
0
)
-
1
)
Φ
1
+
∑
β
‍
∂
Φ
1
∂
ϕ
β
∂
ϕ
β
(
z
,
t
,
t
-
τ
)
∂
t
=
∑
α
‍
R
α
0
(
z
α
-
1
)
Φ
1
+
∑
β
‍
∂
Φ
1
∂
ϕ
β
(
b
β
(
ϕ
(
z
,
t
,
t
-
τ
)
)
+
∂
ϕ
β
(
z
,
u
,
t
-
τ
)
∂
u
|
u
=
t
)
=
∑
α
‍
R
α
0
(
z
α
-
1
)
Φ
1
+
∑
β
‍
∂
Φ
1
∂
ϕ
β
∑
γ
‍
M
γ
β
b
γ
(
z
,
t
)
=
∑
α
‍
R
α
0
(
z
α
-
1
)
Φ
1
+
∑
β
‍
∂
Φ
1
∂
ϕ
β
∑
γ
‍
∂
ϕ
β
∂
z
γ
b
γ
(
z
,
t
)
=
∑
α
‍
R
α
0
(
z
α
-
1
)
Φ
1
+
∑
γ
‍
∂
Φ
1
∂
z
γ
b
γ
(
z
,
t
)
,
after retracing our steps from the previous Lemma.The result presented in Lemma9 was constructed using the method of variation of constants, which provides an alternative (constructive) proof.


### 3.2.  Basic Properties of the Solution

Definition 10.
A nonzero vectorv such that ∑
α
‍
δ
j
α
v
α
=
0, for j
=
1
,
…
,
N, is called a structural zero of F. Other zero eigenvalues of F are called nonstructural.Lemma 11.
(a) There are structural zeroes if and only ifE
>
N. (b) The only structurally stable zero eigenvalues of the matrix F are the structural zeroes.Proof.
(a) The “if” part is immediate sincedim
⁡
(
ker
⁡
(
δ
)
)
≥
E
-
N. The “only if” part follows from the assumption of independent population increments (and therefore independent populations). For (b) put F
=
r
δ. If F
w
=
0 for w nonzero and δ
w
=
p
≠
0, then r
p
=
0. This relation is not structurally stable since by adding c
≠
0 arbitrarily small to all diagonal elements of r we can assure r
′
p
≠
0 for this modified r
′ and for any p
≠
0.Lemma 12.
First integrals of motion for the characteristic equations. Letv be a structural zero of F, then
(22)
∏
α
‍
(
x
α
(
t
)
)
v
α
=
C
is a constant of motion for the characteristic equations (12).Proof.
Rewrite the first equation (12) as (indices renamed)
(23)
d
d
t
ln
⁡
x
α
=
∑
β
‍
F
β
α
(
x
β
-
1
)
.
Then,
(24)
∑
α
‍
v
α
d
d
t
ln
⁡
x
α
=
∑
β
‍
∑
α
‍
(
v
α
F
β
α
)
(
x
β
-
1
)
=
0
.
Since ∑
α
‍
(
v
α
F
β
α
)
=
0, the expression in the lemma is the exponential of the integral obtained from
(25)
∑
α
‍
v
α
d
d
t
ln
⁡
x
α
=
0
.Given relationship (1) between initial conditions, events, and populations, from the probability generating function in event space,
(26)
Ψ
(
z
1
,
…
,
z
E
,
s
,
t
)
=
∑
{
n
α
}
‍
(
∏
α
‍
z
α
n
α
)
P
(
n
1
⋯
n
E
,
s
,
t
)
,
a corresponding generating function in population space can be computed, namely,
(27)
Ψ
¯
(
y
1
,
…
,
y
N
,
s
,
t
)
=
∑
{
X
k
}
‍
(
∏
k
‍
y
k
X
k
)
P
(
X
1
,
…
,
X
N
,
s
,
t
)
,
letting {
n
α
}, {
X
k
} be shorthand for (
n
1
⋯
n
E
), (
X
1
,
…
,
X
N
).Lemma 13 (projection lemma).
Equationz
α
=
∏
k
=
1
N
‍
y
k
δ
k
α realises the transformation from event space to population space. Moreover,
(28)
Ψ
¯
(
y
1
,
…
,
y
N
,
s
,
t
)
=
(
∏
k
‍
y
k
X
k
0
)
Ψ
(
z
1
,
…
,
z
E
,
s
,
t
)
=
(
∏
k
‍
y
k
X
k
0
)
Ψ
(
∏
k
=
1
N
‍
y
k
δ
k
1
,
…
,
∏
k
=
1
N
‍
y
k
δ
k
E
,
s
,
t
)
.Proof.
Taking logarithms onz and y it is clear that whenever there exist structural zeroes the transformation is noninvertible. In this sense we call it a projection. Concerning the statement, we have
(29)
(
∏
k
‍
y
k
X
k
0
)
Ψ
(
z
1
,
…
,
z
E
,
s
,
t
)
=
(
∏
k
‍
y
k
X
k
0
)
∑
{
n
α
}
‍
(
∏
α
‍
(
∏
k
=
1
N
‍
y
k
δ
k
α
)
n
α
)
P
(
{
n
α
}
,
s
,
t
)
=
∑
{
n
α
}
‍
(
∏
k
=
1
N
‍
y
k
X
k
0
+
∑
α
δ
k
α
n
α
)
P
(
{
n
α
}
,
s
,
t
)
.
The statement is proved using (1) and rearranging the sums as
(30)
P
(
X
1
,
…
,
X
N
,
s
,
t
)
=
∑
{
n
α
}
‍
′
P
(
{
n
α
}
,
s
,
t
)
,
where the prime denotes sum over {
n
α
} such that X
j
0
+
∑
α
‍
δ
j
α
n
α
=
X
j.Corollary 14.
The integrals of motion of Lemma12 project under the projection Lemma to C
=
1.Proof.
The above substitution operated on∏
‍
z
α
v
α
=
C gives
(31)
C
=
∏
k
=
1
N
‍
y
k
(
∑
α
‍
v
α
δ
k
α
)
≡
∏
k
=
1
N
‍
y
k
0
=
1
.


### 3.3. Reduced Equation

Assumption 15.
There exist integersm
β, β
=
1
,
…
,
E such that (cf. (9))
(32)
R
α
0
=
∑
β
‍
F
α
β
m
β
.This assumption amounts to considering that there exists a way to mimic the initial condition on the rates with the same matrixF involved in the time evolution.Theorem 16.
Under Assumptions2, 4, and 15, equation (10) with natural initial condition Ψ
(
z
1
,
…
,
z
E
,
s
,
s
,
)
=
Φ
0
(
z
1
,
…
,
z
E
)
=
1 has the solution
(33)
Ψ
(
z
1
,
…
,
z
E
,
s
,
t
)
=
∏
β
‍
f
β
m
β
,
where z
β
f
β
(
z
,
s
,
t
)
=
ϕ
β
(
z
,
t
,
t
-
s
) in terms of the solutions of (12) of Lemma 8. Hence, f
(
z
,
s
,
t
) satisfies (after eliminating the auxiliary time τ)
(34)
d
log
⁡
(
f
β
(
z
,
s
,
t
)
)
d
s
=
-
∑
α
‍
F
α
β
(
s
)
(
z
α
f
α
(
z
,
s
,
t
)
-
1
)
,
f
β
(
z
,
t
,
t
)
=
1
.Proof.
The differences betweenf
(
) and ϕ
(
) are (a) factoring out the initial condition ϕ
β
(
z
,
t
,
0
)
=
z
β and (b) rearranging the time arguments since the argument in ϕ
(
) is natural to the ad hoc autonomous ODE, (12), while the argument in f
(
) is natural to the (nonautonomous) PDE, (10).
The proof is just a matter of rewriting the previous equations under the present assumptions. Because of the natural initial condition, the desired solution to (10) is given by (20) in Lemma 9 which can be rewritten using Assumption 15, as (we use z as shorthand for z
1
,
…
,
z
E)
(35)
Ψ
(
z
,
s
,
t
)
=
Φ
1
(
z
,
s
,
t
)
=
∏
β
‍
(
exp
⁡
∫
s
t
‍
∑
α
‍
F
α
β
(
τ
)
(
ϕ
α
(
z
,
t
,
t
-
τ
)
-
1
)
d
τ
)
m
β
.
Recalling that we defined z
α
f
α
(
z
,
s
,
t
)
=
ϕ
(
z
,
t
,
t
-
s
), the parenthesis can be recognised as the formal solution of (34) integrated from t to s.Corollary 17.
Forv a structural zero of Lemma 12, the generating function is unmodified upon the substitution m
α
→
m
α
+
k
v
α (k
∈
ℝ).Remark 18.
Having structural zeroes, it is a matter of choice to address the problem (a) using the variablesf
α and imposing the additional constraints given by structural zeroes or (b) shifting to population coordinates where these constraints are automatically incorporated. Because of Assumption 2 and Lemma 12, the transformation mapping one description to the other corresponds to
(36)
f
α
=
∏
j
‍
w
j
δ
j
α
.
However, option (b) proves to be more efficient when it comes to deal with approximate solutions (see Section 5).


## 3.1. The Method of Characteristics Revisited

Lemma 8.
Equation (10) with R
α
0
=
0 and with the initial condition Ψ
(
z
1
,
…
,
z
E
,
s
,
s
)
=
Φ
0
(
z
1
,
…
,
z
E
) admits the solution
(11)
Ψ
(
z
1
,
…
,
z
E
,
s
,
t
)
=
Φ
0
(
ϕ
(
z
,
t
,
t
-
s
)
)
,
where ϕ
(
z
,
t
,
τ
) is the flow of the following ODE with initial condition x
(
0
)
=
z, s
(
0
)
=
t:
(12)
d
x
β
d
τ
=
∑
α
‍
F
α
β
(
x
α
-
1
)
x
β
≡
b
β
(
x
,
s
)
,
d
s
d
τ
=
-
1
.Proof.
It is convenient to renames as s
=
x
E
+
1, introduce b
E
+
1
=
-
1, and recall that b
β depends only on the extended variable array x. However, when necessary, we will write explicitly x
(
0
) as (
z
,
t
). We have that
(13)
∂
Ψ
∂
t
|
s
=
∑
β
E
+
1
‍
∂
Φ
0
∂
x
β
|
ϕ
(
x
,
t
-
s
)
d
ϕ
β
(
z
,
t
,
t
-
s
)
d
t
=
∑
β
E
+
1
‍
∂
Φ
0
∂
x
β
|
ϕ
(
x
,
t
-
s
)
×
(
b
β
(
ϕ
(
z
,
t
,
t
-
s
)
)
+
∂
ϕ
β
(
z
,
u
,
t
-
s
)
∂
u
|
u
=
t
)
.
Further, consider the monodromy matrix [40]
(14)
M
γ
β
=
∂
ϕ
β
(
x
,
τ
)
∂
x
γ
.
M transforms the initial velocity into the final velocity; hence
(15)
∑
γ
‍
M
γ
β
b
γ
(
z
,
t
)
=
b
β
(
ϕ
(
z
,
t
,
τ
)
)
.
Thus, we get
(16)
∂
Ψ
∂
t
|
s
=
∑
β
E
+
1
‍
∂
Φ
0
∂
x
β
|
ϕ
(
x
,
t
-
s
)
(
∑
γ
E
+
1
‍
M
γ
β
b
γ
(
z
,
t
)
+
∂
ϕ
β
(
z
,
u
,
t
-
s
)
∂
u
|
u
=
t
)
=
∑
β
E
+
1
‍
∂
Φ
0
∂
x
β
|
ϕ
(
x
,
t
-
s
)
×
(
∑
γ
E
‍
M
γ
β
b
γ
(
z
,
t
)
-
M
E
+
1
β
+
∂
ϕ
β
(
z
,
u
,
t
-
s
)
∂
u
|
u
=
t
)
=
∑
γ
E
‍
∂
Ψ
∂
z
γ
b
γ
(
z
,
t
)
.
From the definition of monodromy matrix above it follows that
(17)
-
M
E
+
1
β
+
∂
ϕ
β
(
z
,
u
,
t
-
s
)
∂
u
|
u
=
t
=
0
,
and by (11) and the chain rule,
(18)
∑
β
‍
∂
Φ
0
∂
x
β
|
x
=
ϕ
(
x
,
t
-
s
)
M
γ
β
=
∑
β
‍
∂
Φ
0
∂
x
β
|
x
=
ϕ
(
x
,
t
-
s
)
∂
ϕ
β
(
z
,
t
,
t
-
s
)
∂
z
γ
=
∂
Ψ
∂
z
γ
.Lemma 9.
Equation (10) with R
α
0
≠
0 and with the initial condition Ψ
(
z
1
,
…
,
z
E
,
s
,
s
)
=
Φ
0
(
z
1
,
…
,
z
E
) admits the solution
(19)
Ψ
(
z
1
,
…
,
z
E
,
s
,
t
)
=
Φ
0
(
ϕ
(
z
,
t
,
t
-
s
)
)
Φ
1
(
z
1
,
…
,
z
E
,
s
,
t
)
,
where Φ
0
(
ϕ
(
z
,
t
,
t
-
s
)
) is a solution of the homogeneous problem with R
α
0
=
0 and
(20)
Φ
1
(
z
1
,
…
,
z
E
,
s
,
t
)
=
exp
⁡
∫
s
t
‍
(
∑
α
‍
R
α
0
(
τ
)
(
ϕ
α
(
z
,
t
,
t
-
τ
)
-
1
)
d
τ
)
.Proof.
The proof is nothing more than an application of the method of the variation of the constants. The details are as follows.
Introducing the proposed solutionΨ into (10) and using Lemma 8 to eliminate Φ
0 we obtain that Φ
1 satisfies (10) with initial condition Φ
1
(
z
,
s
,
s
)
=
1. Computing directly from (20) we have that
(21)
∂
Φ
1
∂
t
|
s
=
∑
α
‍
R
α
0
(
ϕ
α
(
z
,
t
,
0
)
-
1
)
Φ
1
+
∑
β
‍
∂
Φ
1
∂
ϕ
β
∂
ϕ
β
(
z
,
t
,
t
-
τ
)
∂
t
=
∑
α
‍
R
α
0
(
z
α
-
1
)
Φ
1
+
∑
β
‍
∂
Φ
1
∂
ϕ
β
(
b
β
(
ϕ
(
z
,
t
,
t
-
τ
)
)
+
∂
ϕ
β
(
z
,
u
,
t
-
τ
)
∂
u
|
u
=
t
)
=
∑
α
‍
R
α
0
(
z
α
-
1
)
Φ
1
+
∑
β
‍
∂
Φ
1
∂
ϕ
β
∑
γ
‍
M
γ
β
b
γ
(
z
,
t
)
=
∑
α
‍
R
α
0
(
z
α
-
1
)
Φ
1
+
∑
β
‍
∂
Φ
1
∂
ϕ
β
∑
γ
‍
∂
ϕ
β
∂
z
γ
b
γ
(
z
,
t
)
=
∑
α
‍
R
α
0
(
z
α
-
1
)
Φ
1
+
∑
γ
‍
∂
Φ
1
∂
z
γ
b
γ
(
z
,
t
)
,
after retracing our steps from the previous Lemma.The result presented in Lemma9 was constructed using the method of variation of constants, which provides an alternative (constructive) proof.


## 3.2.  Basic Properties of the Solution

Definition 10.
A nonzero vectorv such that ∑
α
‍
δ
j
α
v
α
=
0, for j
=
1
,
…
,
N, is called a structural zero of F. Other zero eigenvalues of F are called nonstructural.Lemma 11.
(a) There are structural zeroes if and only ifE
>
N. (b) The only structurally stable zero eigenvalues of the matrix F are the structural zeroes.Proof.
(a) The “if” part is immediate sincedim
⁡
(
ker
⁡
(
δ
)
)
≥
E
-
N. The “only if” part follows from the assumption of independent population increments (and therefore independent populations). For (b) put F
=
r
δ. If F
w
=
0 for w nonzero and δ
w
=
p
≠
0, then r
p
=
0. This relation is not structurally stable since by adding c
≠
0 arbitrarily small to all diagonal elements of r we can assure r
′
p
≠
0 for this modified r
′ and for any p
≠
0.Lemma 12.
First integrals of motion for the characteristic equations. Letv be a structural zero of F, then
(22)
∏
α
‍
(
x
α
(
t
)
)
v
α
=
C
is a constant of motion for the characteristic equations (12).Proof.
Rewrite the first equation (12) as (indices renamed)
(23)
d
d
t
ln
⁡
x
α
=
∑
β
‍
F
β
α
(
x
β
-
1
)
.
Then,
(24)
∑
α
‍
v
α
d
d
t
ln
⁡
x
α
=
∑
β
‍
∑
α
‍
(
v
α
F
β
α
)
(
x
β
-
1
)
=
0
.
Since ∑
α
‍
(
v
α
F
β
α
)
=
0, the expression in the lemma is the exponential of the integral obtained from
(25)
∑
α
‍
v
α
d
d
t
ln
⁡
x
α
=
0
.Given relationship (1) between initial conditions, events, and populations, from the probability generating function in event space,
(26)
Ψ
(
z
1
,
…
,
z
E
,
s
,
t
)
=
∑
{
n
α
}
‍
(
∏
α
‍
z
α
n
α
)
P
(
n
1
⋯
n
E
,
s
,
t
)
,
a corresponding generating function in population space can be computed, namely,
(27)
Ψ
¯
(
y
1
,
…
,
y
N
,
s
,
t
)
=
∑
{
X
k
}
‍
(
∏
k
‍
y
k
X
k
)
P
(
X
1
,
…
,
X
N
,
s
,
t
)
,
letting {
n
α
}, {
X
k
} be shorthand for (
n
1
⋯
n
E
), (
X
1
,
…
,
X
N
).Lemma 13 (projection lemma).
Equationz
α
=
∏
k
=
1
N
‍
y
k
δ
k
α realises the transformation from event space to population space. Moreover,
(28)
Ψ
¯
(
y
1
,
…
,
y
N
,
s
,
t
)
=
(
∏
k
‍
y
k
X
k
0
)
Ψ
(
z
1
,
…
,
z
E
,
s
,
t
)
=
(
∏
k
‍
y
k
X
k
0
)
Ψ
(
∏
k
=
1
N
‍
y
k
δ
k
1
,
…
,
∏
k
=
1
N
‍
y
k
δ
k
E
,
s
,
t
)
.Proof.
Taking logarithms onz and y it is clear that whenever there exist structural zeroes the transformation is noninvertible. In this sense we call it a projection. Concerning the statement, we have
(29)
(
∏
k
‍
y
k
X
k
0
)
Ψ
(
z
1
,
…
,
z
E
,
s
,
t
)
=
(
∏
k
‍
y
k
X
k
0
)
∑
{
n
α
}
‍
(
∏
α
‍
(
∏
k
=
1
N
‍
y
k
δ
k
α
)
n
α
)
P
(
{
n
α
}
,
s
,
t
)
=
∑
{
n
α
}
‍
(
∏
k
=
1
N
‍
y
k
X
k
0
+
∑
α
δ
k
α
n
α
)
P
(
{
n
α
}
,
s
,
t
)
.
The statement is proved using (1) and rearranging the sums as
(30)
P
(
X
1
,
…
,
X
N
,
s
,
t
)
=
∑
{
n
α
}
‍
′
P
(
{
n
α
}
,
s
,
t
)
,
where the prime denotes sum over {
n
α
} such that X
j
0
+
∑
α
‍
δ
j
α
n
α
=
X
j.Corollary 14.
The integrals of motion of Lemma12 project under the projection Lemma to C
=
1.Proof.
The above substitution operated on∏
‍
z
α
v
α
=
C gives
(31)
C
=
∏
k
=
1
N
‍
y
k
(
∑
α
‍
v
α
δ
k
α
)
≡
∏
k
=
1
N
‍
y
k
0
=
1
.


## 3.3. Reduced Equation

Assumption 15.
There exist integersm
β, β
=
1
,
…
,
E such that (cf. (9))
(32)
R
α
0
=
∑
β
‍
F
α
β
m
β
.This assumption amounts to considering that there exists a way to mimic the initial condition on the rates with the same matrixF involved in the time evolution.Theorem 16.
Under Assumptions2, 4, and 15, equation (10) with natural initial condition Ψ
(
z
1
,
…
,
z
E
,
s
,
s
,
)
=
Φ
0
(
z
1
,
…
,
z
E
)
=
1 has the solution
(33)
Ψ
(
z
1
,
…
,
z
E
,
s
,
t
)
=
∏
β
‍
f
β
m
β
,
where z
β
f
β
(
z
,
s
,
t
)
=
ϕ
β
(
z
,
t
,
t
-
s
) in terms of the solutions of (12) of Lemma 8. Hence, f
(
z
,
s
,
t
) satisfies (after eliminating the auxiliary time τ)
(34)
d
log
⁡
(
f
β
(
z
,
s
,
t
)
)
d
s
=
-
∑
α
‍
F
α
β
(
s
)
(
z
α
f
α
(
z
,
s
,
t
)
-
1
)
,
f
β
(
z
,
t
,
t
)
=
1
.Proof.
The differences betweenf
(
) and ϕ
(
) are (a) factoring out the initial condition ϕ
β
(
z
,
t
,
0
)
=
z
β and (b) rearranging the time arguments since the argument in ϕ
(
) is natural to the ad hoc autonomous ODE, (12), while the argument in f
(
) is natural to the (nonautonomous) PDE, (10).
The proof is just a matter of rewriting the previous equations under the present assumptions. Because of the natural initial condition, the desired solution to (10) is given by (20) in Lemma 9 which can be rewritten using Assumption 15, as (we use z as shorthand for z
1
,
…
,
z
E)
(35)
Ψ
(
z
,
s
,
t
)
=
Φ
1
(
z
,
s
,
t
)
=
∏
β
‍
(
exp
⁡
∫
s
t
‍
∑
α
‍
F
α
β
(
τ
)
(
ϕ
α
(
z
,
t
,
t
-
τ
)
-
1
)
d
τ
)
m
β
.
Recalling that we defined z
α
f
α
(
z
,
s
,
t
)
=
ϕ
(
z
,
t
,
t
-
s
), the parenthesis can be recognised as the formal solution of (34) integrated from t to s.Corollary 17.
Forv a structural zero of Lemma 12, the generating function is unmodified upon the substitution m
α
→
m
α
+
k
v
α (k
∈
ℝ).Remark 18.
Having structural zeroes, it is a matter of choice to address the problem (a) using the variablesf
α and imposing the additional constraints given by structural zeroes or (b) shifting to population coordinates where these constraints are automatically incorporated. Because of Assumption 2 and Lemma 12, the transformation mapping one description to the other corresponds to
(36)
f
α
=
∏
j
‍
w
j
δ
j
α
.
However, option (b) proves to be more efficient when it comes to deal with approximate solutions (see Section 5).


## 4. Examples of Exact Solutions

### 4.1. Example I: Pure Death Process

In the caseN
=
E
=
1, δ
=
-
1, we have r
>
0, F
=
-
r
(
s
), and m
=
-
M
<
0 (M is the initial population value). Equation (34) reads d
f
/
d
s
=
f
r
(
z
f
-
1
), to be integrated in [
t
,
s
], with f
(
z
,
t
,
t
)
=
1. Hence, for constant death rate r we obtain f
(
z
,
s
,
t
)
=
[
z
+
(
1
-
z
)
e
-
r
(
t
-
s
)
]
-
1 and
(37)
Ψ
(
z
,
s
,
t
)
=
[
e
-
r
(
t
-
s
)
+
z
(
1
-
e
-
r
(
t
-
s
)
)
]
M
.


### 4.2.  Example II: Linear Birth and Death Process

We haveN
=
1, E
=
2, and δ
1
1
=
-
δ
1
2
=
1, while r
1 denotes the birth rate and r
2 denotes the death rate. X
0
=
m
1
-
m
2 is the initial population. There exists one structural zero, with vector v
=
(
1,1
)
T, while the corresponding constant of motion results in the relations z
1
z
2
=
C and f
1
f
2
=
1. The differential equation (34) for f
1 (to be integrated in [
t
,
s
]) becomes
(38)
d
f
1
d
s
=
(
r
1
z
1
f
1
2
-
(
r
1
+
r
2
)
f
1
+
C
r
2
z
1
)
,
f
1
(
z
,
t
,
t
)
=
1
,
with solution (for C
=
1 which will be the case after projection and assuming time-independent birth rates)
(39)
f
1
(
z
,
s
,
t
)
=
1
z
1
(
r
2
-
r
1
z
1
)
W
-
r
2
(
1
-
z
1
)
(
r
2
-
r
1
z
1
)
W
-
r
1
(
1
-
z
1
)
,
where W
=
exp
⁡
(
r
2
-
r
1
)
(
t
-
s
). Finally, Ψ
(
z
,
s
,
t
)
=
(
f
1
(
z
,
s
,
t
)
)
m
1
-
m
2 which after projecting in population space using Lemma 13 and Corollary 14 gives (cf. [41])
(40)
Ψ
¯
(
y
,
s
,
t
)
=
[
(
r
2
-
r
1
y
)
W
-
r
2
(
1
-
y
)
(
r
2
-
r
1
y
)
W
-
r
1
(
1
-
y
)
]
m
1
-
m
2
.Remark 19.
A pure death linear process corresponds to a binomial distribution, a pure birth process corresponds to a negative binomial, and a birth-death process corresponds to a combination of both as independent processes.


## 4.1. Example I: Pure Death Process

In the caseN
=
E
=
1, δ
=
-
1, we have r
>
0, F
=
-
r
(
s
), and m
=
-
M
<
0 (M is the initial population value). Equation (34) reads d
f
/
d
s
=
f
r
(
z
f
-
1
), to be integrated in [
t
,
s
], with f
(
z
,
t
,
t
)
=
1. Hence, for constant death rate r we obtain f
(
z
,
s
,
t
)
=
[
z
+
(
1
-
z
)
e
-
r
(
t
-
s
)
]
-
1 and
(37)
Ψ
(
z
,
s
,
t
)
=
[
e
-
r
(
t
-
s
)
+
z
(
1
-
e
-
r
(
t
-
s
)
)
]
M
.


## 4.2.  Example II: Linear Birth and Death Process

We haveN
=
1, E
=
2, and δ
1
1
=
-
δ
1
2
=
1, while r
1 denotes the birth rate and r
2 denotes the death rate. X
0
=
m
1
-
m
2 is the initial population. There exists one structural zero, with vector v
=
(
1,1
)
T, while the corresponding constant of motion results in the relations z
1
z
2
=
C and f
1
f
2
=
1. The differential equation (34) for f
1 (to be integrated in [
t
,
s
]) becomes
(38)
d
f
1
d
s
=
(
r
1
z
1
f
1
2
-
(
r
1
+
r
2
)
f
1
+
C
r
2
z
1
)
,
f
1
(
z
,
t
,
t
)
=
1
,
with solution (for C
=
1 which will be the case after projection and assuming time-independent birth rates)
(39)
f
1
(
z
,
s
,
t
)
=
1
z
1
(
r
2
-
r
1
z
1
)
W
-
r
2
(
1
-
z
1
)
(
r
2
-
r
1
z
1
)
W
-
r
1
(
1
-
z
1
)
,
where W
=
exp
⁡
(
r
2
-
r
1
)
(
t
-
s
). Finally, Ψ
(
z
,
s
,
t
)
=
(
f
1
(
z
,
s
,
t
)
)
m
1
-
m
2 which after projecting in population space using Lemma 13 and Corollary 14 gives (cf. [41])
(40)
Ψ
¯
(
y
,
s
,
t
)
=
[
(
r
2
-
r
1
y
)
W
-
r
2
(
1
-
y
)
(
r
2
-
r
1
y
)
W
-
r
1
(
1
-
y
)
]
m
1
-
m
2
.Remark 19.
A pure death linear process corresponds to a binomial distribution, a pure birth process corresponds to a negative binomial, and a birth-death process corresponds to a combination of both as independent processes.


## 5. Consistent Approximations

When (34) of Theorem 16 cannot be solved exactly, we have to rely on approximate solutions. It is desirable to work with consistent approximations, namely, those where basic properties of the problem are fulfilled exactly, rather than “up to an error of size ϵ”. For example, since our solutions should be probabilities, we want the coefficients of Ψ
(
z
1
,
…
,
z
E
,
t
) to be always nonnegative and sum up to one. Ideally, we want approximations to be constrained to satisfy Lemma 12 and to preserve positive invariance in population space (no jumps into negative populations). Whenever consistency is not automatic, it is mandatory to analyse the approximate solution before jumping into conclusions.


### 5.1. Poisson Approximation

Let us replacef
β
≡
1 in the RHS of (34) (this amounts to propose that f
β is constant and equal to its initial condition). The solution of the modified equation will be a good approximation to the exact solution for sufficiently short times such that f is not significantly different from one.For constant ratesr
α the rhs of (34) does not depend on time. The approximate solution reads f
α
(
z
,
s
,
t
)
=
exp
⁡
(
(
t
-
s
)
∑
β
‍
F
β
α
(
z
β
-
1
)
) and
(41)
Ψ
(
z
1
,
…
,
z
E
,
s
,
t
)
=
∏
α
E
‍
e
m
α
(
t
-
s
)
∑
β
F
β
α
(
z
β
-
1
)
=
e
(
t
-
s
)
∑
β
R
β
0
(
z
β
-
1
)
.In other words, the generating function is approximated as a product of individual Poisson generating functions.This approximation has an error ofo
(
t
). Simple as it looks, this approach gives a good probability generating function that respects structural zeroes, but it does not guarantee positive invariance: a Poisson jump could throw the system into negative populations. However, the probability of such a jump is also o
(
t
). Fortunately, there is a natural way to impose positive invariance and also a natural way to improve the error up to o
(
t
)
3
/
2 in the broader framework of general transition rates [6]. In this way, the approximation can be safely used, with error control, to address problems that resist exact computation one way or the other.


### 5.2.  Consistent Multinomial Approximations

A systematic approach to obtain a chain of approximations of increasing accuracy to (34) can be devised via Picard iterations. Firstly, it is convenient to recast the problem in population coordinates (cf. Remark 18), so that the structural zeroes are automatically taken into account, regardless of other properties of the adopted approximation. From f
α
=
∏
j
‍
w
j
δ
j
α we obtain f
˙
α
/
f
α
=
∑
k
‍
δ
k
α
w
˙
k
/
w
k and (34) can be recast as
(42)
∑
k
‍
δ
k
α
(
w
˙
k
w
k
+
∑
β
‍
r
β
k
(
z
β
(
∏
j
‍
w
j
δ
j
β
)
-
1
)
)
=
0
,
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
w
k
(
z
,
t
,
t
)
=
1
.
Because of the linear independence of population increments (and hence of populations), Lemma 3(b), the counterpart of (34) in population coordinates reads
(43)
w
˙
k
=
-
w
k
∑
β
‍
r
β
k
(
z
β
(
∏
j
‍
w
j
δ
j
β
)
-
1
)
,
w
k
(
z
,
t
,
t
)
=
1
.Remark 20.
In the case of linearlydependent population increments, (34) still holds, while it is possible to elaborate a more general form in terms of the w
j instead of (43). This is, however, outside the scope of this paper.Since∑
α
‍
δ
j
α
m
α
=
X
j
(
s
), we obtain Ψ
(
z
1
,
…
,
z
E
,
s
,
t
)
=
∏
j
‍
w
j
X
j
(
s
), which is, so far, an exact result. Let us produce some approximated values for the w
k’s, and correspondingly for Ψ
(
z
1
,
…
,
z
E
,
s
,
t
).To avoid further notational complications, we will assume in the sequel thatr
β
k is time independent. Hence, w
k will depend only on the time difference t
-
s. Equation (43) can be recast as
(44)
w
˙
k
=
w
k
∑
β
‍
r
β
k
(
z
β
(
∏
j
‍
w
j
δ
j
β
)
-
1
)
,
w
k
(
z
,
0
)
=
1
,
now to be integrated along the interval [
0
,
t
-
s
]. We will even set s
=
0 in what follows.Remark 21.
For the case of time-independent proportionality factorsr
β
k, a corresponding modification can be introduced in (34), namely, overall change of sign, integration in [
0
,
t
-
s
], and setting s
=
0 to reduce the notational burden.Equations (44) satisfy the Picard-Lindelöf theorem since the RHS is a Lipschitz function. Then the Picard iteration scheme converges to the (unique) solution. Taking advantage of the initial condition, the nth order truncated version of the Picard iterations reads
(45)
w
k
(
0
)
(
h
)
=
1
,
w
k
(
1
)
(
h
)
=
1
+
∫
0
h
‍
w
k
(
0
)
(
t
)
∑
β
‍
r
β
k
(
z
β
∏
j
‍
(
w
j
(
0
)
(
t
)
)
δ
j
β
-
1
)
d
t
w
k
(
1
)
(
h
)
=
1
+
h
∑
β
‍
r
β
k
(
z
β
-
1
)
,
w
k
(
n
)
(
h
)
=
1
+
∫
0
h
‍
d
t

[
(
z
β
∏
j
‍
(
w
j
(
n
-
1
)
(
t
)
)
δ
j
β
-
1
)
w
k
(
n
-
1
)
(
t
)
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
×
∑
β
‍
r
β
k
(
z
β
∏
j
‍
(
w
j
(
n
-
1
)
(
t
)
)
δ
j
β
-
1
)
]
n
-
1
,
where [
⋯
]
m denote the Maclaurin polynomial of order m in t.Lemma 22.
Fort
∈
[
0
,
h
], h
>
0 being sufficiently small, and for each order of approximation n
≥
0, w
k
(
n
)
(
t
) in (44) is a polynomial generating function for one individual of the population; that is, it has the following properties:(a)
w
k
(
n
)
(
t
) is a polynomial of degree n in z
1
,
…
,
z
E, where the coefficient of z
1
n
1
,
…
,
z
E
n
E is O
(
t
p
), with p
=
n
1
+
⋯
+
n
E;
(b)
w
k
(
n
)
(
1
,
…
,
1
,
t
)
=
1;
(c)
Δ
k
(
n
)
(
h
)
=
w
k
(
n
)
(
h
)
-
w
k
(
n
-
1
)
(
h
)
=
O
(
h
n
);
(d)
the coefficients ofw
k
(
n
)
(
t
) regarded as a polynomial in {
z
α
} are nonnegative functions of time.Proof.
(a) The statement holds forn
=
0 and n
=
1. Assuming that it holds up to order n
-
1, we have that [
w
k
(
n
-
1
)
(
t
)
∏
j
‍
(
w
j
(
n
-
1
)
(
t
)
)
δ
j
β
]
n
-
1 is also a polynomial of degree n
-
1 in z
1
,
…
,
z
E with time-dependent coefficients of type O
(
t
p
) since each power p of t carries along a power p in {
z
α
} (as well as lower z-powers). A Picard iteration step gives
(46)
w
k
(
n
)
(
t
)
=
1
+
∑
β
‍
r
β
k
z
β
∫
0
h
‍
[
w
k
(
n
-
1
)
(
t
)
∏
j
‍
(
w
j
(
n
-
1
)
(
t
)
)
δ
j
β
]
n
-
1
d
t
-
(
∑
β
‍
r
β
k
)
∫
0
h
‍
w
k
(
n
-
1
)
(
t
)
d
t
.
The first sum contributes with one time integration and one z-power, which exactly preserves the polynomial structure and the order of the coefficients. The result is a polynomial of degree n in z
1
,
…
,
z
E, where the coefficients are polynomials of highest degree at most n in h, each one of lowest order O
(
h
p
) (as in the statement). The second integral does not alter these properties since it generates again a polynomial of degree n
-
1 in z
1
,
…
,
z
E, with time-dependent coefficients of order O
(
h
p
+
1
) or higher, up to h
n (they are just the coefficients of w
k
(
n
-
1
)
(
t
) integrated once in time).
(b) Lettingz
1
=
z
2
=
⋯
=
z
E
=
1, Picard iteration gives trivially w
k
(
n
)
(
t
)
=
1 to all orders, again by induction.
(c) The statement holds forn
=
1. Assuming it holds up to n
-
1, we compute
(47)
w
k
(
n
)
(
h
)
-
w
k
(
n
-
1
)
(
h
)
=
∑
β
‍
r
β
k
z
β
×
∫
0
h
‍
(
[
w
k
(
n
-
2
)
(
t
)
∏
j
‍
(
w
j
(
n
-
2
)
(
t
)
)
δ
j
β
]
n
-
2
[
(
w
j
(
n
-
2
)
(
t
)
+
Δ
j
(
n
-
1
)
(
t
)
)
δ
j
β
(
w
k
(
n
-
2
)
(
t
)
+
Δ
k
(
n
-
1
)
(
t
)
)
×
∏
j
‍
(
w
j
(
n
-
2
)
(
t
)
+
Δ
j
(
n
-
1
)
(
t
)
)
δ
j
β
]
n
-
1
-
[
w
k
(
n
-
2
)
(
t
)
∏
j
‍
(
w
j
(
n
-
2
)
(
t
)
)
δ
j
β
]
n
-
2
)
d
t
-
(
∑
β
‍
r
β
k
)
∫
0
h
‍
(
w
k
(
n
-
1
)
(
t
)
-
w
k
(
n
-
2
)
(
t
)
)
d
t
.
Both expressions in the first integral have the same Maclaurin polynomial up to order (
n
-
2
). Hence, the difference in the first integral contains only terms of order O
(
t
n
-
1
). This is also the case for the second integral. After integration, we obtain Δ
k
(
n
)
(
h
)
=
O
(
h
n
).
(d) The statement holds forn
=
0 and n
=
1 and in general for the zero-order coefficient, because it is unity for t
=
0, as a consequence of the initial condition. For the higher-order coefficients, we use induction again. If all z-coefficients in w
k
(
n
-
1
)
(
t
) are nonnegative, then the same holds automatically for w
k
(
n
)
(
h
) up to order n
-
1, since the O
(
h
n
) contribution adding to each inherited coefficient from the previous step does not alter its sign, for h sufficiently small. It remains to show that the coefficients corresponding to p
=
n are nonnegative. These contributions arise from the (
n
-
1
)th order of
(48)
r
β
k
z
β
[
w
k
(
n
-
1
)
(
t
)
∏
j
‍
(
w
j
(
n
-
1
)
(
t
)
)
δ
j
β
]
n
-
1
after time integration. If all δ
j
β are positive, then the expression above is a product of z-polynomials with nonnegative coefficients and the statement follows. Suppose that some δ
j
β
=
-
1. Then by Lemma 7, r
β
k
=
0 for k
≠
j. Hence,
(49)
r
β
k
z
β
w
k
(
n
-
1
)
(
t
)
∏
j
‍
(
w
j
(
n
-
1
)
(
t
)
)
δ
j
β
=
{
0
,
j
≠
k
,
r
β
k
z
β
∏
j
≠
k
‍
(
w
j
(
n
-
1
)
(
t
)
)
δ
j
β
,
j
=
k
.
However, δ
j
β
≠
-
1 for j
≠
k; that is, they are nonnegative, and hence the expression has only nonnegative coefficients.


#### 5.2.1. Basic Properties of the First- and Second-Order Approximations

Let us analyse the generating functionΨ
(
z
1
,
…
,
z
E
,
t
)
=
∏
j
‍
w
j
X
j
(
0
). The simplest nontrivial applications of the approximation scheme are the first- and second-order approximations. Expanding explicitly up to second-order, we get
(50)
w
k
(
h
)
=
1
+
h
∑
β
‍
r
β
k
(
z
β
-
1
)
+
h
2
2
(
∑
β
‍
r
β
k
(
z
β
-
1
)
)
×
(
∑
γ
‍
r
γ
k
(
z
γ
-
1
)
)
+
h
2
2
∑
β
‍
r
β
k
z
β
∑
γ
‍
F
γ
β
(
z
γ
-
1
)
.We begin by noticing that the first-order approximation is obtained by neglecting all terms inh
2 in the above equation. Whenever w
j is linear in z, which is always the case for the first-order approximation, the resulting generating function Ψ corresponds to the product of (independent) multinomial distributions for the events affecting each subpopulation.Let us now considerΨ
(
z
1
,
…
,
z
E
,
t
) computed with the second-order approximated w
j’s to perform one simulation step of size h.Counting powers ofz
γ on each w
k we have (a) the probability of no events occurring in time-step h
(51)
P
k
(
h
,
0
)
=
1
-
h
(
∑
β
‍
r
β
k
)
+
h
2
2
(
∑
β
‍
r
β
k
)
2
,
(b) the probability of only one (β) event occurring in time-step h
(52)
P
k
(
h
,
I
,
β
)
=
r
β
k
{
h
-
h
2
∑
α
‍
r
α
k
-
h
2
2
∑
α
‍
F
α
β
}
,
and (c) the probability of another event (α) occurring after the first one:
(53)
P
k
(
h
,
I
I
,
β
,
α
)
=
h
2
2
r
β
k
(
r
α
k
+
F
α
β
)
.


#### 5.2.2. Simulation Algorithm

A simulation step can be operated in the following way. Let(54)
∑
α
(
P
k
(
h
,
I
I
,
β
,
α
)
)
+
P
k
(
h
,
I
,
β
)
=
r
β
k
{
h
-
h
2
2
∑
α
‍
r
α
k
}
=
P
¯
k
(
h
,
I
,
β
)
be the probability of one (β) event occurring alone or followed by a second event. It follows that
(55)
P
k
(
h
,
I
I
,
β
,
α
)
=
P
¯
k
(
h
,
I
,
β
)
h
2
(
r
k
α
+
F
α
β
)
(
1
-
(
h
/
2
)
∑
α
r
α
k
)
+
o
(
h
2
)
.
The second factor above is the conditional probability P
k
(
h
,
α
/
β
) for the second event being α given that there was a first (β) event.The algorithm proceeds in two steps. First a multinomial deviate is computed with probabilities{
P
¯
k
(
h
,
I
,
β
)
} (the probability of no event occurring is still P
k
(
h
,
0
)). This gives the number n
β of occurrences of each event. In the second step, for each n
β we compute the deviates for the occurrence of a second event or no more events, with the multinomial given by the array of probabilities {
P
k
(
h
,
α
/
β
)
} defined above (the probability of no second event occurring is here 1
-
(
h
/
2
)
(
∑
α
‍
(
r
k
α
+
F
α
β
)
/
(
1
-
(
h
/
2
)
∑
α
‍
r
α
k
)
)
).


## 5.1. Poisson Approximation

Let us replacef
β
≡
1 in the RHS of (34) (this amounts to propose that f
β is constant and equal to its initial condition). The solution of the modified equation will be a good approximation to the exact solution for sufficiently short times such that f is not significantly different from one.For constant ratesr
α the rhs of (34) does not depend on time. The approximate solution reads f
α
(
z
,
s
,
t
)
=
exp
⁡
(
(
t
-
s
)
∑
β
‍
F
β
α
(
z
β
-
1
)
) and
(41)
Ψ
(
z
1
,
…
,
z
E
,
s
,
t
)
=
∏
α
E
‍
e
m
α
(
t
-
s
)
∑
β
F
β
α
(
z
β
-
1
)
=
e
(
t
-
s
)
∑
β
R
β
0
(
z
β
-
1
)
.In other words, the generating function is approximated as a product of individual Poisson generating functions.This approximation has an error ofo
(
t
). Simple as it looks, this approach gives a good probability generating function that respects structural zeroes, but it does not guarantee positive invariance: a Poisson jump could throw the system into negative populations. However, the probability of such a jump is also o
(
t
). Fortunately, there is a natural way to impose positive invariance and also a natural way to improve the error up to o
(
t
)
3
/
2 in the broader framework of general transition rates [6]. In this way, the approximation can be safely used, with error control, to address problems that resist exact computation one way or the other.


## 5.2.  Consistent Multinomial Approximations

A systematic approach to obtain a chain of approximations of increasing accuracy to (34) can be devised via Picard iterations. Firstly, it is convenient to recast the problem in population coordinates (cf. Remark 18), so that the structural zeroes are automatically taken into account, regardless of other properties of the adopted approximation. From f
α
=
∏
j
‍
w
j
δ
j
α we obtain f
˙
α
/
f
α
=
∑
k
‍
δ
k
α
w
˙
k
/
w
k and (34) can be recast as
(42)
∑
k
‍
δ
k
α
(
w
˙
k
w
k
+
∑
β
‍
r
β
k
(
z
β
(
∏
j
‍
w
j
δ
j
β
)
-
1
)
)
=
0
,
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
w
k
(
z
,
t
,
t
)
=
1
.
Because of the linear independence of population increments (and hence of populations), Lemma 3(b), the counterpart of (34) in population coordinates reads
(43)
w
˙
k
=
-
w
k
∑
β
‍
r
β
k
(
z
β
(
∏
j
‍
w
j
δ
j
β
)
-
1
)
,
w
k
(
z
,
t
,
t
)
=
1
.Remark 20.
In the case of linearlydependent population increments, (34) still holds, while it is possible to elaborate a more general form in terms of the w
j instead of (43). This is, however, outside the scope of this paper.Since∑
α
‍
δ
j
α
m
α
=
X
j
(
s
), we obtain Ψ
(
z
1
,
…
,
z
E
,
s
,
t
)
=
∏
j
‍
w
j
X
j
(
s
), which is, so far, an exact result. Let us produce some approximated values for the w
k’s, and correspondingly for Ψ
(
z
1
,
…
,
z
E
,
s
,
t
).To avoid further notational complications, we will assume in the sequel thatr
β
k is time independent. Hence, w
k will depend only on the time difference t
-
s. Equation (43) can be recast as
(44)
w
˙
k
=
w
k
∑
β
‍
r
β
k
(
z
β
(
∏
j
‍
w
j
δ
j
β
)
-
1
)
,
w
k
(
z
,
0
)
=
1
,
now to be integrated along the interval [
0
,
t
-
s
]. We will even set s
=
0 in what follows.Remark 21.
For the case of time-independent proportionality factorsr
β
k, a corresponding modification can be introduced in (34), namely, overall change of sign, integration in [
0
,
t
-
s
], and setting s
=
0 to reduce the notational burden.Equations (44) satisfy the Picard-Lindelöf theorem since the RHS is a Lipschitz function. Then the Picard iteration scheme converges to the (unique) solution. Taking advantage of the initial condition, the nth order truncated version of the Picard iterations reads
(45)
w
k
(
0
)
(
h
)
=
1
,
w
k
(
1
)
(
h
)
=
1
+
∫
0
h
‍
w
k
(
0
)
(
t
)
∑
β
‍
r
β
k
(
z
β
∏
j
‍
(
w
j
(
0
)
(
t
)
)
δ
j
β
-
1
)
d
t
w
k
(
1
)
(
h
)
=
1
+
h
∑
β
‍
r
β
k
(
z
β
-
1
)
,
w
k
(
n
)
(
h
)
=
1
+
∫
0
h
‍
d
t

[
(
z
β
∏
j
‍
(
w
j
(
n
-
1
)
(
t
)
)
δ
j
β
-
1
)
w
k
(
n
-
1
)
(
t
)
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
×
∑
β
‍
r
β
k
(
z
β
∏
j
‍
(
w
j
(
n
-
1
)
(
t
)
)
δ
j
β
-
1
)
]
n
-
1
,
where [
⋯
]
m denote the Maclaurin polynomial of order m in t.Lemma 22.
Fort
∈
[
0
,
h
], h
>
0 being sufficiently small, and for each order of approximation n
≥
0, w
k
(
n
)
(
t
) in (44) is a polynomial generating function for one individual of the population; that is, it has the following properties:(a)
w
k
(
n
)
(
t
) is a polynomial of degree n in z
1
,
…
,
z
E, where the coefficient of z
1
n
1
,
…
,
z
E
n
E is O
(
t
p
), with p
=
n
1
+
⋯
+
n
E;
(b)
w
k
(
n
)
(
1
,
…
,
1
,
t
)
=
1;
(c)
Δ
k
(
n
)
(
h
)
=
w
k
(
n
)
(
h
)
-
w
k
(
n
-
1
)
(
h
)
=
O
(
h
n
);
(d)
the coefficients ofw
k
(
n
)
(
t
) regarded as a polynomial in {
z
α
} are nonnegative functions of time.Proof.
(a) The statement holds forn
=
0 and n
=
1. Assuming that it holds up to order n
-
1, we have that [
w
k
(
n
-
1
)
(
t
)
∏
j
‍
(
w
j
(
n
-
1
)
(
t
)
)
δ
j
β
]
n
-
1 is also a polynomial of degree n
-
1 in z
1
,
…
,
z
E with time-dependent coefficients of type O
(
t
p
) since each power p of t carries along a power p in {
z
α
} (as well as lower z-powers). A Picard iteration step gives
(46)
w
k
(
n
)
(
t
)
=
1
+
∑
β
‍
r
β
k
z
β
∫
0
h
‍
[
w
k
(
n
-
1
)
(
t
)
∏
j
‍
(
w
j
(
n
-
1
)
(
t
)
)
δ
j
β
]
n
-
1
d
t
-
(
∑
β
‍
r
β
k
)
∫
0
h
‍
w
k
(
n
-
1
)
(
t
)
d
t
.
The first sum contributes with one time integration and one z-power, which exactly preserves the polynomial structure and the order of the coefficients. The result is a polynomial of degree n in z
1
,
…
,
z
E, where the coefficients are polynomials of highest degree at most n in h, each one of lowest order O
(
h
p
) (as in the statement). The second integral does not alter these properties since it generates again a polynomial of degree n
-
1 in z
1
,
…
,
z
E, with time-dependent coefficients of order O
(
h
p
+
1
) or higher, up to h
n (they are just the coefficients of w
k
(
n
-
1
)
(
t
) integrated once in time).
(b) Lettingz
1
=
z
2
=
⋯
=
z
E
=
1, Picard iteration gives trivially w
k
(
n
)
(
t
)
=
1 to all orders, again by induction.
(c) The statement holds forn
=
1. Assuming it holds up to n
-
1, we compute
(47)
w
k
(
n
)
(
h
)
-
w
k
(
n
-
1
)
(
h
)
=
∑
β
‍
r
β
k
z
β
×
∫
0
h
‍
(
[
w
k
(
n
-
2
)
(
t
)
∏
j
‍
(
w
j
(
n
-
2
)
(
t
)
)
δ
j
β
]
n
-
2
[
(
w
j
(
n
-
2
)
(
t
)
+
Δ
j
(
n
-
1
)
(
t
)
)
δ
j
β
(
w
k
(
n
-
2
)
(
t
)
+
Δ
k
(
n
-
1
)
(
t
)
)
×
∏
j
‍
(
w
j
(
n
-
2
)
(
t
)
+
Δ
j
(
n
-
1
)
(
t
)
)
δ
j
β
]
n
-
1
-
[
w
k
(
n
-
2
)
(
t
)
∏
j
‍
(
w
j
(
n
-
2
)
(
t
)
)
δ
j
β
]
n
-
2
)
d
t
-
(
∑
β
‍
r
β
k
)
∫
0
h
‍
(
w
k
(
n
-
1
)
(
t
)
-
w
k
(
n
-
2
)
(
t
)
)
d
t
.
Both expressions in the first integral have the same Maclaurin polynomial up to order (
n
-
2
). Hence, the difference in the first integral contains only terms of order O
(
t
n
-
1
). This is also the case for the second integral. After integration, we obtain Δ
k
(
n
)
(
h
)
=
O
(
h
n
).
(d) The statement holds forn
=
0 and n
=
1 and in general for the zero-order coefficient, because it is unity for t
=
0, as a consequence of the initial condition. For the higher-order coefficients, we use induction again. If all z-coefficients in w
k
(
n
-
1
)
(
t
) are nonnegative, then the same holds automatically for w
k
(
n
)
(
h
) up to order n
-
1, since the O
(
h
n
) contribution adding to each inherited coefficient from the previous step does not alter its sign, for h sufficiently small. It remains to show that the coefficients corresponding to p
=
n are nonnegative. These contributions arise from the (
n
-
1
)th order of
(48)
r
β
k
z
β
[
w
k
(
n
-
1
)
(
t
)
∏
j
‍
(
w
j
(
n
-
1
)
(
t
)
)
δ
j
β
]
n
-
1
after time integration. If all δ
j
β are positive, then the expression above is a product of z-polynomials with nonnegative coefficients and the statement follows. Suppose that some δ
j
β
=
-
1. Then by Lemma 7, r
β
k
=
0 for k
≠
j. Hence,
(49)
r
β
k
z
β
w
k
(
n
-
1
)
(
t
)
∏
j
‍
(
w
j
(
n
-
1
)
(
t
)
)
δ
j
β
=
{
0
,
j
≠
k
,
r
β
k
z
β
∏
j
≠
k
‍
(
w
j
(
n
-
1
)
(
t
)
)
δ
j
β
,
j
=
k
.
However, δ
j
β
≠
-
1 for j
≠
k; that is, they are nonnegative, and hence the expression has only nonnegative coefficients.


### 5.2.1. Basic Properties of the First- and Second-Order Approximations

Let us analyse the generating functionΨ
(
z
1
,
…
,
z
E
,
t
)
=
∏
j
‍
w
j
X
j
(
0
). The simplest nontrivial applications of the approximation scheme are the first- and second-order approximations. Expanding explicitly up to second-order, we get
(50)
w
k
(
h
)
=
1
+
h
∑
β
‍
r
β
k
(
z
β
-
1
)
+
h
2
2
(
∑
β
‍
r
β
k
(
z
β
-
1
)
)
×
(
∑
γ
‍
r
γ
k
(
z
γ
-
1
)
)
+
h
2
2
∑
β
‍
r
β
k
z
β
∑
γ
‍
F
γ
β
(
z
γ
-
1
)
.We begin by noticing that the first-order approximation is obtained by neglecting all terms inh
2 in the above equation. Whenever w
j is linear in z, which is always the case for the first-order approximation, the resulting generating function Ψ corresponds to the product of (independent) multinomial distributions for the events affecting each subpopulation.Let us now considerΨ
(
z
1
,
…
,
z
E
,
t
) computed with the second-order approximated w
j’s to perform one simulation step of size h.Counting powers ofz
γ on each w
k we have (a) the probability of no events occurring in time-step h
(51)
P
k
(
h
,
0
)
=
1
-
h
(
∑
β
‍
r
β
k
)
+
h
2
2
(
∑
β
‍
r
β
k
)
2
,
(b) the probability of only one (β) event occurring in time-step h
(52)
P
k
(
h
,
I
,
β
)
=
r
β
k
{
h
-
h
2
∑
α
‍
r
α
k
-
h
2
2
∑
α
‍
F
α
β
}
,
and (c) the probability of another event (α) occurring after the first one:
(53)
P
k
(
h
,
I
I
,
β
,
α
)
=
h
2
2
r
β
k
(
r
α
k
+
F
α
β
)
.


### 5.2.2. Simulation Algorithm

A simulation step can be operated in the following way. Let(54)
∑
α
(
P
k
(
h
,
I
I
,
β
,
α
)
)
+
P
k
(
h
,
I
,
β
)
=
r
β
k
{
h
-
h
2
2
∑
α
‍
r
α
k
}
=
P
¯
k
(
h
,
I
,
β
)
be the probability of one (β) event occurring alone or followed by a second event. It follows that
(55)
P
k
(
h
,
I
I
,
β
,
α
)
=
P
¯
k
(
h
,
I
,
β
)
h
2
(
r
k
α
+
F
α
β
)
(
1
-
(
h
/
2
)
∑
α
r
α
k
)
+
o
(
h
2
)
.
The second factor above is the conditional probability P
k
(
h
,
α
/
β
) for the second event being α given that there was a first (β) event.The algorithm proceeds in two steps. First a multinomial deviate is computed with probabilities{
P
¯
k
(
h
,
I
,
β
)
} (the probability of no event occurring is still P
k
(
h
,
0
)). This gives the number n
β of occurrences of each event. In the second step, for each n
β we compute the deviates for the occurrence of a second event or no more events, with the multinomial given by the array of probabilities {
P
k
(
h
,
α
/
β
)
} defined above (the probability of no second event occurring is here 1
-
(
h
/
2
)
(
∑
α
‍
(
r
k
α
+
F
α
β
)
/
(
1
-
(
h
/
2
)
∑
α
‍
r
α
k
)
)
).


## 5.2.1. Basic Properties of the First- and Second-Order Approximations

Let us analyse the generating functionΨ
(
z
1
,
…
,
z
E
,
t
)
=
∏
j
‍
w
j
X
j
(
0
). The simplest nontrivial applications of the approximation scheme are the first- and second-order approximations. Expanding explicitly up to second-order, we get
(50)
w
k
(
h
)
=
1
+
h
∑
β
‍
r
β
k
(
z
β
-
1
)
+
h
2
2
(
∑
β
‍
r
β
k
(
z
β
-
1
)
)
×
(
∑
γ
‍
r
γ
k
(
z
γ
-
1
)
)
+
h
2
2
∑
β
‍
r
β
k
z
β
∑
γ
‍
F
γ
β
(
z
γ
-
1
)
.We begin by noticing that the first-order approximation is obtained by neglecting all terms inh
2 in the above equation. Whenever w
j is linear in z, which is always the case for the first-order approximation, the resulting generating function Ψ corresponds to the product of (independent) multinomial distributions for the events affecting each subpopulation.Let us now considerΨ
(
z
1
,
…
,
z
E
,
t
) computed with the second-order approximated w
j’s to perform one simulation step of size h.Counting powers ofz
γ on each w
k we have (a) the probability of no events occurring in time-step h
(51)
P
k
(
h
,
0
)
=
1
-
h
(
∑
β
‍
r
β
k
)
+
h
2
2
(
∑
β
‍
r
β
k
)
2
,
(b) the probability of only one (β) event occurring in time-step h
(52)
P
k
(
h
,
I
,
β
)
=
r
β
k
{
h
-
h
2
∑
α
‍
r
α
k
-
h
2
2
∑
α
‍
F
α
β
}
,
and (c) the probability of another event (α) occurring after the first one:
(53)
P
k
(
h
,
I
I
,
β
,
α
)
=
h
2
2
r
β
k
(
r
α
k
+
F
α
β
)
.


## 5.2.2. Simulation Algorithm

A simulation step can be operated in the following way. Let(54)
∑
α
(
P
k
(
h
,
I
I
,
β
,
α
)
)
+
P
k
(
h
,
I
,
β
)
=
r
β
k
{
h
-
h
2
2
∑
α
‍
r
α
k
}
=
P
¯
k
(
h
,
I
,
β
)
be the probability of one (β) event occurring alone or followed by a second event. It follows that
(55)
P
k
(
h
,
I
I
,
β
,
α
)
=
P
¯
k
(
h
,
I
,
β
)
h
2
(
r
k
α
+
F
α
β
)
(
1
-
(
h
/
2
)
∑
α
r
α
k
)
+
o
(
h
2
)
.
The second factor above is the conditional probability P
k
(
h
,
α
/
β
) for the second event being α given that there was a first (β) event.The algorithm proceeds in two steps. First a multinomial deviate is computed with probabilities{
P
¯
k
(
h
,
I
,
β
)
} (the probability of no event occurring is still P
k
(
h
,
0
)). This gives the number n
β of occurrences of each event. In the second step, for each n
β we compute the deviates for the occurrence of a second event or no more events, with the multinomial given by the array of probabilities {
P
k
(
h
,
α
/
β
)
} defined above (the probability of no second event occurring is here 1
-
(
h
/
2
)
(
∑
α
‍
(
r
k
α
+
F
α
β
)
/
(
1
-
(
h
/
2
)
∑
α
‍
r
α
k
)
)
).


## 6. Implementation: Developmental Cascades without Mortality

In this and the following Sections we apply the tools developed up to now to the description of a subprocess in insect development, namely, the development of immature larvae, which is, in principle, a highly individual subprocess. We assume that the evolution of an insect from egg to pupa occurs via a sequence ofE
>
0 maturation stages. Biologists recognise some substages by inspection in the development of larvae (instars) and even more stages when observed by other methods [35]. The immature individual may die along the process or otherwise complete it and exit as pupa or adult. The actual value of apparent maturation stages E will depend on environmental and experimental conditions and will ultimately be specified when analying experimental data in Section 7.1.Empirical evidence based on existing and new experimental results as well as biological insight produced by the present description will be the subject of an independent work [42]. We discuss however a couple of examples from the literature in Section 7.1.The events for this process are maturation events promoting individuals from one subpopulation into to the next and death events for all subpopulations. In other words, each subpopulation represents a given (intermediate) level of development and maturation events correspond to progress in the development. Development at levelj promotes one individual to level j
+
1 and the event rate is linear in the j-subpopulation. Hence, in R
α
=
∑
j
‍
r
α
k
X
j, the matrix r is square and diagonal, for 0
≤
α
≤
E
-
1 (we have exactly E
+
1 subpopulations, counted from 0 to E). All of the environmental influence is, at this point, incorporated through r.Moreover, levelE is the matured pupa and we may set r
E
=
0 since at pupation stage the present mechanism ends and new processes take place. Subpopulation E is the exit point of the dynamics and we assume it to be determined by the stochastic evolution of the immature stages 0
≤
α
≤
E
-
1.In the same spirit, the action of each event on the populations modifies just the actual population and the next one in the chain. Hence,(56)
δ
j
α
=
{
-
1
,
α
=
j
,
+
1
,
α
=
j
-
1
.
Following Section 3, [
R
α
0
]
T
=
(
r
0
M
,
0
,
…
,
0
) and
(57)
F
α
β
=
{
-
r
α
,
β
=
α
,
+
r
α
,
β
=
α
-
1
,
0
,
otherwise
.
Finally, letting m
T
=
-
M
(
1,1
,
…
,
1
) we can write R
α
0
=
∑
β
‍
F
α
β
m
β. In this framework the index pair (
j
,
α
) is highly interconnected and we can achieve a complete description using just one index. We use 0
≤
j
≤
E
-
1 throughout this section.The procedure of Section3 gives Ψ
(
{
z
}
,
t
)
=
∏
j
‍
(
f
j
)
m
j
=
∏
j
‍
(
f
j
-
M
), while the dynamical system for f
j reads, recalling the recasting in Remark 21;
(58)
f
˙
j
=
f
j
(
-
r
j
(
z
j
f
j
-
1
)
+
r
j
+
1
(
z
j
+
1
f
j
+
1
-
1
)
)
,
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
j
=
0
,
…
,
E
-
2
,
and f
˙
E
-
1
=
f
E
-
1
(
-
r
E
-
1
(
z
E
-
1
f
E
-
1
-
1
)
). The overall initial condition is f
j
(
0
)
=
1. Let θ
j
=
log
⁡
(
f
j
) and
(59)
ψ
j
=
∑
β
=
j
E
-
1
‍
θ
β
;
W
j
=
exp
⁡
(
-
ψ
j
)
.
Note that W
j
f
j
=
W
j
·
exp
⁡
(
θ
j
)
=
W
j
+
1. Rewriting the equations,
(60)
θ
˙
j
=
-
r
j
(
z
j
exp
⁡
(
θ
j
)
-
1
)
+
r
j
+
1
(
z
j
+
1
exp
⁡
(
θ
j
+
1
)
-
1
)
,
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
j
<
E
-
1
,
ψ
˙
j
=
-
r
j
(
z
j
exp
⁡
(
θ
j
)
-
1
)
=
θ
˙
j
+
ψ
˙
j
+
1
;
j
<
E
-
1
,
θ
˙
E
-
1
=
ψ
˙
E
-
1
=
-
r
E
-
1
(
z
E
-
1
exp
⁡
(
θ
E
-
1
)
-
1
)
,
with initial condition θ
j
(
0
)
=
0 for all j. With this notation the generating function reads Ψ
(
{
z
}
,
t
)
=
W
0
M, and we formulate hence an ODE for the W
j’s:
(61)
W
˙
j
+
r
j
W
j
=
r
j
z
j
W
j
+
1
,
W
j
(
0
)
=
1
;
j
=
0
,
…
,
E
-
2
,
W
˙
E
-
1
=
r
E
-
1
(
z
E
-
1
-
W
E
-
1
)
,
W
E
-
1
(
0
)
=
1
.
By defining the auxiliary quantity W
E
(
t
)
=
1, the solution for the W
j-differential equations can be written in recursive form:
(62)
W
j
(
t
)
=
e
-
r
j
t
+
z
j
r
j
∫
0
t
‍
e
r
j
(
s
-
t
)
W
j
+
1
(
s
)
d
s
,
j
=
0
,
…
,
E
-
1
.
In Appendix B we give the details of the explicit solution of this problem in terms of Laplace transforms. Finally, we obtain
(63)
W
0
(
t
)
=
(
∏
m
=
0
E
-
1
‍
z
m
)
(
1
-
∑
k
=
0
E
-
1
‍
exp
⁡
(
-
r
k
t
)
∏
j
=
0
,
j
≠
k
E
-
1
‍
r
j
r
j
-
r
k
)
+
∑
k
=
1
E
‍
∏
j
=
0
E
-
k
z
j
z
E
-
k
r
E
-
k
∑
m
=
0
E
-
k
‍
r
m
exp
⁡
(
-
r
m
t
)
∏
l
=
0
,
l
≠
m
E
-
k
‍
r
l
r
l
-
r
m
.The solution looks more compact when allr
k are assumed to take the common value r:
(64)
W
0
(
t
)
=
∏
m
=
0
E
-
1
z
m
(
E
-
1
)
!
∫
0
r
t
‍
exp
⁡
(
-
s
)
s
E
-
1
d
s
+
∑
k
=
1
E
‍
∏
m
=
0
E
-
k
z
m
z
E
-
k
(
r
t
)
E
-
k
exp
⁡
(
-
r
t
)
(
E
-
k
)
!
.


## 7. Developmental Cascades with Mortality

We add to the previous descriptionE subpopulation specific death events, Q
α
=
q
α
X
α, where α
=
0
,
…
,
E
-
1. We do not consider death of the pupa as an event for the same reasons as before: the pupa subpopulation X
E leaves the present framework of study. We have actually event pairs for each subpopulation and hence the analysis can still be done with just one index α
=
0
,
…
,
E
-
1, since each event is population specific. The calculations get however more involved. It is practical to organise the events as follows: Q
0
,
R
0
,
Q
1
,
R
1
,
…
,
Q
E
-
1, and R
E
-
1. Then, (
(
Q
,
R
)
0
)
T
=
(
q
0
M
,
r
0
M
,
0
,
…
,
0
). We have
(65)
δ
α
2
α
=
-
1
,
(
death
)
,
α
<
E
,
δ
α
2
α
+
1
=
-
1
δ
α
+
1
2
α
+
1
=
+
1
}
,
(
development
)
,
α
<
E
.
F
2
α
β
=
{
-
q
α
,
β
=
2
α
,
-
q
α
,
β
=
2
α
+
1
,
q
α
,
β
=
2
α
-
1
,
F
2
α
+
1
β
=
{
-
r
α
,
β
=
2
α
+
1
,
-
r
α
,
β
=
2
α
,
r
α
,
β
=
2
α
-
1
,
α
<
E
.
Recall that the indices in F run from 0 to 2
E
-
1. It is also practical to order the variables f
α and the z
α in pairs (
g
,
f
) and (
z
,
y
) as (
g
0
,
f
0
,
…
,
g
E
-
1
,
f
E
-
1
) and (
z
0
,
y
0
,
…
,
z
E
-
1
,
y
E
-
1
). Also for notational convenience we will assume that there exist q
E
=
r
E
=
0 when necessary. Defining
(66)
H
α
(
g
α
,
f
α
)
=
q
α
(
z
α
g
α
-
1
)
+
r
α
(
y
α
f
α
-
1
)
,
we have (recall again Remark 21 and note that H
E
≡
0),
(67)
g
˙
α
=
-
H
α
(
g
α
,
f
α
)
g
α
,
f
˙
α
=
(
-
H
α
(
g
α
,
f
α
)
+
H
α
+
1
(
g
α
+
1
,
f
α
+
1
)
)
f
α
,
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
α
=
0
,
…
,
E
-
2
,
g
˙
E
-
1
g
E
-
1
=
f
˙
E
-
1
f
E
-
1
=
-
H
E
-
1
(
g
E
-
1
,
f
E
-
1
)
.
The overall initial condition is still g
α
(
0
)
=
1
=
f
α
(
0
).Since there are more events than subpopulations, according to Definition1 and Corollary 14, F has structural zeroes corresponding to constants of motion in the associated dynamical system. The zero eigenvectors of F are V
α
T
=
(
0
,
…
,
-
1,1
,
1,0
,
…
), nonzero in positions 2
α, 2
α
+
1, and 2
α
+
2, for α
=
0
,
…
,
E
-
1 (for the last one, recall that there is no position “2
E"). Definitions similar to those used in the previous section will be useful as well; namely, θ
α
=
log
⁡
g
α and ϕ
α
=
log
⁡
f
α, changing H
α and the initial conditions correspondingly. Hence,
(68)
ϕ
˙
α
=
θ
˙
α
-
θ
˙
α
+
1
⟹
ϕ
α
=
θ
α
-
θ
α
+
1
,
α
=
0
,
…
,
E
-
2
,
θ
˙
α
=
-
H
α
(
θ
α
,
θ
α
-
θ
α
+
1
)
,
α
=
0
,
…
,
E
-
2
,
ϕ
˙
E
-
1
=
θ
˙
E
-
1
=
-
H
E
-
1
(
θ
E
-
1
,
θ
E
-
1
)
.
The structural zeroes and constants of motion allow solving all equations related to the ϕ-variables and only a chain of θ-equations is left, very much in the spirit of the previous section. m
T
=
-
M
(
1,0
,
…
,
0,0
) satisfies R
0
=
F
m. Finally, Ψ
(
{
z
,
y
}
,
t
)
=
g
0
-
M. We now compute g
0
(
t
).LetW
α
=
e
-
θ
α
=
g
α
-
1, 0
≤
α
≤
E
-
1. Defining the auxiliary variable W
E
(
t
)
=
1, the dynamical equations for θ in terms of W read
(69)
W
˙
α
+
(
q
α
+
r
α
)
W
α
=
q
α
z
α
+
r
α
y
α
W
α
+
1
,
W
α
(
0
)
=
1
,
with solution
(70)
W
α
(
t
)
=
e
-
(
q
α
+
r
α
)
t
+
∫
0
t
‍
e
-
(
q
α
+
r
α
)
(
t
-
s
)
(
q
α
z
α
+
r
α
y
α
W
α
+
1
(
s
)
)
d
s
.
In Appendix B we give the details of the explicit solution of this problem in terms of Laplace transforms. Finally, we obtain
(71)
W
0
(
t
)
=
∏
m
=
0
E
-
1
‍
r
m
y
m
q
m
+
r
m
-
∑
m
=
0
E
-
1
‍
r
m
y
m
e
-
(
q
m
+
r
m
)
t
q
m
+
r
m
×
∏
k
=
0
,
k
≠
m
E
-
1
‍
r
k
y
k
q
k
-
q
m
+
r
k
-
r
m
+
∑
k
=
1
E
‍
1
[
r
y
]
E
-
k
×
[
∏
j
=
0
,
j
≠
m
E
-
k
‍
∏
m
=
0
E
-
k
‍
r
m
y
m
q
m
+
r
m
[
q
z
]
E
-
k
+
∑
m
=
0
E
-
k
‍
(
1
-
[
q
z
]
E
-
k
q
m
+
r
m
)
x
x
x
x
x
x
x
x
x
x
x
x
×
r
m
y
m
e
-
(
q
m
+
r
m
)
t
×
∏
j
=
0
,
j
≠
m
E
-
k
‍
r
j
y
j
q
j
-
q
m
+
r
j
-
r
m
]
.
Again we have that Ψ
(
{
z
,
y
}
,
t
)
=
W
0
M. Specialising for the case where all r and all q coefficients are equal, the solution reads
(72)
W
0
(
t
)
=
∏
m
=
0
E
-
1
r
y
m
(
E
-
1
)
!
∫
0
t
‍
e
-
(
q
+
r
)
x
x
E
-
1
d
x
+
∑
k
=
0
E
-
1
‍
∏
m
=
0
k
r
y
m
r
y
k
k
!
[
t
k
e
-
(
q
+
r
)
t
+
q
z
k
∫
0
t
‍
e
-
(
q
+
r
)
x
x
k
d
x
]
.


### 7.1. Examples from the Literature

In [19], Section 3, the drug delivery process by oral ingestion of a medicine consisting of microspheres containing the active principle is discussed. The process is modeled considering that microspheres enter the stomach (E
=
0) and either dissolve (passing subsequently to the circulatory system) or proceed to the next stage of the GI tract, namely, the duodenum. Both processes are assumed to obey linear transition rates. In the same way, the undissolved duodenal microspheres pass to the next intestinal stage—the jejunum—and later to the ilium. The remaining undissolved microspheres entering the colon are eliminated without further dissolution. This is a neat example of a cascade, where dissolution corresponds to mortality in (71), while maturation corresponds to progress through the different stages. Remaining particles exit the system upon arrival to the fifth and last stage (E
=
4). Using the data listed in page 239 of [19] for the dissolution and passage coefficients, the results of that paper can be directly read from (71). The probabilities of staying or dissolving at a given stage, function of time, appear as coefficients of W
0 for different powers of z
k and y
k. As an example, we compute the probabilities p
i
(
t
) of permanence at a given stage (as a function of time), using the λ- and μ-notation from page  239 in [19] (corresponding in (71) to r and q, resp.). Define first P
(
s
)
=
∏
m
=
0
s
‍
λ
m, s
=
0,1
,
2; Q
(
s
,
m
)
=
∏
k
=
0
,
k
≠
m
s
‍
(
λ
k
-
λ
m
+
μ
k
-
μ
m
), s
=
1,2
,
3, for each integer m
∈
[
0
,
s
] and finally R
(
s
)
=
λ
s
+
μ
s, s
=
0
,
…
,
3. With these ingredients the probabilities are p
0
(
t
)
=
exp
⁡
(
-
t
R
0
), while for i
=
1,2
,
3 we have p
i
(
t
)
=
P
(
i
-
1
)
∑
m
=
0
i
‍
(
exp
⁡
(
-
t
R
m
)
/
Q
(
i
,
m
)
).


## 7.1. Examples from the Literature

In [19], Section 3, the drug delivery process by oral ingestion of a medicine consisting of microspheres containing the active principle is discussed. The process is modeled considering that microspheres enter the stomach (E
=
0) and either dissolve (passing subsequently to the circulatory system) or proceed to the next stage of the GI tract, namely, the duodenum. Both processes are assumed to obey linear transition rates. In the same way, the undissolved duodenal microspheres pass to the next intestinal stage—the jejunum—and later to the ilium. The remaining undissolved microspheres entering the colon are eliminated without further dissolution. This is a neat example of a cascade, where dissolution corresponds to mortality in (71), while maturation corresponds to progress through the different stages. Remaining particles exit the system upon arrival to the fifth and last stage (E
=
4). Using the data listed in page 239 of [19] for the dissolution and passage coefficients, the results of that paper can be directly read from (71). The probabilities of staying or dissolving at a given stage, function of time, appear as coefficients of W
0 for different powers of z
k and y
k. As an example, we compute the probabilities p
i
(
t
) of permanence at a given stage (as a function of time), using the λ- and μ-notation from page  239 in [19] (corresponding in (71) to r and q, resp.). Define first P
(
s
)
=
∏
m
=
0
s
‍
λ
m, s
=
0,1
,
2; Q
(
s
,
m
)
=
∏
k
=
0
,
k
≠
m
s
‍
(
λ
k
-
λ
m
+
μ
k
-
μ
m
), s
=
1,2
,
3, for each integer m
∈
[
0
,
s
] and finally R
(
s
)
=
λ
s
+
μ
s, s
=
0
,
…
,
3. With these ingredients the probabilities are p
0
(
t
)
=
exp
⁡
(
-
t
R
0
), while for i
=
1,2
,
3 we have p
i
(
t
)
=
P
(
i
-
1
)
∑
m
=
0
i
‍
(
exp
⁡
(
-
t
R
m
)
/
Q
(
i
,
m
)
).


## 8. Concluding Remarks

The description of stochastic population dynamics at event level has several advantages with respect to the description at population level. While the projection onto population space is straightforward, the description of the number of events carries biological information that can be lost in the projection. From a practical, mathematical, point of view, the description in terms of events is more regular and makes room for a general discussion, since the number of events always increases by one.When individual populations' processes are considered, the event rates become linear and, in addition, if the result of the event is to decrease one (sub)population, then they can only be proportional to the decreased population.The probability generating function for the number of events,{
n
β
}, conditioned to an initial state described by population values {
X
i
(
0
)
}, has been obtained as a function of the solutions of a set of ODE known as the equations of the characteristics but presented in a form oriented towards numerical calculus rather than the usual presentation oriented towards geometrical methods.We have introduced short-time approximations that, to the lowest order, are in coincidence with previously considered heuristic proposals. Our presentation is systematic and can be carried out up to any desired order. In particular, we have shown how the second-order approximation can be implemented.The general solution represents independent individuals identically distributed within their (sub)population compartment.Additionally, we believe that obtaining general solutions of the linear problem is a sound starting point when the more complicated approximations for nonlinear rates are considered.Further, we have solved a stochastic individual developmental model, that is, a model that at the population level is described by rates that are linear in the (sub)populations. The model not only describes just the populations but also describes the events that change the populations.



---

*Source: 873624-2014-02-17.xml*