Approximations of Cumulants of the Stochastic Power Law Logistic Model

Asymptotic approximations of the first three cumulants of the quasi-stationary distribution of the stochastic power law logistic model are derived. The results are based on a system of ODEs for the first three cumulants. We deviate from the classical moment closure approach by determining approximations without closing the system of equations. The approximations are explicit in the model’s parameters, conditions for validity of the approximations are given, magnitudes of approximation errors are given, and spurious solutions are easily detected and eliminated. In these ways, we provide improvements on previous results for this model.


Introduction
We study a stochastic version of the power law logistic model. Logistic models are used as models for the size of a population with density-dependent growth. This means that the net growth rate per individual is a decreasing function of the population size. The simplest decreasing function is the linear one. It leads to the classical logistic model, whose deterministic version was studied by Verhulst already in 1838. The stochastic version of this model has been studied extensively, as shown by Nåsell (2017). More general deterministic power law logistic models have been discussed by Banks (1994) and by Tsoularis and Wallace (2002), while stochastic versions of such models are treated by Matis et al. (1998), by Renshaw (2011), and by Bhowmick et al. (2016).
The model that we deal with is a Markov Chain with finite state space and continuous time, and with an absorbing state at the origin. It can also be described as a finitestate birth-death process. Two model formulations are given in Sect. 2. They are B Ingemar Nåsell ingemar@kth.se 1 Department of Mathematics, The Royal Institute of Technology, 100 44 Stockholm, Sweden mathematically equivalent, but with different parameters. The maximum population size N is an important parameter in the second one of the two formulations. This formulation is basic for our work, since N serves the role of being the large parameter for which asymptotic approximations of various quantities can be derived.
The state variable X (t) of either of the two formulations is interpreted as the number of individuals in a population. Ultimate extinction of this variable is a fact for this model. This means in particular that the stationary distribution of X (t) is degenerate with probability one at the origin. A related and important random variable X Q (t) is defined by conditioning X (t) on non-extinction. The stationary distribution of this conditioned random variable is the so-called quasi-stationary distribution (QSD). The QSD is important for the information it gives about the long-term behavior of a surviving population. Our main goal is to derive approximations of the first three cumulants of the QSD.
Our approach toward this goal resembles the well-known moment closure method in the sense that both are based on a system of ODEs for the first few (three) cumulants. This system of ODEs is not closed: the number of unknown functions (cumulants) is larger than the number of equations. This fact has caused all practitioners of moment closure to introduce an approximating assumption that leads to closure of the system of ODEs. The method is used widely, as shown by the recent review paper by Kuehn (2016). However, it is also known to possess several weaknesses. The first of these weaknesses is that no condition for validity of the resulting approximation is available, a second one is that the magnitude of the approximation error can not be evaluated, a third one is that it often leads to spurious solutions that require large efforts to eliminate, and a fourth one is that the dependence of the approximating expressions on the model's parameters is not known. An early basic paper in the area was written by Whittle (1957). The desirability of closing the system of ODEs for the cumulants cannot be denied. It is indeed a necessity if the goal is to derive exact expressions for the cumulants. However, it is important to realize that closure is not necessary if one is satisfied with approximate results, as we are in this paper.
We do not use moment closure methods at all. Instead, we apply the alternative mathematical method introduced in Nåsell (2017). Thus, we base our work on a system of ODEs for cumulants, just as is done by practitioners of moment closure. However, we avoid the step that leads to closure of the system of ODEs, since the ad hoc nature of this step is the cause of several weaknesses of the moment closure method. Our alternative is to search for asymptotic approximations. This requires the reparameterization that accompanies the second one of the two model formulations of Sect. 2. To proceed with our approach, we also need information about the orders of magnitude (in terms of N ) of the cumulants that appear in the system of ODEs. We motivate our assumptions in this regard with numerical evaluations. The method that we use avoids the weaknesses associated with moment closure methods.
Methods for numerical evaluations of the QSD and of the associated cumulants are discussed in Sect. 3. Our approach is different from what has been used by previous workers on this model. Results of numerical evaluations are used to motivate forms of basic assumptions that are made for determining asymptotic approximations of the first few cumulants of the QSD and also for checking the orders of magnitude of the errors of the approximations that we derive.
Sections 4 and 5 study the two random variables X (t) and X Q (t), respectively. We give a system of first-order ODEs for the first 3 cumulants of the unconditioned random variable X (t) in Sect. 4. The same system of equations is shown in Sect. 5 to provide an approximation for the system of ODEs of the first 3 cumulants of the conditioned random variable X Q (t) in a particular parameter region. This means that approximations of the cumulants of the QSD in this region are found as critical points of the system of ODEs. Asymptotic approximations of the coordinates of the critical points of this system of ODEs are derived in Sect. 6. They serve to give asymptotic approximations of the first 3 cumulants of the QSD. These results hold for small positive integer values of the parameter s that describes the power of the population size n that gives the decreasing function of the population size that reflects the density dependence of the net birth rate of the model. An extension of these results to arbitrary positive values of s, both integer and non-integer, is proposed in Sect. 7. Comparisons are made in Sect. 8 with published results that also deal with the cumulants of the QSD of the same model. The paper ends with some concluding comments in Sect. 9.

