Multinomial approximation to the Kolmogorov Forward Equation for jump (population) processes

Abstract: We develop a simulation method for Markov Jump processes with finite time steps based in a quasilinear approximation of the process and in multinomial random deviates. The second-order approximation to the generating function, Error = O(dt), is developed in detail and an algorithm is presented. The algorithm is implemented for a Susceptible-Infected-Recovered-Susceptible (SIRS) epidemic model and compared to both the deterministic approximation and the exact simulation. Special attention is given to the problem of extinction of the infected population which is the most critical condition for the approximation.


Introduction
In a broad sense, stochastic population processes correspond to the time-evolution of countable sets (and subsets) of individuals in interaction. The description is performed by tracking the number of members of each relevant subpopulation over time. The classification scheme used for the populations is dictated by the problem. For example, in an epidemic problem, it is often sensible to group the population in at least two sets: susceptible, S and infectious, I. In the description of the lifecycle of an insect, we will often refer to subsets (compartments) corresponding to different developmental stages such as eggs, larvae, pupae and adults. Such models can be used in a large range of problems in chemistry and biology. The quality of such models largely ABOUT THE AUTHOR The authors participate in several long-term collaborations concerning the propagation of epidemic diseases in different populations (ranging from human groups to crops) with emphasis in those propagated by vectors, thus involving ecological aspects. The aim of this cross-disciplinary group is to develop the understanding of general population dynamics problems combining mathematical and biological expertise in a way that goes beyond the individual disciplines.

PUBLIC INTEREST STATEMENT
The overall goal of population dynamics is to understand the collective time-evolution of populations in terms of their reciprocal and environmental influences. The description is made in terms of the number of individuals that share the same relevant characteristics and results in a intrinsic (unavoidable) stochasticity. This manuscript describes a new implementation of the core stochastic methods. The new method takes advantage of particularities that more often than not appear in biological populations such as in the description of biological cycles and/or the progress in time of epidemic outbreaks. The present method is more efficient than a general approach, while keeping accuracy within desirable levels. Understanding, modelling and predicting population processes is an important ingredient in the development of public policies such as health or farming, to mention just a few.
depends on the quality of the transcription between natural science and mathematics. For example, stochastic compartmental models applied to vector-borne diseases such as yellow fever perform very satisfactorily when compared to real (historic) data (Fernández, Otero, Schweigmann, & Solari, 2013). The vector component of the model can be used to forecast the expected range of the vector (in this case, the mosquito Aedes aegypti) (Otero, Solari, & Schweigmann, 2006), yielding results that agree well with field studies (Zanotti et al., 2015). However, if the modelling strategy is going to be useful, it is a practical requirement that it shall use a reasonable amount of computational resources. The general family of population processes described in this work is associated with the probabilities ruled by the Kolmogorov Forward Equation (KFE) for jump processes (Kolmogoroff, 1931), since populations update by jumps (events) such as the hatching of an egg that decreases the egg-population by one and increases the larvae population by one.
The description of population problems by means of stochastic methods has a modern history of over 100 years (McKendrick, 1914) (though its origins can be traced further back to Malthus, 1798). This approach has several advantages over other existent approaches. In the first place, it recognises the individual level and the intrinsic impossibility of describing this level in a deterministic way. Further, it deals with integer variables, which is the only consistent way of dealing with individuals at all levels of presence (from extinction to groups of billions). Further, for systems where the actual evolution time is relevant, the use of Markov Jump processes (Ethier & Kurtz, 1986) (continuous-time Markov chains) provides a rationale for organising the dynamics. The relation between the stochastic description of population problems and thepossibly more populardifferential equations approach has been studied and established by Kurtz (see Ethier & Kurtz, 1986), where the validity and limitation of the so-called deterministic limit is discussed (see also Aparicio, Natiello, & Solari, 2012).
The "jumps" in the Markov Jump process are the symbolic description of the dynamical evolution, identifying situations where the population suffers drastic changes at a point in time (e.g., an individual is born). This change is called an event. The list of events will control the outcome possibilities of a model. Events come in different flavours, characteristic times, and so on and the more accurate the description in biological terms, the larger the number of events that need to be considered.
The origin of this approach goes back to Kolmogorov (Kolmogoroff, 1931) and, more specifically, to Feller's work on stochastic processes (Feller, 1949) developing from the equations that he named Kolmogorov equations. Subsequently, Kendall (1949Kendall ( , 1950 devised an algorithm to implement Feller's approach. Event-based models, though powerful and successful (see above), are computationally timeconsuming. The approach requires to track all individual events in their given time order. On the contrary, experimental data are usually collected and grouped within reasonable time intervals (in an epidemic outbreak, we speak of number of cases per week, while e.g., municipal housing policies consider only the balance between birth, death, emigration and immigration over a year). Eventually, we arrive to a situation where approximations to Kendall's algorithm are required. Rather than computing a large number of individual events, it would be desirable to estimate the overall event-count within a fixed time step, drastically shortening computation resources without loosing accuracy.
Approximating a stochastic process with another stochastic process is an interesting issue by itself. Consider families of processes such as e.g., KðNÞ, a Kendall problem for a population of size N and AðNÞ its approximation. If the approximation should have a chance to be satisfactory, we should have, in some abstract sense, that KðNÞ ¼ AðNÞ þ EðNÞ, where EðNÞ is the error of the approximation. A minimum, necessary, requirement is of course that lim N!1 EðNÞ ¼ 0, but this is not enough. In all situations, simulations and modelling involve a finite value (or range) of N. A way to estimate and control EðNÞ for any N is needed. Things do not get easier when one realises that since EðN 0 ÞÞ0 for any finite N 0 , the statistics of KðN 0 Þ and AðN 0 Þ will differ, eventually, when the number of repetitions is very large. Approximations involve hence different tolerances: (a) small EðNÞ so that it makes sense to attempt to approximate K with A, and (b) medium large statistics so that while A satisfactorily uncovers the dynamical effects of the system in study, its stochastic difference with K can still be considered negligible.
In previous works, we have developed a Poisson approximation (Solari & Natiello, 2003) to Kendall's algorithm as well as a general view of approximation methods (Solari & Natiello, 2014). These approximation methods are organised building upon the concept of linear events (those where the probability rate depends linearly on the involved populations) along with a consistent multinomial approximation for the linear situation with constant (in time) coefficients.
The goal of this work is to develop a new consistent multinomial approximation to Kendall processes based on a linear approximation to event rates (called quasilinear approximation), allowing for general time-dependencies in the event rates. The manuscript is organised closely following Solari and Natiello (2014), of which it is an extension and in some form a natural development. After a section refreshing the idea of KFE, we state the quasilinear approximation in Section 3. This section starts by presenting and defining the approximation, subsequently formulating the computation of averages within it. Then, the error estimates of the approximation, which are necessary to control the accuracy of the implementation are computed in Section 3.2 and finally the generating function is computed. The section ends with the statement of the stochastic dynamical equations, in Lemmas 3 and 4, for the second-order approximation. These lemmas are the central tools in the following sections dealing with examples. Section 4 deals with a well-known example, rendered technically difficult by the time-dependency. Section 5 dwells in a more elaborated example, which is the source for numerical tests. Some concluding remarks are found in Section 6.