Model Formulation
The stochastic power law logistic model is a model for the size of a population with density-dependent growth. It is formulated as a birth-death process {X (t), t ≥ 0}. The hypotheses of the model are summarized in the descriptions of the population birth rate λ n and the population death rate μ n as functions of the state n of the process. The formulation given by Matis et al. (1998) takes the following form: the population birth rate equals and the population death rate is The state space appears to be unbounded and equal to {0, 1, 2, . . . }. The parameters a 1 , b 1 , a 2 , b 2 , s are all assumed to be positive. This model formulation takes its form from the similar model introduced by Bartlett (1957) for the special case of the Verhulst logistic model where s = 1. We shall use a different formulation, where The state space of the process in this formulation is finite and equal to {0, 1, 2, . . . , N }. The parameter space contains the five parameters N , R 0 , α, μ, and s. Among these, N is a large positive integer that represents the maximum population size, R 0 is a positive dimensionless threshold parameter, α is a nonnegative dimensionless parameter, μ is a positive death rate with dimension inverse time, and s is a positive dimensionless number. This model formulation is an extension from the formulation of the logistic Verhulst model, corresponding to s = 1, given by Nåsell (2017). An advantage of this formulation over the previous one is that dimensions and biological interpretations of the parameters are known. Under the assumption that the initial distribution is supported on the state space, the process remains there for all time, since μ 0 = 0 and λ N = 0. The origin is seen to be an absorbing state, since both λ 0 and μ 0 are equal to zero. Absorption at the origin corresponds to extinction of the population. To study a surviving population, we introduce the conditioned random variable X (Q) (t) by conditioning X (t) on the event X (t) > 0. Thus, The state space of the conditioned random variable X (Q) (t) differs from that of X (t) in one respect: The origin can be reached by X (t), but not by the conditioned random variable X (Q) (t). Thus, the state space of the latter one of these random variables is equal to {1, 2, . . . , N }. Stationary distributions of the two random variables are vastly different. The stationary distribution of X (t) is degenerate with probability one at the origin, while the stationary distribution of the conditioned random variable X (Q) (t) is the important quasi-stationary distribution (QSD), supported on the state space {1, 2, . . . , N }. We shall use the second of the two model formulations in this paper. An important reason for this is that the second formulation contains a parameter N that can take large values. The presence of such a parameter is essential for the formulation of asymptotic approximations. Another advantage of the second formulation is that it gives information about dimensions of parameters and state variables. Such information is of value for biological interpretations of all quantities that appear in the study of the model.
The two formulations of the stochastic power law logistic model are essentially equivalent if we require in the first formulation that (a 1 /b 1 ) 1/s = N is an integer, that the state space is finite and equal to {0, 1, 2, . . . , N }, and that b 2 ≥ 0. We assume in what follows in this paper that these three requirements are met. The second model formulation can then be seen as a reparameterization of the first one. We find then that the four parameters a 1 , a 2 , b 1 , b 2 of the first formulation can be expressed in terms of the five parameters μ, R 0 , α, N , s of the second formulation as follows: Obviously, s takes the same value in both formulations. We note also that some results that we shall refer to below use the following notations for the sums and differences of the parameters a i and b i : As already mentioned, the logistic Verhulst model can be seen as a special case of the power law logistic model, corresponding to s = 1. Asymptotic approximations of the first three cumulants of the QSD for the Verhulst model are derived in Nåsell (2017). The goal of the present paper is to derive corresponding approximations of the first three cumulants of the QSD of the power law logistic model for positive values of s. We proceed by first deriving asymptotic approximations of the first three cumulants of the QSD for the s-values 2, 3, and 4. These results are then extended in Sect. 7 to both integer and non-integer positive values of s.

Numerical Evaluations
Numerical evaluations have been used in earlier work on this model, and we shall also use them here. However, there are differences between our approach and those of earlier workers both in the methods used for deriving numerical results, and in the type of evaluations that are carried out. Bartlett et al. (1960) argue incorrectly that a stationary distribution cannot exist for the model that they are concerned with, since λ 0 = 0, while we claim that the process in this case has a degenerate stationary distribution with probability one at the origin. Bartlett et al. claim furthermore that "the probability distribution for the stationary (or quasi-stationary) distribution must satisfy the recurrence relation" μ n p n = λ n−1 p n−1 .
We find it disturbing that they refer to both stationary and quasi-stationary distributions here, without clarifying which of these two distributions that they mean. The statement can certainly not hold for both of them. The fact is that it is incorrect for both of them. Instead, the relation (14) holds for the stationary distribution p (0) of a related auxiliary process X (0) (t) that is useful for studying the quasi-stationary distribution. The transition rates of this auxiliary process are equal to those of the process X (t), with the one exception that the rate μ 1 for transition from the state 1 to the absorbing state 0 is replaced by zero. The stationary distribution p (0) is non-degenerate and can be evaluated explicitly. It has been used to study the QSD for a SIS model, which is the special case of the model that we are concerned with here that corresponds to the parameter values s = 1 and α = 0, in Nåsell (2011). It follows from this study that p (0) is a good approximation of the QSD in its body, but not in the left tail, in the parameter region where R 0 > 1, but also that p (0) is not a good approximation of the QSD if R 0 < 1, nor in the parameter region near R 0 = 1. Similar relations between p (0) and the QSD are expected to hold for α > 0 and for integer values of s larger than 1. Both Bartlett et al. (1960) and Matis et al. (1998) use the recurrence relation (14) as a basis for numerical evaluations of the QSD. Our approach is different, as shown below. However, we expect the cumulants of the stationary distribution p (0) to be close to the cumulants of the QSD in the parameter region where R 0 > 1, since the tail probabilities have a minor influence on the cumulant values. We note that Renshaw (2011) uses p (0) as a "working definition" of quasi-stationarity. We do not follow this approach. Our results require a sharp line to be drawn between exact and approximate relations.
To describe our method for numerical evaluations, we need a second auxiliary process X (1) (t). The birth rates of this process are equal to those of the original process X (t), while the death rates are slightly smaller than those of X (t), and equal to μ n−1 . This means that the second auxiliary process can be interpreted as the original process X (t) with one surviving and immortal individual. Like the first of the two auxiliary processes, the second one lacks absorbing state, and its stationary distribution p (1) is non-degenerate and can be evaluated explicitly.
Our numerical evaluations of the QSD make use of the restart map Ψ analyzed by Ferrari et al. (1995), and also discussed by Nåsell (2011). This restart map is defined as follows: Let a probability vector ν = (ν 1 , ν 2 , . . . , ν N ) be given. Define a process related to the process X (t) that we are studying by the requirement that whenever the original process reaches the state zero, it is immediately restarted at some state j with probability ν j . The restarted process has the state space {1, 2, . . . , N }, and a unique stationary distribution p. The map Ψ is then defined by Ψ (ν) = p. The quasi-stationary distribution q is a fixed point of the mapping Ψ , since Ψ (q) = q. Furthermore, Ferrari et al. (1995) show that an iteration scheme can be established by repeated applications of the map Ψ to an arbitrary initial probability vector, and that it converges to the quasi-stationary distribution q. Suitable initial probability vectors are given by the stationary distributions p (0) and p (1) of the two auxiliary processes X (0) (t) and X (1) (t). It turns out that that the map Ψ can be described explicitly with the aid of these two stationary distributions. A derivation is given by Nåsell (2011). By denoting the components of the vector p = Ψ (ν) by p n , we find that p n = π n S n p 1 , n = 1, 2, . . . , N , where π 1 = 1, π n = λ 1 λ 2 · · · λ n−1 μ 2 μ 3 · · · μ n , n = 2, 3, . . . , N , and with ρ 1 = 1, ρ n = λ 1 λ 2 · · · λ n−1 μ 1 μ 2 · · · μ n−1 , n = 2, 3, . . . , N , and The stationary distributions p (0) and p (1) of the two auxiliary processes are determined from the sequences π n and ρ n as follows: and p (1) n = ρ n p (1) Our numerical method for determining the quasi-stationary distribution consists in applying the restart map Ψ to a suitable initial distribution and continues iterations until successive iterates are sufficiently close. In case R 0 > 1, which is the parameter region of main interest for the study of cumulants in this paper, the stationary distribution p (0) is recommended over the distribution p (1) as initial distribution.
Numerical evaluations are used by Matis et al. (1998) for showing that the errors that their approximations lead to are small. With our different parameterizations, we have access to a parameter N that can take large values. We can then use numerical evaluations to extract additional information both about the model and about the approximations that we derive. One goal is to derive information about the orders of magnitude (in terms of N ) of various cumulants, and another one concerns the probability of taking the value 1 in the quasi-stationary distribution. Information of these two types is important for the formulation of hypotheses as a basis for finding asymptotic approximations of the first few cumulants for large values of N , as will be shown below. Another use of numerical evaluations is to derive information about the orders of magnitude (in terms of N ) of the errors that are committed in using our approximations. Such results will be used below to support the forms of asymptotic approximations that we derive.