Kolmogorov Forward Equation
Let ðn 1 ; n 2 ; Á Á Á ; n E Þ ¼ n denote a state of the system with the count of how many events of each type has occurred up to a given time. Denote the product z n 1 1 z n 2 2 Á Á Á z n E E by z n and the generating function by Φ ¼ ∑ n z n Pðn; t; X 0 Þ, where Pðn; t; X 0 Þ is the probability of having n events at time t given the initial population state X 0 . Finally, X is the population array ðX 1 ; X 2 ; Á Á Á ; X N Þ, where the subpopulations satisfy X j ¼ X 0 j þ ∑ α n α δ α j . The transition matrix δ α j denotes the modification acted by event α on population j (Solari & Natiello, 2003) and it has integer entries.
For a general Markov Jump process, the KFE for the generating function reads (Solari & Natiello, 2003), where W α ðXðnÞÞ is the transition rate associated with event α, in other words, the probability per unit time for the occurrence of event α given the state XðnÞ.
The population array X ¼ ðX 1 ; X 2 ; Á Á Á ; X N Þ denotes the count of individuals X i in each subpopulation i (at a given time t or after the occurrence of n events) which is necessarily non-negative. The following definition is therefore proper: Definition 1. A set of states U, with the property that if X 2 U at t 0 then X 2 U for all t ! t 0 is called positively invariant.
Positively invariant sets are often called trapping sets. For example, the extinction state where all subpopulations are zero is a trapping set. The set of states with non-negative populations is a positively invariant set as well.

Quasilinear approximation
A common feature in all population problems involving higher organisms is that the events modifying the system can be regarded as belonging to only three classes: (1) Individual death events where the relevant consequence is to reduce some subpopulation in one unit, while no other subpopulation is modified.
(2) Development events where one subpopulation increases by one and another decreases by the same amount (e.g., a transition from pupa to adult in insect development).
(3) Creation events where some populations decrease by a total of k members and some other populations increase by a larger number k 0 > k. Quite often k ¼ 1. Consider for example oviposition in insects, where the number of fecundated females decreases by one unit and the number of eggs increases by k 0 ! 1.
An interesting fact in biological processes is that rarely, if ever, an event can exist that does not decrease any population. Biological "creation" processes always involve the modification of some previous entity (usually one individual of another subpopulation which in the act of creation ceases to exist). In this poetical sense, becoming adult decreases the young population in one, being born decreases the pregnant female population in one, becoming ill decreases the healthy (often called susceptible) population in one, etc., and no other subpopulation is decreased in each case. Hence, we will hereafter consider processes where the events α ¼ 1; Á Á Á ; E may be sorted according to the subpopulation JðαÞ that they decrease. Clearly, not all population problems are linear as in (Solari & Natiello, 2014). However, similar arguments to those used to prove Lemma 7 in (Solari & Natiello, 2014) apply here, this is, Lemma 1. Let the polynomial transition rate W α be such that it only decreases the subpopulation JðαÞ in one unit, i.e., δ α JðαÞ ¼ À1 and δ α i ! 0. Population space is positively invariant only if W α ¼ X JðαÞ f α ðXÞ: (2) Proof. Positive invariance demands that for any sufficiently short time h the probability P α ðhÞ of occurrence of an event α such that X JðαÞ ðt þ hÞ ¼ X JðαÞ ðtÞ À 1 < 0 is zero (we are just saying that no event can push a subpopulation into negative counts). Under the Markov Jump Process assumptions (Durrett, 2001;Feller, 1940), this probability may be written as P α ðhÞ ¼ W α ðXðtÞÞh þ oðhÞ. Hence, W α ðX 1 . . . 0 Jα . . .Þ ¼ 0 and we can extract a factor of X JðαÞ thus proving the relation.
Polynomial smoothness gives the desired result. ◽ Therefore, we introduce the quasilinear approximation, namely we let the transition rates be approximated by where JðαÞ is the subpopulation decreased by α. We shall denote by ψðz; tÞ ¼ ∑ n z n Pðn; tÞ the quasilinear approximant of Φ and P the quasilinear approximant of the probabilities. The quasilinear approximation consists in replacing f α ðXÞ by a self-consistent average < f α ðXÞ > ψ (to be specified later in detail), i.e., an average based on the generating function ψ. In other words, we extract a linear part associated with the decreasing subpopulation on each W α (in the situation where W is actually linear, f is just a constant). The quasilinear approximation allows for further computation. Indeed, cf. (Solari & Natiello, 2014;Equation 10).
An important feature of this work is that the quasilinear approach is an approximation. Hence, error bounds relative to the full-problem must be considered. We discuss this issue in subsection 3.2.
Since the approximated problem is linear, most of the results in Solari and Natiello (2014) can be transported or generalised to the present work, once the appropriate changes for possible timedependencies in the quasilinear coefficients r k β are introduced. For the present problem, the main specification reads (the δ-symbol with two population indices is the Kronecker delta) and it will introduce a timedependency through the expectation values.

Dynamical equation for the average populations
Starting from the population values where it is sufficient to sum over all events β such that δ β k Þ0. Recalling Equations (3) and (4) (JðαÞ is the population decreased by event α), we obtain expressions within the quasilinear approximation: This result has an equivalent in (Solari & Natiello, 2003;Equation 48), the large N limit: Lemma 2. If lim Ω!1 < X k Ω > ¼ x k exists and W β ðXÞ is polynomial in X, then under the generalised mass action law (Ethier & Kurtz, 1986;Solari & Natiello, 2003) being Ω a typical size for the system) the deterministic limit of the population dynamics Proof. From Equation (6), we obtain By the generalised mass action law, Since the right-hand side is finite, we can also interchange the limits in the left-hand side.

Error estimates
We want to reconsider the error estimates of (Solari & Natiello, 2003) in order to assess the accuracy of the approximation. The procedure in (Solari & Natiello, 2003; Sections 3 and 4) is general enough and it can be reproduced here, taking care of the proper expression for the W's and the generating function.
Bounds to the error of the approximation Δ ¼ Φ À ψ can be computed by integrating the difference between the exact and approximate time-evolutions from the same initial condition. We start by defining the auxiliary operator L formally solving the exact time-evolution as This is achieved by expanding the right-hand side of Equation (1) using Equation (2) and reordering the sum over events. B k denotes the set of events that decrease the subpopulation k. Note that for β 2 B k , δ β JðβÞ ¼ À1 and r j β ¼ 0 for j Þ JðβÞ (from Equation 5), being then JðβÞ ¼ k.
The approximate time-evolution in the quasilinear approximation is recovered putting together Equations (1) and (3), finally leading to We can now compute the error estimate as: given that ψðz; 0; X 0 Þ ¼ Φðz; 0; X 0 Þ (Φ being the exact generating function). Note that all operators in the square bracket depend on s:. Hence, for suitable bounded positive constants A and B, using Grönwall's inequality it also holds that ðz α ψðz; t À s; X 0 þ δ α k Þ À ψðz; t À s; X 0 Þ h i Aðt À sÞe BðtÀsÞ : Finally, for a suitable bounded positive constant, C it follows that Δðz; t; X 0 Þ C B 2 1 þ e Bt ðBt À 1Þ The present approach has smaller upper bounds than the method developed in (Solari & Natiello, 2003), since the constant C is proportional to B: Thus, only the nonlinear processes contribute to the error.
The conditions that optimise the estimation of average values and dispersions require that (Solari & Natiello, 2003; Equations 26 and 51) where ψðz; t; X 0 Þ is the approximated generating function, to be considered in the next section. The later expression can be read where the averages are taken with the approximated generating function with initial condition