ODEs for Cumulants of the Unconditioned Random Variable X(t)
The starting point for our derivations of asymptotic approximations of the first three cumulants of the QSD is ODEs for the first three cumulants of the conditioned random variable X Q (t). It turns out that these ODEs are closely related to the ODEs for the corresponding cumulants of the unconditioned random variable X (t) in one important region of parameter space, namely where R 0 > 1. We start therefore by giving these latter ODEs.
The derivation of these ODEs makes use of the moment generating function M(θ, t) and the cumulant generating function K (θ, t) of the random variable X (t). They are defined by and where the state probabilities p n (t) = P{X (t) = n} satisfy the Kolmogorov forward equations p n (t) = λ n−1 p n−1 (t) − (λ n + μ n ) p n (t) + μ n+1 p n+1 (t), n = 0, 1, . . . , N .
(Here, we agree to put λ −1 = μ N +1 = p −1 (t) = p N +1 (t) = 0 so that (24) makes sense formally for all n values indicated.) The definitions of the cumulant generating function K and of the transition rates λ n and μ n can be used to derive a partial differential equation (PDE) for K . It can be written as follows: The same result follows by applying the reparameterization in (6)-(9) to the expression for the PDE given by Renshaw (2011), and numbered (3.4.11). The cumulants κ i (t) are obtained from a power series expansion of the cumulant generating function K (θ, t) as follows: We use this result for determining the ODEs of the first few cumulants of the random variable X (t). They are found by determining the first few terms of the expansion of the PDE of the cumulant generating function K (θ, t) in terms of θ and then identifying terms with equal powers of θ .
To express our results, we use the capital letters A, B, and C to denote the time derivatives of the first three cumulants of the random variable X (t), with superscripts indicating the s-values, as follows: Expressions for the functions A (s) , B (s) , C (s) for the s-values 2, 3 and 4 can then be written as follows: Derivations of these results using Maple are given in Nåsell (2018). We note that the 3 derivatives A (2) , A (3) , A (4) in (30), (33), (36) of the first cumulant κ 1 for the 3 s-values 2, 3, and 4 are found from the expressions (27), (28), (29) in Matis et al. (1998) by reparameterization using (10)-(13) after introducing the correction that the first term (a − bκ 4 1 ) in the right-hand side of (29) is written (a − bκ 4 1 )κ 1 . Similarly, the 2 derivatives B (2) and B (3) in (31) and (34) of the second cumulant κ 2 for the 2 s-values 2 and 3 are found from the expressions (33) and (34) in Matis et al. by the same reparameterization after noting that the minus sign of the term −2aκ 2 in the right-hand side of their formula (34) is incorrect and should be changed to a plus sign.