Generating function
The quasilinear generating function ψðz; t 0 ; tÞ can be computed taking advantage of the linearity properties of the approximation (Solari & Natiello, 2014;Equation 35). For the present work, rendering time-dependency explicit, we have where φ α ðz; t; τÞ satisfies the equation to be integrated over 0 τ t À t 0 , and the coefficients m β express the initial condition in terms of the event matrix δ: X k ðt 0 Þ ¼ ∑ β δ β k m β . Recasting the result in population space (Solari & Natiello, 2014;Equation 36), we may write d dτ ln φ α ðz; t; τÞ ð Þ¼∑ k δ α k d dτ ln w k ðz; t; τÞ ð Þ , (where w k is the generating function for one individual of the population k, see Lemma 5 in Appendix A) thus arriving to a differential equation for the generating function for each subpopulation to be integrated over 0 τ t À t 0 (cf. Solari & Natiello, 2014;Equation 43): d dτ ln w k ðz; t; τÞ ð Þ¼∑ β r k β ðt À τÞ z β Y j w δ β j j ðz; t; τÞ À 1 0 @ 1 A ; w k ðz; t; 0Þ ¼ 1: We recall that for β 2 B k , δ β JðβÞ ¼ À1 and r j β ¼ 0 for j Þ JðβÞ. Hence, for β 2 B k we have JðβÞ ¼ k, ðδ β k þ 1Þ ¼ 0 and d dτ w k ðz; t; τÞ þ ∑ β 2 B k r k β ðt À τÞw k ðz; t; τÞ ¼ ∑ β 2 B k r k β ðt À τÞz β Q jÞk w δ β j j ðz; t; τÞ; w k ðz; t; 0Þ ¼ 1: We rewrite Equation (10) as an integral equation. We set first H k ðt; t 0 ; aÞ ¼ expðÀ Hence, The first term is the probability of no events, while the second term is built with two objects, The first contribution gives the probability of the first event taking place precisely in the interval x; x þ dx ½ while the second factor is the generating function for an individual after the occurrence of an event β at time u. The factor z β counts the event β.

Ordinary differential equation form of the approximation
The approximation lemma in Appendix A establishes the relevant properties of the generating functions w k within the approximation. Picard's iteration scheme starts with w ð0Þ k ðz; t; t À t 0 Þ ¼ 1. The next-order approximation obeys (always with the initial condition w k ðz; t; 0Þ ¼ 1 at any order) with the solution w ð1Þ k ðz; t; t À t 0 Þ ¼ 1 þ ðz β À 1Þ 1 À # ðnÞ ðz; t; τÞ: When computing specific problems, the coefficients r k β may actually depend on the solutions w k and a further approximation scheme will be required to approximate the time integrals. A more manageable expression is obtained by changing variables in Equation (12)  Then, Equation (12) reads,

Second-order approximation
From the approximation lemma in Appendix A, we know that w ðnÞ k is a n-th degree polynomial in the variables z β , with time-dependent coefficients. We propose the following approximated expression for w k : since we will only need event strings of length at most two to implement the second-order approximation.
Lemma 3. The following relations hold, for β 2 B k and α 2 B j : Proof. Straightforward substitution of Equation (14) into Equation (13) gives the result by separating terms after powers of z. ◽ Note also that P kj βα can be decomposed as the sum of δ β j identical probabilities per generated individual, Each member of subpopulation k undergoes an event β with probability P k β . In terms of the population, n β events are produced, with n β a random variable distributed with Multinomial P k β n o ; X k ð0Þ .
The occurrence of this event produces an increment δ β j n β ¼ Δ j β in the subpopulations j : δ β j > 0. The newly arrived members of the population j may in turn undergo a new transition triggered by event α (or no transition at all). The probability of this sequence βα is Q βα ¼ P kj βα δ β j P k β , i.e., the probability of event α following event β for each of the new individuals δ β j , conditioned to the actual occurrence of β. Thus, there will be n α additional events subtracting individuals from the j subpopulation as a consequence of event α 2 B j . This will modify in turn the subpopulations in the set l δ l α > 0 È É by an amount Δ l βα ¼ δ l α n α . The value of n α is drawn from Multinomial Q βα Lemma 4. The approximated generating function Equation (14) corresponds to a concatenation of multinomial processes.
Proof. Consider a process described by strings of length two Ξ ¼ βα such that β 2 B k ; α 2 B j ; δ β j >1. Then, where ψ =β ðz α Þ is given by the conditional probability and corresponds to a multinomial The generating function is The generating function of the marginal probability is itself a multinomial distribution resulting from the competing exponential events. Finally, 14 corresponds to the first order in Δt of the development of ð1 þ ∑ α 2 B j ððz α À 1ÞQ βα Þ The proof actually applies to any order of the truncation of 12 which results always in the concatenation of multinomial distributions arising from the generating function for the marginal probabilities (which are themselves multinomial since they are the consequence of competing exponential events).

Example: SIR
Equipped with the lemmas of the previous section, we are in a position to compute an approximate generating function and approximate dynamics keeping the steps of the time-evolution sufficiently low to guarantee error control.
Let us consider a SIR epidemic system. There are three subpopulations: Susceptible, Infected and Recovered. At a given time t 0 we have the initial values S 0 , I 0 and R 0 . We are interested in the dynamical outcome at a later time t and hence τ ¼ t À t 0 . The transitions are as follows: (1) Contagion δ 1 S ¼ À1; There are no events decreasing the recovered population, hence w R ¼ 1. Further, at time t we have: ¼ < S 0 À n 1 >< I 0 þ n 1 À n 2 > À < ðn 1 À < n 1 >Þðn 1 À n 2 À < n 1 À n 2 >Þ >: The expected average values at time t for a process initiated at time t 0 read, The last line can be computed in two equivalent ways. Either as stated, using that the individual expectation values correspond to adequate z-derivatives of the quasilinear generating function evaluated at z ¼ 1, or noting that for this problem n 2 can be splitted as a contribution independent of event 1, plus a contribution following the occurrence of event 1, i.e., n 2 ¼ n 2 þ n 12 . While n 2 is independent of n 1 , n 12 is distributed with Binomial P SI 12 P S 1 ; n 1 . Hence, < n 1 n 2 > À < n 1 >< n 2 > ¼ < n 1 n 12 > À < n 1 >< n 12 > ¼ < n 2 1 > À < n 1 > 2 À Á P SI 12 ðt;tÀt 0 Þ P S 1 ðt;tÀt 0 Þ since from the joint distribution Pðn 1 ; n 12 Þ ¼ Pðn 1 ÞP n 12 n 1 we have that the expectation value of n 12 conditioned to n 1 is < n 12 > n 1 ¼ n 1 P SI 12 P S 1 and < n 1 n 12 > ¼ < n 1 < n 12 > n 1 > ¼ < n 2 1 > P SI 12 P S
The stochastic number of events occurring in a short time-interval h starting at t 0 can be obtained as follows: (1) Solve 15 and obtain the coefficients P k ðhÞ.
(3) For the case of contagion events, we note that p ¼ P SI 12 δ 1 I P S 1 ¼ P SI 12 P S 1 and compute the stochastic number of events following a contagion as 1 À pð1 À zÞ ð Þ n1δ I 1 ¼ ð1 À pð1 À zÞÞ n1 .

Example: SIRS
In this example, we extend the previous results including loss of immunity of recovered individuals and provide numerical computations within the quasilinear approximation. While in the SIR model of the previous section an individual will suffer no other events after following the sequence S ! I ! R, the loss of immunity in the present example allows for the occurrence of arbitrarily long sequences of S ! I ! R events, being the computational situation more involved.
The generating function reads, With these elements we can compute Hence, The remaining rates are linear, hence r I 2 ¼ c and r R 3 ¼ d: As for the differential equations, The stochastic number of events occurring in a short time-interval h can be obtained as follows: (1) Solve 16 and obtain the coefficients P k ðhÞ.
(3) Compute the stochastic number of following events. For example, from p ¼ P SI 12 δ I I P S 1 ¼ P SI 12 P S 1 the following events distribute as ð1 À pð1 À z 2 ÞÞ n1δ I 1 ¼ ð1 À pð1 À z 2 ÞÞ n1 . Similar equations hold for the other two cases, modifying the indices cyclically.