ODEs for Cumulants of the Conditioned Random Variable X (Q) (t)
The goal of our study is to determine approximations of the first three cumulants of the QSD for the stochastic variable X (t). This leads us to a study of the stationary values of the system of ODEs for the first three cumulants of the conditioned random variable X (Q) (t). These ODEs turn out to be closely related to the ODEs for the cumulants of the unconditioned random variable X (t) in one important parameter region, namely where R 0 > 1. Identification of this parameter region is of high importance, since it gives the validity condition for the approximations that we derive.
In similarity to the case in the previous section, we use the cumulant generating function K (Q) (θ, t) of the random variable X (Q) (t) that is of concern here to derive the ODEs of the first few cumulants of this random variable. To derive the PDE of this cumulant generating function, we proceed as in the previous section. The main difference from the case in the previous section is that the expression for the time derivative of the state probability p (Q) n (t) is different from the counterpart of (24). By using the relation in (5), we find that Differentiation and use of the Kolmogorov forward equations in (24) gives The PDE for the cumulant generating function K (Q) (θ, t) of the conditioned random variable X Q (t) turns out to be quite similar to the PDE for K (θ, t). It can be written as follows: As in the previous section, we determine the ODEs of the first few cumulants by expanding this PDE in terms of θ and identifying terms of equal powers of θ . We find then that the last term of the right-hand side of every such ODE contains p (Q) 1 (t) as a factor. The stationary value of this probability is equal to the quasi-stationary probability q 1 . We claim that this probability is exponentially small in N for any α ≥ 0 and any positive value of s, integer or not, when R 0 > 1. Arguments that support this are given below. The last term can therefore be ignored when we are searching for asymptotic approximations of the critical points. An important consequence of this is that asymptotic approximations of the first few cumulants of the QSD for R 0 > 1 are found as stationary solutions of the ODEs in Sect. 4 for the corresponding cumulants of the unconditioned random variable X (t), with obvious rules for exclusion of spurious solutions. Derivations of these results are given below in Sect. 6.
We show first that q 1 is exponentially small in N when R 0 > 1, α ≥ 0, and s = 1. It follows from Nåsell (2011) that this result is true for the SIS model, defined by s = 1 and α = 0. Furthermore, it follows from Nåsell (2001a, b) that the same result holds also for the Verhulst model, with s = 1 and α > 0. Support for the wider claim that q 1 is exponentially small in N for R 0 > 1, α ≥ 0, and any s > 0 is given in Table 1, Results are shown for R 0 = 3, α = 1, the s-values 0.5, 1, 2, 3, and the N -values 100, 200, and 400. The results indicate that q 1 is exponentially small. Derivations of these results using Maple are given in Nåsell (2018) which shows numerically determined values of q 1 for R 0 = 3, α = 1, the s-values 0.5, 1, 2, 3, and the N -values 100, 200, and 400. The table shows that q 1 decreases strongly with both s and N for these values of R 0 and α. It shows also that a doubling of the N -value for fixed values of R 0 , α, and s leads to an approximate squaring of the value of q 1 , as would be expected when q 1 is exponentially small in N .
The condition that we have found under which q 1 is exponentially small in N , namely R 0 > 1, turns out to be the condition for the validity of the approximations of the first three cumulants of the QSD that we derive later. We note that a corresponding condition of validity has been missing from earlier results for the same model, based on cumulant closure.

Cumulants of the QSD: First Results
This section is used to derive asymptotic approximations of the first three cumulants of the quasi-stationary distribution of the stochastic power law logistic model for the s-values 2, 3, and 4 and large N -values in the parameter region where R 0 > 1. The method that we use is similar to the one introduced in the study of the cumulants of the QSD of the Verhulst logistic model in Nåsell (2017). This latter study actually represents the special case with s = 1 of the model that we study here. Our results are based on assumptions about the forms of asymptotic approximations of the first few cumulants for large values of N . For the Verhulst model with s = 1 and R 0 > 1, we give arguments in Nåsell (2017) for assuming that the first four cumulants of the QSD are O(N). The results make it easy to conjecture that all finite order cumulants of the QSD for s = 1 are O(N ) in the parameter region R 0 > 1. It is tempting to guess that this holds also for integer values of s larger than 1.
In what follows in this section, we shall assume that the first s + 3 cumulants of the QSD are O(N) when s = 1, 2, 3, 4, α ≥ 0, and R 0 > 1. As a basis for this, we refer to the numerical results in Table 2. It shows numerically determined values of the first seven cumulants of the QSD for the parameter values R 0 = 10 and α = 1, with the two s-values 1 and 4, and with the three N -values 100, 200, and 400. It is easy to verify from the table that all cumulants experience an approximate doubling both when N is doubled from 100 to 200 and also when N is doubled from 200 to 400. We take this as a strong indication that all of the cumulants considered are O(N ). The results in Table 2 are derived using Maple in Nåsell (2018), where also similar results are given for the s-values 2 and 3.
We study first the case s = 2. The assumptions that the first five cumulants κ 1 −κ 5 are O(N ) in the parameter region R 0 > 1 for s = 2 are basic for our further assumptions that the first five cumulants have the following asymptotic behaviors for R 0 > 1 and s = 2: The reason for including different numbers of terms in these asymptotic approximations will become apparent shortly. Thus, we have introduced eight unknowns, namely x 1 , x 2 , x 3 , y 1 , y 2 , z 1 , u 1 , u 2 . We proceed to determine the first six of them. By inserting the above asymptotic approximations of the first five cumulants into the expressions (30)-(32) for the functions A (2) , B (2) , C (2) , we find that asymptotic approximations for them can be written as follows: where The reason for including different numbers of terms in the assumed asymptotic approximations of the cumulants κ 1 − κ 5 can be understood with reference to the resulting asymptotic approximations of the quantities A (2) , B (2) , and C (2) . We note for example from the expression (30) for A (2) that an additional term in the assumed asymptotic approximation of κ 2 would contribute a term of the order of 1/N 2 to the asymptotic approximation of A (2) . Inclusion of such an additional term would therefore be absorbed in the error term for the asymptotic approximation of A (2) . We note from the expressions in (50)-(55) that only the first six of the eight coefficients introduced above appear in these expressions. Our assumptions for s = 2 that the first five cumulants have the asymptotic approximations given in (42)-(46) lead to the asymptotic approximations of A (2) , B (2) , C (2) in (47)-(49). These expressions are equal to the derivatives with respect to time of the first three cumulants. Setting them equal to zero gives conditions for the stationary points of the corresponding system of ODEs. Since we are working with approximations instead of exact results, we transform these conditions into the three requirements that . These requirements are satisfied by the basic conditions that the six expressions 1 in (50)-(55) are equal to zero. We note that each of these expressions is a polynomial of degree three in the six unknown coefficients x 1 , x 2 , x 3 , y 1 , y 2 , z 1 . Our basic problem is to solve the six equations that are formed by setting the corresponding expressions equal to zero for the six unknown coefficients. Each equation is a polynomial of degree s + 1 = 3 in the six unknown coefficients. Solution appears possible, since the number of equations equals the number of unknowns. Further inspection of the six equations reveals that considerable simplification can be achieved in the solution by determining the six unknowns in a definite order, as described below. The first advantage is that each of the equations to be solved contains only one unknown coefficient, and the second one is that all equations except the first one are linear in the coefficient to be solved for. The very first equation is of degree s + 1 = 3 in the case we are considering here with s = 2. All of the roots of the first equation except one are spurious solutions that turn out to be easy to identify and to delete. Whenever a solution is found for one of the coefficients, its value is immediately inserted into the remaining unsolved equations. The order of solution that leads to these pleasant results is as follows: First, the equation A (2) 1 = 0 is solved for x 1 . This equation has in this case, with s = 2, s + 1 = 3 solutions. Among them, we exclude x 1 = 0 and x 1 < 0 as the only spurious solutions that appear in this method. The solution found for x 1 is then inserted into the remaining five unsolved equations. After finding x 1 , we proceed to solve the equation B (2) 1 = 0 for y 1 , the equation A (2) 2 = 0 for x 2 , the equation C (2) 1 = 0 for z 1 , the equation B (2) 2 = 0 for y 2 , and the equation A (2) 3 = 0 for x 3 . It is elementary to use these rules to solve for the six unknown coefficients. Alternative derivations using Maple are given in Nåsell (2018). The solutions are as follows: The method used here for derivation of asymptotic approximations of the first three cumulants of the QSD for s = 2 can in principle be followed for higher integer values of s. We proceed to consider the case when s = 3. We begin by assuming that the first six cumulants have the following asymptotic behaviors for R 0 > 1 and s = 3: Insertions of these asymptotic approximations of the first six cumulants into the expressions (33)-(35) for the functions A (3) , B (3) , C (3) lead to the following asymptotic approximations for them: where These six expressions are all polynomials of degree four in the six unknowns x 1 , x 2 , x 3 , y 1 , y 2 , z 1 . We solve the six equations formed by putting each of these expressions equal to zero for the six unknowns x 1 , x 2 , x 3 , y 1 , y 2 , z 1 . First, the equation A 1 = 0 is solved for x 1 . Among the four roots, we exclude the one that is equal to zero, and the two that are complex conjugate as the only spurious solutions that appear in this method. The remaining five equations are then solved for the remaining five unknowns, using an order of solution similar to what is described for the case s = 2 above. Derivations using Maple are given in Nåsell (2018). The results are as follows: x 2 = − 2 3 The last problem addressed in this section is the derivation of asymptotic approximations of the first three cumulants of the QSD for s = 4. We begin by assuming that the first seven cumulants have the following asymptotic behaviors for R 0 > 1 and s = 4: Insertions of these asymptotic approximations of the first seven cumulants into the expressions (36)-(38) for the functions A (4) , B (4) , C (4) lead to the following asymptotic approximations for these three functions: where A (4) These six expressions are all polynomials of degree five in the six unknown quantities x 1 , x 2 , x 3 , y 1 , y 2 , z 1 . As above, we solve the six equations formed by putting each of these expressions equal to zero for the six unknowns x 1 , x 2 , x 3 , y 1 , y 2 , z 1 . First, the equation A (4) 1 = 0 is solved for x 1 . Among the five roots, we exclude four of them as spurious solutions, namely one that is equal to zero, one that is negative, and two that are imaginary. The remaining five equations are then solved for the remaining five unknowns, using the same order as above for the cases s = 2 and s = 3. Derivations using Maple are given by Nåsell (2018). The results are as follows:

Cumulants of the QSD: Extensions to Positive Values of s
The results derived in the previous section give asymptotic approximations of the first three cumulants of the QSD for the integer s-values 2, 3, and 4. The present section is used to extend these approximations to positive values of s, both integer and noninteger. The extension to positive integer values of s starts out by giving asymptotic approximations of the first three cumulants of the QSD for the integer s-values from 1 to 10. For this, we use the method spelled out in detail in the previous section. It turns out that the s-dependence of these results can be described with the aid of five functions of s that all can be expressed as ratios of polynomials in s of degree two. After these five functions have been determined, we indicate that they can be used to extend our results to arbitrary positive integer values of s. No further evaluations are necessary for the extension of these results to positive non-integer values of s. We argue that small changes in the parameter s will lead to small changes in the first few cumulants of the QSD. The corresponding continuity argument lies behind our claim that the results that have been derived for positive integer values of s are valid for all positive values of s. For easy reference we quote first from Nåsell (2017), the expressions for the six coefficients x 1 , x 2 , x 3 , y 1 , y 2 , z 1 in the case s = 1: The asymptotic approximations of the first three cumulants of the QSD for the stochastic power law logistic model are given in terms of the following expressions for the coefficients x 1 , x 2 , x 3 , y 1 , y 2 , z 1 for any s > 0:

I. Nåsell
where the functions h 1 − h 5 are defined as follows: h 4 (s) = s 2 + 12s + 11 24s 2 , (120) To show that the functions h 1 − h 5 are determined by these expressions for s > 0, we determine first the values that they take for the integer s-values 1-10, and given in Table 3. The values taken by these five functions for the integer s-values 1-4 are readily found from the expressions found for the 6 coefficients x 1 , x 2 , x 3 , y 1 , y 2 , z 1 given in the previous and present sections. The further results for the s-values from 5 to 10 are based on asymptotic approximations of the first three cumulants of the QSD of the power law logistic model that have been derived using Maple, and are reported in Nåsell (2018).
It is straightforward to use the entries in Table 3 for the functions h 1 , h 2 , and h 5 to confirm that the expressions (117), (118), (121) for these functions are valid for all integer values of s from 1 to 10.
We turn now to deal with the functions h 3 and h 4 . To find expressions for them, we assume that each of them can be written as a quotient of two polynomials in s of degree 2. To be specific, we assume that To determine the values of the six coefficients a 0 , a 1 , a 2 , b 0 , b 1 , b 2 , we establish six equations that are found by equating the values of h 3 (s) for the integer s-values from 1 to 6 according to (122) to the corresponding values in Table 3. By solving these equations using Maple, as seen in Nåsell (2018), we find that h 3 (s) equals the expression in (119). An entirely similar treatment of the function h 4 leads to the expression (120) for h 4 (s). So far, we can conclude that the expressions (119) and (120) (119) and (120) and the corresponding values in Table 3. We take this as a strong indication that the domain of validity of the expressions (119) and (120) can be further extended to all positive integer values of s. As a last step of extension, we use a continuity argument to claim that the domain of validity for the five functions h 1 − h 5 can be extended from all positive integers s to all positive values of s. We illustrate the result for non-integer s-values by giving the numerical values of the error terms for the three cumulants κ 1 , κ 2 , κ 3 in case R 0 = 10 and α = 1, for the two s-values 0.5 and 3.5, and the 3 N -values 100, 200, and 400 in Table 4. The table indicates that the error term for κ 1 is divided by approximately 4 for each doubling of N , the error term for κ 2 is divided by approximately 2 for each doubling of N , and the error term for κ 3 is approximately constant when N is doubled. This is consistent with the claims that the error terms of κ k are O(1/N 3−k ) for the k-values 1, 2, and 3. Evaluations of the error terms in Table 4 are derived using Maple in Nåsell (2018).
The results in (111)-(116) show the way in which our asymptotic approximations of the first three cumulants of the QSD of the stochastic power law logistic model depend on the four parameters N , R 0 , α, and s. The one-term asymptotic approximations of the first three cumulants κ 1 , κ 2 , and κ 3 are useful for the information that they give about the behaviors of these cumulants. They can be written as follows: where expressions for the three coefficients x 1 , y 1 , z 1 are given in (111), (114), and (116), respectively. We use these expressions to derive some properties of the first three cumulants, valid for sufficiently large values of N .  We note first that the expectation κ 1 is an increasing function of the power s. This follows from the following expression for the derivative of x 1 with respect to s:

I. Nåsell
Our second observation concerns the variance κ 2 . For sufficiently large values of N , we find that it has a maximum as a function of s at s = s 2 , where This follows from the following expression for the derivative of y 1 with respect to s: We conclude in particular that κ 2 is increasing in s for s < s 2 and decreasing in s for s > s 2 , provided N is large enough. Our third observation concerns the third cumulant κ 3 . It is useful for determining the skewness of the QSD. Actually, the QSD has negative skewness if κ 3 is negative, and positive skewness if κ 3 is positive. It follows from the expression (116) that z 1 is positive if s < s 3 and negative if s > s 3 , where It is easy too see that the same intervals of s lead to positive, respectively negative, skewness of the QSD if N is sufficiently large.
Our brief study of the s-dependence of the first three cumulants shows that one can identify two cases with rather different behaviors. The first case occurs when s < min(s 2 , s 3 ), and the second one when s > max(s 2 , s 3 ). In the first case, we conclude that both κ 1 and κ 2 are increasing functions of s, and that the QSD has Fig. 1 The QSD for the stochastic power law logistic model with N = 100, R 0 = 5, α = 1, and the 5 s-values 0.2, 0.5, 1, 3, and 10 positive skewness. In the second case, we conclude as in the first case that κ 1 is increasing in s, while κ 2 now is a decreasing function of s and the skewness of the QSD has changed from positive to negative.
We illustrate the s-dependence of the first three cumulants by plotting the QSDs for five different s-values with the constant parameter values N = 100, R 0 = 5, and α = 1. In these cases, we find s 2 = 0.4055 and s 3 = 0.3846. Numerically determined QSDs are shown in Fig. 1 for the five s-values 0.2, 0.5, 1, 3, and 10. The figure shows that the first cumulant (the expectation) increases as a function of s, that the second cumulant (the variance) increases as a function of s for small s, but decreases as a function of s for large s, and finally that the third cumulant (a measure of skewness) is positive for small s-values, and negative for large s-values. (For the interpretations of the plots in Fig. 1, it is useful to note that the probabilities q n are positive for all n in the state space {1, 2, . . . , N }. Recall that the QSDs are discrete, and that the individual probabilities q n are determined from the plotted curves by reading off the values at each abscissa n.) The asymptotic approximations that we have given of the first three cumulants are valid for sufficiently large N -values. It is seen from the plot in Fig. 1 that the requirement that q 1 is exponentially small is not satisfied for s = 0.2 when R 0 = 5, N = 100 and α = 1. Clearly, a larger value of N is needed here for q 1 to be exponentially small in N .