Comments
Let us assume that we have the SIRS system S ! ri I ! rr R ! rs S and our population values are S; I; R ð Þ¼ S 0 ; 1; R 0 ð Þ . In case, the race between recovery of infected and new infection (i.e., the next event influencing I) is won by the recovery, the epidemic stops at that point, and the infected population becomes extinct. We know then that, with this initial condition, extinction is bounded by P I¼0 f g > rr rrþS 0 ri . We know as well that the probability of extinction increases with time. We further know that the probability of extinction in a time-interval Dt is bounded as P I¼0 f g ðDtÞ > rr rrþS 0 ri 1 À expðÀðrr þ S 0 riÞDtÞ ð Þ to a very good approximation (the parenthesis expresses the probability of the occurrence of an event in ½0; Dt).
We now consider the case in our approximation up to first order. Since, IðDtÞ ¼ 1 À n 2 þ n 1 and n 2 1 we need n 1 ¼ 0 and n 2 ¼ 1.
The probability Pðn 1 ¼ 0Þ ¼ ð1 À P 1 Þ S 0 ¼ expðÀriS 0 DtÞ must be multiplied by the probability Pðn 2 ¼ 1Þ ¼ P 2 ¼ 1 À expðÀrrDtÞ since both events are assumed to be independent in the approximation. Thus, P I¼0 f g ¼ expðÀriS 0 DtÞ 1 À expðÀrrDtÞ ð Þ . This function is convex and it is zero at t ¼ 0 while also lim t!1 P I¼0 f g ðtÞ ¼ 0. Hence, unlike the exact case, the probability is not non-decreasing.
The observation imposes a qualitative limit for the possible time steps, namely Dt ( logð rr ri S 0 þ 1Þ=rr. For example, S 0 ¼ 6000, rr ¼ 0:8, ri ¼ 0:0002 we have Dt ( 0:64 Another interesting feature that generalises to all systems is that the expression < S > that appears in the denominator of r ¼ ri < SI > < S > can be written as < S > ¼ < S 0 À n 1 þ n 3 > ¼ S 0 þ À @ψ @z 1 þ @ψ @z 3 n o z¼1 and has the form < S > ¼ S 0 a 1 þ R 0 a 2 þ I 0 a 3 . The factors a i are Oðh iÀ1 Þ and according to the approximation lemma (in Appendix A), they admit a power-series expansion. If only the lowest orders are to be kept, it is needed that rrDt; riDt; rsDt f g ( 1. In the case of our example, considering S 0 ; I 0 ; R 0 f g%10 4 this requires Dt < 0:01, but the requirement to neglect R 0 a 2 in front of S 0 a 1 in all cases require to consider the worst scenario S 0 ¼ 1; R 0 ¼ 10 4 . This translates into Dt ( 10 À4 . Thus, further approximations should be considered with caution and it is advisable not to perform them. In the same form, < SI > will have 12 terms when parsed by the first and second powers of S 0 ; I 0 ; R 0 f g .

Simulations
We use the SIRS model to illustrate the differences between different simulations. The model corresponds exactly to a Feller-Kendall process; hence, we visually compare the exact process, the deterministic model and the present approach.
In Figure 1, it is shown that the approximated method closely follows the exact simulation in contrast with the deterministic method that presents large departures after reaching the minimum number of infected individuals (upper right panel in the figure). Both the total number of runs (out of 10 4 ) where the infected population goes into extinction and the distribution in time of these events is well represented by the approximation method (no extinction is possible in the deterministic simulation).

Concluding remarks
The goal of this manuscript has been to discuss, develop and implement a quasilinear approximation to the Feller-Kendall algorithm for Markov Jump processes. The work follows closely the ideas in Solari and Natiello (2014), where linear and constant transition rates are discussed. In the present work, we extend those ideas to the more useful situation where rates are not necessarily linear andquite importantthey are also time-dependent. The relevance of the use of approximations ultimately lies in the judgement of the user. However, for population problems with realistic aspirations, the Feller-Kendall algorithm soon becomes too slow and reliable approximations are a reasonable way out. The alternatives available for choosing a proper approximation level are therefore also discussed. Lemma 5. For t 2 ½0; h, h ¼ t À t 0 > 0 sufficiently small and for each order of approximation