Discussion of Other Approaches
In this section, we discuss five published results that all deal with the problem in this paper, namely to determine approximations of the first three cumulants of the QSD of the stochastic power law logistic model. All these results are formally different from each other and also different from the results that we have derived and presented in Sect. 7 of this paper. We use the powerful method of asymptotic approximations to analyze these results. In this discussion, we are mainly limited to the case s = 1.
The first published results are those due to Bartlett (1957) and Bartlett et al. (1960). The latter results are referred to as the BGL approximations. The model they analyze is the first one formulated in Sect. 2. Their derivation is based on a study of the changes of the stochastic variable X (t), and of its second and third powers, during an infinitesimal time interval. This requires approximations for which validity and error magnitudes have not been analyzed. The resulting approximations of the first three cumulants of the QSD are denoted by κ (B) 1 , κ (B) 2 , and κ (B) 3 , respectively. They are given as follows in Bartlett et al. (1960): In the quoted paper, these three approximations are denoted by the sign ∼ instead of ≈. This would indicate that the approximations have the desirable property of being asymptotic under the condition that some identified variable becomes large. However, no such variable appears in the derivations. Because of this, the implicit claim by Bartlett et al. (1960) that the approximations are asymptotic has not been proved by them. It is, however, easy to use the results that we have derived to show that the BGL approximations of the first three cumulants are asymptotic. All one has to do is to use the reparameterization in (6)-(9). This leads to the following results: where x 1 , x 2 , y 1 , z 1 are given in (105), (106), (108), (110). This shows that the approximation (131) of the first cumulant of the QSD is equal to the sum of the first two terms of its asymptotic approximation for large N -values, while the approximations (132) and (133) of the second and third cumulants of he QSD are equal to the first terms of the corresponding asymptotic approximations. All three approximations are therefore asymptotic for large N . We note, however, that the BGL method does not by itself allow this important conclusion. Additional information about properties of the method will be uncovered in our study below of published results number three and four, referred to as BR1 and BR2. The second published result that we comment on was presented by Matis and Kiffe (1996) in case s = 1, and by Matis et al. (1998) when s is a positive integer. A valuable contribution of these papers is that they show that a system of ordinary differential equations (ODEs) for the first few cumulants of the unconditioned random variable X (t) can be derived from the PDE for the cumulant generating function K (θ, t) of X (t). As already mentioned, we use a slight variation of this approach for deriving ODEs of the first few cumulants of the conditioned random variable X Q (t). It turns out that the first system of ODEs (for the cumulants of X (t)) serves as an approximation of the second system of ODEs (for the cumulants of X Q (t)) under the condition that R 0 > 1. The critical points of the two systems of ODEs correspond to stationary distributions of the two random variables X (t) and X Q (t), respectively. The two systems of ODEs for the cumulants are not closed, in the sense that the number of cumulants is larger than the number of equations. This may appear as an undesirable property of the problem. Closure is clearly necessary if one wants or needs exact solutions. However, for our purposes it is important to realize that closure is not needed to find approximate solutions. To achieve closure, Matis et al. (1996Matis et al. ( , 1998 assume that all cumulants of sufficiently large order are equal to zero. This assumption is clearly at odds with our finding that all cumulants of the QSD are of order O(N ). It turns out, however, that it leads to good numerical approximations of the cumulants. This behavior can be understood from our results. We note for example from Eqs. (36)-(38) that the right-hand sides of the ODEs for the first three cumulants when s = 4 depend upon all seven cumulants κ 1 − κ 7 . But we note also in this case from (93)-(98) that our assumption that the cumulants κ 4 − κ 7 are all O(N ) leads us to conclude that they have no influence on the first three terms of the asymptotic approximation of κ 1 , nor on the first two terms of the asymptotic approximation of κ 2 , or on the first term of the asymptotic approximation of κ 3 . The same conclusion can be drawn from the assumption made by Matis and Kiffe that the cumulants κ k are equal to zero for k ≥ 4. This argument shows that cumulant approximations based on cumulant closure and an assumption of cumulant neglect can lead to asymptotic approximations if they are acceptable at all, as shown numerically in several cases by Matis et al. (1998). But the assumption that some cumulant κ k equals zero is incorrect and can lead to errors if one studies, e.g., the QSD instead of its cumulants.
The remaining three published results that we discuss here are all based on the BGL method. Renshaw (2011) has used the BGL method in an effort to improve the results by Bartlett et al. (1960) by giving additional terms in the approximations of the first three cumulants of the QSD. He proceeds to formulate two approaches that we refer to as BR1 and BR2, respectively, where the letters B and R are used to refer to Bartlett and Renshaw. In each of these two approaches, he proceeds to derive three relations between cumulants that turn out to be similar to relations between critical points of the system of ODEs for the first three cumulants. To describe his findings, we quote first from Renshaw (2011) the ODEs for the first three cumulants in the case s = 1. They are derived from the PDE for the cumulant generating function and do not involve any approximations. For consistency in notation, we denote the time derivatives of the first three cumulants by A (1) , B (1) , and C (1) , respectively, when s = 1. We get The same equations are found in Matis and Kiffe (1996), and, after using the reparameterization in (10)-(13), in Nåsell (2017). The two approaches discussed by Renshaw (2011) are the third and fourth of the published results that we discuss. It turns out that the first two of the cumulant relations derived by Renshaw are in each of his two approaches equal to the expressions found by setting the cumulant derivatives κ 1 = A (1) and κ 2 = B (1) quoted above equal to zero. However, the third cumulant relation that he derives is different from what is found by setting κ 3 = C (1) equal to zero in each of the two approaches. In the approach BR1 the third cumulant relation is given by relation (3.5.21) in Renshaw's book, while it is given by the expression following (3.5.38) in the same book in the approach BR2. We use the notations C (a) = 0 and C (b) = 0 to refer to these relations. For the first one of these, we find that C (a) is written as follows after using the reparameterization in (10)- (13): We proceed to derive asymptotic approximations of the first three cumulants of the QSD for this case. As in Sect. 6, our basic assumption is that the first four cumulants have the following asymptotic behaviors: By inserting these asymptotic approximations of the first four cumulants into the expressions for A (1) , B (1) , C (a) in (137), (138), (140), we find that asymptotic approximations for them can be written as follows: where

B
(1) The basic mathematical problem at this point is to determine the six coefficients x 1 , x 2 , x 3 , y 1 , y 2 , z 1 so that the following three conditions are satisfied: These conditions are satisfied by setting the eight expressions (148)-(155) equal to zero. We proceed to solve the resulting eight equations for the six unknown coefficients x 1 , x 2 , x 3 , y 1 , y 2 , z 1 . At this point it appears that there are more equations than unknowns, but the apparent problem that this causes will readily be solved. The equation A (1) 1 = 0 is first solved for x 1 , and the spurious solution x 1 = 0 is excluded. As above, as soon as the solution is found for any of the coefficients x 1 , x 2 , x 3 , y 1 , y 2 , z 1 , it is inserted into the remaining unsolved equations. The equation B (1) 1 = 0 is then solved for y 1 , and the equation A (1) 2 = 0 is solved for x 2 . After this we find that the two expressions C We proceed to derive asymptotic approximations of the first three cumulants of the QSD for this case, using the three cumulant relations A (1) = 0, B (1) = 0, and C (b) = 0. As above, we assume that the first four cumulants have the following asymptotic behaviors: By inserting these asymptotic approximations of the first four cumulants into the expressions for A (1) , B (1) , and C (b) , we find that asymptotic approximations for them can be written as follows: Expressions for the five quantities A (1) 3 , B (1) 2 are given in (148) We set these seven expressions equal to zero, and solve the resulting seven equations for the six unknown coefficients x 1 , x 2 , x 3 , y 1 , y 2 , z 1 . The value of x 1 is found by solving the equation A (1) 2 = 0 for x 2 . The values of x 1 , x 2 , and y 1 are the same as the ones found in (159), (160), (162). Furthermore, we insert these values of x 1 , x 2 , y 1 into (173) and (174). This shows that C (b) 2 = 0 and that therefore the number of equations available for solving for the six coefficients is reduced from seven to six. We also get To determine the three coefficients x 3 , y 2 , z 1 we solve the equation C (b) 1 = 0 for z 1 , the equation B (1) 2 = 0 for y 2 , and the equation A (1) 3 = 0 for x 3 . The rather surprising result is that the values of the six coefficients x 1 , x 2 , x 3 , y 1 , y 2 , z 1 are found to be the same as in (105)-(110). This means that even though the explicit results from the BR2 method are different from the results from our preferred method, the asymptotic approximations agree.
In his development of BR2, Renshaw (2011) uses the equation C By using the reparameterization in (11) and (13) we find that this relation was shown to hold already by Bartlett et al. (1960). It is seen from (108) and (110) that this relation is asymptotic for large N . However, the arguments used by Renshaw in his derivation can be criticized. As also pointed out by Bhowmick et al. (2016), it does not make sense to assume that κ 1 is large in comparison with a 1 , b 1 , a 2 , b 2 , since these four parameters are rates whose values depend on the unit of time, which is arbitrary, while κ 1 is independent of the time unit. It is also unrealistic to assume that κ 2 is small compared with κ 1 , since it contradicts our finding that all cumulants are O(N ).
Renshaw's two efforts to extend the BGL approximation beyond what was derived by Bartlett et al. (1960) leads to two different conclusions. In one case (BR2) the results provide improvements, while they do not in the other case (BR1). This is enough to conclude that the BGL method cannot be trusted to give improved results without further investigations.
We turn now to consider the results reported by Bhowmick et al. (2016). Actually, they study a more general model than the one that we are concerned with here. The population birth rate λ n in their model is given by while their population death rate μ n equals They allow δ to be a positive number, not necessarily an integer. This leads to difficulties in determining dimensions and biological interpretations of the parameters whenever δ = 1. To avoid such difficulties, we limit ourselves to a study of their model in the special case when δ = 1. To agree with our model formulation we furthermore put β = s = 1, α 1 = 1 and α 2 = α. By using results derived by Bhowmick et al. (2016) and reported in Sect. 4 of their paper, we find that their approximation κ (B B) 1 of the first cumulant κ 1 equals where By including three terms in the asymptotic approximation of κ 1 in (179), we find that By including two terms in the asymptotic approximation of κ (B B) 2 , we find that Bhowmick et al. use four quantities A, B, C, and D in a formula that determines their approximation κ (B B) 3 of the third cumulant κ 3 as follows: These four quantities are defined in terms of the parameters that are used to describe the first of the two model formulations in Sect. 2. After reparameterization they can be expressed as follows: Thus, all four of these quantities are determined as functions of κ (B B) 1 . We use the asymptotic approximation of κ (B B) 1 in (181) to determine one-term asymptotic approximations of each of A, B, C, D. Two terms of the asymptotic approximation of κ (B B) 1 are needed for A, while one term suffices for the remaining three quantities. The results are It follows from this that the one-term asymptotic approximation of the ratio (A + B)/(C + D) equals By inserting this result and the one-term asymptotic approximation of κ (B B) 2 from (183) into (184), we get The resulting asymptotic approximations of the first three cumulants derived by Bhowmick et al. are found in (181),(183),and (194). These can be compared with the results derived by using our preferred method and found in (105)-(110). The comparisons show that the first two terms of the asymptotic approximation of κ (B B) 1 agree with the corresponding terms of the asymptotic approximation of the true first cumulant κ 1 , but there is disagreement between the third terms. Moreover, the first term of the asymptotic approximation of κ (B B) 2 agrees with the first term of the asymptotic approximation of the true second cumulant κ 2 , but there is disagreement between the second terms. Finally, there is disagreement between the one-term asymptotic approximation of κ (B B) 3 and the corresponding one-term asymptotic approximation of the true third cumulant κ 3 . Independent supports for these conclusions are found from the numerical evaluations of the error terms of the results presented by Bhowmick et al. and given in Table 5. We see here that the error term for the BB approximation of κ 1 is reduced by the factor two for each doubling of N . This is indicative of an error that is of magnitude of O(1/N ). In addition, the error term for the BB approximation of κ 2 is found to be approximately independent of N , which indicates that it is of the magnitude of O(1). The error term for the BB approximation of κ 3 is, finally, found to be approximately proportional to N , which makes the corresponding approximation useless. The reason for the larger errors in the results arrived at by Bhowmick et al. is that one or more of the approximation steps taken in the application of BGL method brings in errors. This is similar to the reason for the large errors of the BR1 result, commented on above.
It is interesting to note that the study by Bhowmick et al. (2016) is the only one among the five studies commented on in this section that uses the parameterization that accompanies the second of the two model formulations described in Sect. 2. However, we note also that Bhowmick et al. do not make use of the important property that this model formulation offers, namely the availability of the large parameter N that can be used for deriving asymptotic approximations.

Concluding Comments
The method that we have used here to derive asymptotic approximations of the first few cumulants of the QSD of the stochastic power law logistic model emphasizes the importance of the second model formulation in Sect. 2. It gives access to two parameters that are basic for our study, namely the maximum population size N and the threshold parameter R 0 . We have shown that the condition R 0 > 1 for large N is required for the results that we have presented here, while both the moment closure method and the three approaches BR1, BR2, and BB based on the BGL method have the serious weaknesses that they cannot produce conditions for validity of the approximations that are derived. In addition, we note that magnitudes of approximation errors are easy to establish in our method, as they are for any asymptotic approximation, while they can not be produced in the moment closure method, nor in approaches based on the BGL method. We note furthermore that spurious solutions that could require large efforts to eliminate appeared in early studies based on moment closure, while the spurious solutions that appear in the method that we use here are easy to identify and eliminate. Our results show that the number of spurious solutions is equal to the parameter s whenever s is a positive integer. As a further comment, we note that the dependence of the new approximations on the model parameters is explicit in our approach, while they are unknown in results based on moment closure and in BR1 and BR2.
We have found that our method based on the second model formulation of Sect. 2, and followed by a search for asymptotic approximations, provides a powerful approach for determination of approximations of the first few cumulants of the QSD for the power law logistic model. We have also shown that the method of determining asymptotic approximations can be used to study other approaches to the same problem.
We conclude that the method of determining ODEs for the first few cumulants of the QSD introduced by Matis and Kiffe (1996) is preferred over the BGL method for deriving relations between the cumulants. An important reason for this is that the ODEs for the first few cumulants of the QSD are exact, while the BGL method can cause errors of unknown magnitude. We emphasize the obvious fact that whenever exact solutions of a mathematical problem are difficult to establish, then one should search for approximations. Furthermore, approximation methods are then definitely preferred that give both conditions for validity of the approximations and magnitudes of approximation errors. Because of this, we conclude that our method for determining asymptotic approximations of the first few cumulants is preferred over any other method that has been used on this problem.