Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Mathematical Models for Population Growth with Variable Carrying Capacity: Analytical Solutions

Mathematical Models for Population Growth with Variable Carrying Capacity: Analytical Solutions Article Mathematical Models for Population Growth with Variable Carrying Capacity: Analytical Solutions 1, 2 M. Rodrigo * and D. Zulkarnaen School of Mathematics and Applied Statistics, University of Wollongong, Wollongong, NSW 2522, Australia Department of Mathematics, Universitas Islam Negeri Sunan Gunung Djati, Bandung 40614, Indonesia * Correspondence: marianit@uow.edu.au Abstract: A general population model with variable carrying capacity consisting of a coupled system of nonlinear ordinary differential equations is proposed, and a procedure for obtaining analytical solutions for three broad classes of models is provided. A particular case is when the popula- tion and carrying capacity per capita growth rates are proportional. As an example, a generalised Thornley–France model is given. Further examples are given when the growth rates are not propor- tional. A criterion when inflexion may occur is also provided, and results of numerical simulations are presented. Keywords: population model; variable carrying capacity; exact solution; power-law logistic; logistic; Gompertz MSC: 92D25; 34A12; 34A34 1. Introduction The logistic model has been extensively used to study the cause and effect relationship Citation: Rodrigo, M.; Zulkarnaen, D. between the ‘carrying capacity’ (i.e., the population size that available resources can support) Mathematical Models for Population and the population size. See, for instance, the references in Brauer & Castillo-Chávez [1], Growth with Variable Carrying Capacity: Analytical Solutions. Gotelli [2], Pastor [3], as well the seminal papers by Verhulst [4] and Pearl & Reed [5]. AppliedMath 2022, 2, 466–479. It is typically assumed that the carrying capacity does not vary with time and, therefore, https://doi.org/10.3390/ the logistic model exhibits a sigmoidal shape when the population size is plotted as a appliedmath2030027 function of time. However, many phenomena such as human population growth exhibit more complex behaviours, unlike population species grown in laboratory cultures, for Academic Editor: Leonid Shaikhet example. Thus, it is of theoretical and practical interest to investigate mathematical models Received: 6 July 2022 that incorporate variable carrying capacities. One of the main objectives of studying such Accepted: 18 August 2022 models is to explore how different functional forms of the carrying capacity influence the Published: 29 August 2022 dynamics of the population variable and its long-term behaviour. Meyer [6] and Meyer & Ausubel [7] studied a bi-logistic model derived from a logistic Publisher’s Note: MDPI stays neutral model but with a sigmoidal time-dependent carrying capacity. Cohen [8] proposed a with regard to jurisdictional claims in published maps and institutional affil- human population growth model with a variable carrying capacity that is a function of the iations. population size itself. The aforementioned models make the case that the inclusion of a variable carrying capacity is more reflective of the human condition. Safuan et al. [9] proposed a coupled system of ordinary differential equations (ODEs) to describe the interaction between the population and its carrying capacity. The model they Copyright: © 2022 by the authors. considered does not require prior knowledge of the carrying capacity, nor does it impose Licensee MDPI, Basel, Switzerland. constraints on the initial conditions. Assuming a special form of the carrying capacity in This article is an open access article the logistic equation, the same authors obtained an analytical solution in series form [10]. distributed under the terms and As pointed out by Cohen [8], there is no consensus with regards to appropriate models conditions of the Creative Commons for human carrying capacity. However, most would accept that human carrying capacity Attribution (CC BY) license (https:// is influenced by food availability, amongst other factors. Hopfenberg [11] postulated that creativecommons.org/licenses/by/ food production data are the sole variable that influences human carrying capacity, and 4.0/). AppliedMath 2022, 2, 466–479. https://doi.org/10.3390/appliedmath2030027 https://www.mdpi.com/journal/appliedmath AppliedMath 2022, 2 467 a simple linear relationship between human carrying capacity and the food production index is assumed. More recently, Zulkarnaen & Rodrigo [12] proposed three classes of human population dynamical models of logistic type where the carrying capacity is a function of the food production index. They employed an integration-based parameter estimation technique [13] to derive explicit formulas for the model parameters. Using actual population and food production index data, the results of numerical simulations of their models suggested that an increase in food availability implies an increase in carrying capacity, but the carrying capacity is ‘self-limiting’ and does not increase indefinitely. Thornley & France [14] proposed an ‘open-ended’ form of the logistic equation by considering a system of two ODEs representing the coupled processes of growth and development. Their model is ‘open-ended’ in the sense that dynamic changes in nutrition and environment can influence growth and development, which, in turn, may affect the asymptotic carrying capacity value. Subsequently, Thornley et al. [15] found an analytical solution to the Thornley–France model in the case of constant parameters. The solution of the system of ODEs is expressed in terms of the solution of a single ODE of power-law logistic type, also referred to as the q-logistic model, which frequently arises in ecology and elsewhere [16]. A related article by Wu et al. [17] formulated the variable carrying capacity by exploring a resource dynamic-based feedback mechanism underlying the population growth models. The inclusion of variable carrying capacities in interacting species such as in predator–prey models have been considered in Al-Moqbali et al. [18], Al-Salti et al. [19]. Power-law logistic models have been investigated by von Bertalanffy [20] and Richards [21]. The principal (nonnegative) parameter of these models is denoted by q. The Gompertz model and the logistic model are recovered when q = 0 and q = 1, respectively. Larger values of q behave like a logistic model but with an increasingly sharper cessation of growth as the asymptote (i.e., the constant carrying capacity limit) is approached [15]. As a fraction of the asymptote, inflexion can occur over the range from (Gompertz) through (‘ordinary’ logistic) and then to 1 (for large q). Determining the point where inflexion takes place can be especially important when fitting the model to actual data that exhibit a sigmoidal trend. See Banks [22] for a detailed analysis of the q-logistic model. More recently, Albano et al. [23] considered a general growth curve that includes the Malthus, Richards, Gompertz and other models. Their generalised model is essentially a Bernoulli ODE, which is well known to be analytically tractable. They investigated the analytical and numerical properties of the solution as the parameters varied. Here we propose a general population model with a variable carrying capacity, which includes the Thornley–France [14], Safuan–Jovanoski–Towers–Sidhu [9] and Meyer– Ausubel [7] models as special cases. Moreover, when the carrying capacity is kept constant, the proposed model system reduces to a single ODE population model and recovers the Gompertz, ‘ordinary’ logistic and q-logistic models, amongst others. The idea is to extract the essential properties of such models without getting ‘bogged down’ by particular cases. We provide a procedure for obtaining, when possible, the analytical solution of this general population model. Different models can then be chosen depending on the particular phe- nomenon being studied. An important tractable special case is when the per capita growth rates of the population and carrying capacity are proportional to each other. We give a criterion for when inflexion may occur. Several illustrative examples are also provided. 2. Population Models with Variable Carrying Capacities Consider the initial value problem (IVP) dN dK = N f (N, K), = K g(N, K), N(0) = N , K(0) = K , (1) 0 0 dt dt where N(t) and K(t) are the population and carrying capacity, respectively, at time t. The initial values N and K are positive and given. Let f and g be sufficiently smooth 0 0 bivariate functions on R  R , where R = (0,¥). Denote by D , where j = 1, 2, the + + + partial derivative with respect to the independent variable in the jth position. We assume AppliedMath 2022, 2 468 1 dN that the population per capita growth rate decreases with increasing population N dt (D f (x, y) < 0) and increases with increasing carrying capacity (D f (x, y) > 0). As 1 2 for the signs of D g(x, y) and D g(x, y), it is not obvious what these should be since 1 2 1 dK the behaviour of the carrying capacity per capita growth rate may depend on the K dt particular population species. When g is identically zero, then K(t) = K and (1) reduces to dN = N f (N, K ), N(0) = N . 0 0 dt A well-known example is the q-logistic model [22] r N f (N, K ) = 1 , q K where r > 0 is the intrinsic growth rate, and q > 0 is a parameter related to the point of inflexion of the solution. The case q = 1 gives the ‘ordinary’ logistic model. In the limit as q ! 0 , the q-logistic model reduces to the Gompertz model [24] f (N, K ) = r log . In the case of a variable carrying capacity, the Thornley–France model [14] takes the form N N (2) f (N, K) = a 1 , g(N, K) = b 1 , K K where a, b > 0. Note that D g(x, y) > 0 and D g(x, y) < 0. Safuan et al. [9] proposed 1 2 the model (3) f (N, K) = a 1 , g(N, K) = b cN, where a, b, c > 0. Here we see that D g(x, y) < 0 and D g(x, y) = 0. Meyer [6] and Meyer & Ausubel [7] assumed that f (N, K) = a 1 , g(N, K) = b cK, where a, b, c > 0. This time, D g(x, y) = 0 and D g(x, y) < 0. 1 2 2.1. Construction of Analytically Tractable Population Models with Variable Carrying Capacities Suppose that there exists a positive C -function F such that y = F(x) solves the IVP yg(x, y) dx x f (x, y) dy = 0, y = K when x = N . (4) 0 0 Define G(x) = dz. (5) z f (z, F(z)) We claim that G(N(t)) = t, K(t) = F(N(t)) (6) gives the formal solution of (1) in implicit form. Indeed, (4) and (6) imply that F(x)g(x, F(x)) K(t)g(N(t), K(t)) 0 0 F (x) = , F (N(t)) = . x f (x, F(x)) N(t) f (N(t), K(t)) AppliedMath 2022, 2 469 Implicitly differentiating G(N(t)) = t with respect to t yields 0 0 0 1 = G (N(t))N (t) = N (t) N(t) f (N(t), K(t)) dN or = N f (N, K). Hence dt K(t)g(N(t), K(t)) 0 0 0 0 K (t) = F (N(t))N (t) = N (t) = K(t)g(N(t), K(t)) N(t) f (N(t), K(t)) dK or = K g(N, K). It is clear that G(N(0)) = G(N ) = 0 and K(0) = F(N(0)) = F(N ) = 0 0 dt K . This proves the claim. Therefore, the task of finding an analytical solution of (1) essentially reduces to solving the first-order ODE (4). 2.2. Determination of Inflexion Points for the Population Species Variable Before we consider some examples, let us first investigate where inflexion may occur for the population species variable. As mentioned previously, this is particularly relevant during model fitting. Using (6) and differentiating the first ODE in (1) with respect to t, we have d N dN = [ f (N, F(N)) + N D f (N, F(N)) + N D f (N, F(N))F (N)] . 1 2 dt dt Therefore, inflexion for the function N may occur at some positive root N of the equation H(z) = f (z, F(z)) + zD f (z, F(z)) + zD f (z, F(z))F (z) = 0. (7) 1 2 If this occurs, it will be when t = t , where t = G(N(t )) = G(N ) = dz. (8) z f (z, F(z)) 2.3. First Class of Analytically Tractable Models: f and g Are Proportional Suppose that there exists a 2 Rnf0g such that g(x, y) = a f (x, y). This basically as- 1 dN 1 dK sumes that the per capita growth rates and are proportional. Then (4) and (5) give N dt K dt a a y = F(x) = K N x , G(x) = dz, N z f (z, K N z ) 0 0 respectively. From (6) we obtain the exact solution N(t) a a dz = t, K(t) = K N [N(t)] (9) a 0 N z f (z, K N z ) 0 0 of the system (1). Example 1. Let q q a x b x f (x, y) = 1 , g(x, y) = 1 , (10) q y q y where a, b > 0 and q  0. Then g(x, y) = a f (x, y), where a = < 0. Thus, (1) becomes q q dN a N dK b N = N 1 , = K 1 , N(0) = N , K(0) = K . (11) 0 0 dt q K dt q K AppliedMath 2022, 2 470 The ‘value’ when q = 0 is meant to be understood as the limit when q ! 0 ; this is related to the Gompertz model in the case of a constant carrying capacity. The Thornley–France model (2) is a special case of (10) if we take q = 1. From the first equation of (9), we obtain Z Z a q q(1a) N(t) (K N ) [N(t)] q 1 1 0 1 t = dz = du. a q(1a) q q(1a) q a a(1 a) u(1 u) N z[1 (K N ) z ] (K N ) N 0 0 0 0 0 (12) Recall the formula 1 u I(u) = du = log , 0 < u < 1. u(1 u) 1 u If I(b ) I(a ) = c , where 0 < a < b < 1 and c > 0, then it is not difficult to show that 0 0 0 0 0 0 b = . 1a 0 c 1 + e q(1a) a q 1 q a q q(1a) Taking a = (K N ) N = (K N ) , b = (K N ) [N(t)] and 0 0 0 0 0 0 0 0 0 c = a(1 a)t, we deduce that the exact solution of (11) is " # q(1a) a q (K N ) 0 a a N(t) = , K(t) = K N [N(t)] . (13) 1 0 1(K N ) 0 a(1a)t 1 + e (K N ) As a, b > 0 and a = < 0, we have from (13) that 1 1 a a 1a 1a N = lim N(t) = (K N ) , K = lim K(t) = (K N ) . ¥ 0 ¥ 0 0 0 t!¥ t!¥ This is consistent with the fact that f (N , K ) = g(N , K ) = 0 in (11) yields that any ¥ ¥ ¥ ¥ equilibrium of (11), and there are infinitely many, lies on the curve K = N . ¥ ¥ Let us now look at the inflexion points of N by studying the roots of H(z) = 0 in (7). Straightforward calculations give q q a x a x D f (x, y) = , D f (x, y) = , xF (x) = aF(x), x y y y xD f (x, y) = a + q f (x, y), yD f (x, y) = a q f (x, y). 1 2 If we define H (x, y) = f (x, y) + xD f (x, y) + ayD f (x, y), 0 1 2 then H(z) = 0 if and only if H (z, F(z)) = 0. However, H (x, y) = [1 + q(1 a)] f (x, y) 0 0 a(1 a), so that H (z, F(z)) = 0 if and only if a z a(1 a) 1 = f (z, F(z)) = . q F(z) 1 + q(1 a) This is equivalent to q 1a 1 z z z 1 1a = = = , K = (K N ) . ¥ 0 1 + q(1 a) F(z) K K N z ¥ Thus H(N ) = 0 if and only if q(1a) N = K . (14) 1 + q(1 a) AppliedMath 2022, 2 471 Finally, we see from (8) that q(1a) q(1a) 1 N N t = I I . (15) a(1 a) K K ¥ ¥ As I is an increasing function of u, and assuming that N < N , it follows that t > 0. It was already noted that q = 1 in (10) reduces to the Thornley–France model (2). Then, if q = 1, (13) recovers the analytical solution found in Thornley et al. [15]. Equation (14) indicates where inflexion occurs as a fraction of the asymptotic carrying capacity, while (15) gives the time of inflexion (compare with (17) and (18), respectively, found in Thornley et al. [15] with an appropriate renaming of parameters). On the other hand, when we let q ! ¥, then (14) shows that tends to unity so that, similar to the Thornley–France model when q = 1, exponential growth is sustained for longer and the inflexion value N moves closer to the asymptotic carrying capacity value K . If q > 0 and a = 0 (corresponding to b = 0 in (10)), then K(t) = K and (1) reduces to a single ODE dN a N = N 1 , dt q K which is the q-logistic model. In particular, if q = 1 and a = 0 (i.e., the ‘ordinary’ logistic model), then (14) shows that the inflexion value is one-half of the asymptotic carrying capacity value, which + N 1 b is well known. If q ! 0 , then we deduce from (14) that tends to for any a = < 0. K e a This is similar in behaviour to the case when a = 0 (i.e., b = 0 in (10)), so that K(t) = K and (1) recovers the Gompertz model dN K = aN log . dt N N 1 + Thus tends to as q ! 0 even for the variable carrying capacity system (11). K e Here we include some results of numerical simulations of (11) (equivalently, (13)). Let a = 0.2, b = 0.1, N = 5 and K = 20. Take four different values of q = 0.01, 1, 5, 20 to see and compare the 0 0 population dynamics, as shown in Figure 1, as well as the carrying capacity behaviour in Figure 2. As can be seen in both figures, the larger the value of q, the longer it takes for the population and carrying capacity to reach equilibrium. 20.000 13.333 6.667 inflexion point 0.000 0.000 50.000 100.000 150.000 Figure 1. Profiles of N(t) in (13) for q = 0.01, 1, 5, 20 (with a = 0.2, b = 0.1, N = 5 and K = 20) and 0 0 inflexion points calculated from (14) and (15). AppliedMath 2022, 2 472 20.000 13.333 6.667 0.000 0.000 50.000 100.000 150.000 Figure 2. Profiles of K(t) in (13) for q = 0.01, 1, 5, 20 (with a = 0.2, b = 0.1, N = 5 and K = 20). 0 0 Figure 3 (q = 0.01) and Figure 4 (q = 1.2) depict some trajectories in the NK-plane and illustrate that the carrying capacity experiences a decline as the population grows larger. This is a consequence of the proportionality of the per capita growth rates with a negative proportionality constant. Moreover, when the initial conditions are varied, the trajectories approach the equilibrium line K = N . This behaviour is typical for other values of q as well. ¥ ¥ 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 Figure 3. Representative solution trajectories of (11) in the NK-plane with q = 0.01, a = 0.2 and b = 0.1. AppliedMath 2022, 2 473 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 Figure 4. Representative solution trajectories of (11) in the NK-plane with q = 1.2, a = 0.2 and b = 0.1. 2.4. Second Class of Analytically Tractable Models: f and g Are Homogeneous Functions Suppose that f and g in (1) are homogeneous functions, i.e., f (lx, ly) = l f (x, y), g(lx, ly) = lg(x, y) for every l > 0. Let v = in (4), so that g(x, y) = g(x, xv) = xg(1, v) and f (x, y) = f (x, xv) = x f (1, v). Then (4) is transformed into the separable ODE dv g(1, v) f (1, v) x = v , dx f (1, v) whose general solution is y f (1, v) F(v) = F = log(x) + log(C), F(v) = dv, (16) x v[g(1, v) f (1, v)] where C > 0 is an arbitrary constant of integration. Hence, assuming that F is invertible, we deduce that the solution of (4) is 1 0 y = F(x) = xF (log(Cx)), log(C) = F log(N ). (17) Next, we determine G. From (5), we see that G(x) = dz. (18) z f (z, zF (log(Cz))) Finally, (6) implies that the exact solution of (1) in the case when both f and g are homoge- neous functions is N(t) dz = t, K(t) = N(t)F (log(C N(t))). (19) z f (z, zF (log(Cz))) Example 2. Let x cx f (x, y) = a 1 , g(x, y) = a , y y AppliedMath 2022, 2 474 where 0 < a < c, so that (1) becomes dN N dK cN = aN 1 , = K a , N(0) = N , K(0) = K ; (20) 0 0 dt K dt K compare with model (3) proposed by Safuan et al. [10]. Observe that f and g are homogeneous functions but are not proportional since a < c. Then (16) gives h  i a y a y y F(v) = [v log(v)], F = log = log(Cx) (21) a c x a c x x and, therefore, c K K K 0 0 0 a 0 1 log(C N ) = log or (C N ) = e . (22) 0 0 a N N N 0 0 0 The second equation in (21) can be rewritten as y a c c x a e = exp log(Cx) = (Cx) . (23) x a Let us recall some properties of the Lambert W-function. Suppose that u, z 2 R. Then z = W(u) 1 1 satisfies the equation ze = u if u  . Furthermore, W(u) < 0 for  u < 0 and e e W(u) 1 W (u) = , u 6= , 0. (24) u[1 + W(u)] e Identifying z = , u = (Cx) , 1 1 ac and assuming that 0 < x  e so that u  , we obtain from (23) that C e y = F(x) = xW((Cx) ). We see that F(x) > 0 since the W-function is negative here because  u < 0. Following a similar argument, we deduce from (21) that ac c 1 z 1 1 a a F (z) = W(e ), F (log(Cz)) = W((Cz) ). Therefore, (18) yields Z Z 1 1 a x (Cx) 1 W((Cz) ) 1 W(u) G(x) = dz = du 1 1 a a c a a u[1 + W(u)] N (C N ) 0 z[1 + W((Cz) )] 0 1 c c 1 1 a a = [W((Cx) ) W((C N ) )], c a where we used (24) in the last step. Equation (22) gives c K K 0 0 W((C N ) ) = W e = N N 0 0 z z since z = W(u) = W(ze ) if u = ze . Hence the exact solution of (20) from (19) is c K c 1 1 a a W([C N(t)] ) = + (c a)t, K(t) = N(t)W([C N(t)] ). (25) Example 3. Suppose that y y f (x, y) = a log , g(x, y) = a log + b, x x AppliedMath 2022, 2 475 where a, b > 0. Therefore, (1) becomes dN K dK K = aN log , = K a log + b , N(0) = N , K(0) = K . (26) 0 0 dt N dt N We see that f and g are homogeneous functions but are not proportional as b > 0. Then (16) implies that h  i a y a y F(v) = [log(v)] , F = log = log(Cx) 2b x 2b x and 2 1 a K 1 a K 0 0 log(C) = log log(N ), [log(C N )] = log . (27) 0 0 2b N 2b N 0 0 Since 1 1 h i 2b 2 2b 2 1 1 F (v) = exp v , F (log(Cv)) = exp log(Cv) , a a Equation (17) yields h i 2b y = F(x) = x exp log(Cx) . Furthermore, h i 2b 1 1 2 2 f (z, zF (log(Cz))) = a log(Cz) = (2ab) [log(Cz)] and from (18), we obtain Z Z x log(Cx) 1 1 1 G(x) = dz = u du 1 1 1 N log(C N ) 2 2 2 (2ab) 0 z[log(Cz)] (2ab) 0 1 1 2 1 2 1 2 2 2 2 = [log(Cx)] [log(C N )] . ab ab 1 1 1 ab 2 2 2 Substituting into (19) gives [log(C N(t))] = [log(C N )] + ( ) t, which simplifies to K ab log(C N(t)) = log(C N ) + at log + t N 2 using (27). Thus, at ab 2 N(t) = N e . Moreover, (19) also implies that h i 2b 2 K(t) = N(t) exp log(C N(t)) , so that h i 2b log(K(t)) = log(N(t)) + log(C N(t)) 1 1 1 2b 1 2b ab 2 2 2 = log(N(t)) + [log(C N )] + t a a 2 = log(N(t)) + log + bt. 0 AppliedMath 2022, 2 476 Hence at at K K K ab 2 K ab 2 0 0 0 0 bt bt t t +bt 2 2 K(t) = e N(t) = e N e = K e . 0 0 N N N N 0 0 0 0 Thus, the exact solution of (26) is at at K ab 2 K ab 2 0 0 t t +bt 2 2 N(t) = N e , K(t) = K e . (28) 0 0 N N 0 0 2.5. Third Class of Analytically Tractable Models: f and g Determine an Exact ODE Assume that f and g are such that ¶ ¶ [x f (x, y)] + [yg(x, y)] = f (x, y) + xD f (x, y) + g(x, y) + yD g(x, y) = 0. (29) 1 2 ¶x ¶y This implies that (4) is an exact ODE whose general solution is Y(x, y) = C = Y(N , K ), 0 0 where C is an arbitrary constant of integration and D Y(x, y) = yg(x, y), D Y(x, y) = x f (x, y). (30) 1 2 Then Z Z Y(x, y) = yg(x, y) dx + Q(y) = x f (x, y) dy + P(x), (31) where P = P(x) and Q = Q(y) are arbitrary functions of integration. Differentiating Y in (31) with respect to y and using the second equation in (30) and then (29), we get Z Z Z Q(y) = x f (x, y) f (x, y) dx xD f (x, y) dx dy. (32) An analogous calculation yields Z Z Z P(x) = yg(x, y) f (x, y) dy xD f (x, y) dy dx. (33) Assuming that Y(x, y) = C = Y(N , K ) can be solved for y = F(x), then we have 0 0 from (5) and the second equation in (30) that Z Z x x 1 1 G(x) = dz = dz. (34) z f (z, F(z)) D Y(z, F(z)) N N 2 0 0 Therefore, (6) expresses the exact solution of (1) as N(t) dz = t, K(t) = F(N(t)). (35) D Y(z, F(z)) N 2 Looking at condition (29) more closely, we conclude that Z Z yg(x, y) = f (x, y) dy xD f (x, y) dy + R(x) (36) 1 AppliedMath 2022, 2 477 for some arbitrary function R = R(x). Hence the exactness condition (29) necessarily implies the above form for g. Thus from (31) and (36), we have Z Z Z Y(x, y) = f (x, y) dy + x D f (x, y) dy R(x) dx Z Z Z + x f (x, y) f (x, y) dx xD f (x, y) dx dy Z Z = x f (x, y) dy R(x) dx. This simplified form for Y is then substituted into the first equation in (35). Example 4. Let a > 0. If 1 1 f (x, y) = a 1 = a axy , D f (x, y) = ay , then (36) evaluates to yg(x, y) = ay + 2ax log(y), taking R(x) = 0 for simplicity. Therefore, (1) becomes dN N dK = aN 1 , = aK + 2aN log(K), N(0) = N , K(0) = K . (37) 0 0 dt K dt We have Y(x, y) = x(a axy ) dy = ax[y x log(y)], where C = Y(N , K ) = aN [K N log(K )]. The equation Y(x, y) = C gives 0 0 0 0 0 0 y C y C 2 1 2 log(y) = x , e = x exp x . x a x a Hence, as in Example 2, C C 1 1 2 1 2 y = F(x) = xW x exp x if x exp x  . a a e Then h i 2 1 D Y(x, y) = ax ax y , D Y(z, F(z)) = az 1 2 2 F(z) and 1 C 2 1 + W(z exp( z ) 1 = . 1 2 F(z) W(z exp( z ) From (34) we obtain 1 2 W(z exp( z ) G(x) = dz, 1 2 z[1 + W(z exp( z )] which cannot be further evaluated analytically (unlike in Example 2). The exact solution of (37) using (35) is 1 C 2 N(t) W(z exp( z ) dz = at, 1 2 z[1 + W(z exp( z )] (38) 1 2 K(t) = N(t)W [N(t)] exp [N(t)] . a AppliedMath 2022, 2 478 3. Concluding Remarks In this article, we proposed population models with variable carrying capacities modelled by a coupled system (1) of two nonlinear ODEs and found their analytical solutions. While it was clear that the assumptions D f (x, y) < 0 and D f (x, y) > 0 are 1 2 reasonable since they describe the behaviour of the population per capita growth rate, we showed through several explicit examples that corresponding assumptions for D g(x, y) and D g(x, y) that describe the carrying capacity per capita growth rate are not obvious and may be model dependent. One possible reason for this is that carrying capacity is not directly observable, unlike the population size. If the per capita growth rates are proportional, as in Example 1, then in addition to the analytical solution, we also found a criterion for the occurrence of inflexion in the population profile as a fraction of the asymptotic carrying capacity. This criterion does not apply to the models in Examples 2 and 3 since they do not have a nontrivial equilibrium. However, it should apply for the model in Example 4 with a nontrivial 1 1 2 2 equilibrium (N , K ) = (e , e ). Further classes of analytically tractable models of ¥ ¥ the form (1), not necessarily in the context of population growth with variable carrying capacity, can, of course, be found by considering analytically tractable cases of (4) following the idea in Rodrigo [25]. Future work would involve a more careful investigation of the assumptions for the carrying capacity per capita growth rate. While this article focused on the analytical aspects of population models with variable carrying capacities, fitting the models to actual population data is the next important step. The estimation of parameters in the models can be undertaken by adapting the arguments in Zulkarnaen & Rodrigo [12], Holder & Rodrigo [13]. Parameter estimation using an integration-based technique in population models with variable carrying capacities that depend on food availability is the subject of an article by the authors that is currently under review [26]. Another research direction is the modelling of interacting population species where the carrying capacities are not fixed anymore but may depend on time and/or space. Results for predator–prey models have been obtained in Al-Moqbali et al. [18], Albano et al. [23], although spatio-temporal models (e.g., chemotaxis [27]) can also be considered. Travelling waves are an important class of biologically relevant solutions for spatio-temporal models described by partial differential equations. A travelling wave coordinate transformation leads to a system of higher dimensional ODEs, for which analytically tractable models can potentially be identified using the techniques in the current article. Author Contributions: Conceptualization, M.R.; Formal analysis, M.R.; Investigation, M.R.; Method- ology, M.R.; Project administration, M.R.; Software, M.R. and D.Z.; Supervision, M.R.; Validation, M.R.; Visualization, M.R. and D.Z.; Writing—original draft, M.R. and D.Z.; Writing—review & editing, M.R. and D.Z. All authors have read and agreed to the published version of the manuscript. Funding: This research received no external funding. Data Availability Statement: Not applicable. Conflicts of Interest: The authors declare no conflict of interest. References 1. Brauer, F.; Castillo-Chávez, C. Mathematical Models in Population Biology and Epidemiology, 2nd ed.; Springer: New York, NY, USA, 2012. 2. Gotelli, N. A Primer of Ecology, 2nd ed.; Sinauer Associates: Sunderland, MA, USA, 1998. 3. Pastor, J. Mathematical Ecology of Populations and Ecosystems; Wiley-Blackwell: London, UK, 2008. 4. Verhulst, P.F. Notice sur la loi que la population suit dans son accroissement. Corresp. Mathématique Phys. 1838, 10, 113–121. 5. Pearl, R.; Reed, L.J. A further note on the mathematical theory of population growth. Proc. Natl. Acad. Sci. USA 1922, 8, 365–368. [CrossRef] 6. Meyer, P.S. Bi-logistic growth. Technol. Forecast. Soc. Chang. 1994, 47, 89–102. [CrossRef] 7. Meyer, P.S.; Ausubel, J.H. Carrying capacity: A model with logistically varying limits. Technol. Forecast. Soc. Chang. 1999, 61, 209–214. [CrossRef] AppliedMath 2022, 2 479 8. Cohen, J.E. Population growth and the Earth’s human carrying capacity. Science 1995, 269, 341–346. [CrossRef] 9. Safuan, H.M.; Jovanoski, Z.; Towers, I.N.; Sidhu, H.S. Coupled logistic carrying capacity. ANZIAM J. 2012, 53, 172–184. [CrossRef] 10. Safuan, H.M.; Jovanoski, Z.; Towers, I.N.; Sidhu, H.S. Exact solution of a non-autonomous logistic population model. Ecol. Model. 2013, 251, 99–102. [CrossRef] 11. Hopfenberg, R. Human carrying capacity is determined by food availability. Popul. Environ. 2003, 25, 109–117. [CrossRef] 12. Zulkarnaen, D.; Rodrigo, M.R. Modelling human carrying capacity as a function of food availability. ANZIAM J. 2020, 62, 318–333. [CrossRef] 13. Holder, A.B.; Rodrigo, M.R. An integration-based method for estimating parameters in a system of differential equations. Appl. Math. Comput. 2013, 219, 9700–9708. [CrossRef] 14. Thornley, J.H.M.; France, J. An open-ended logistic-based growth function. Ecol. Model. 2005, 184, 257–261. [CrossRef] 15. Thornley, J.H.M.; Shepherd, J.J.; France, J. An open-ended logistic-based growth function: analytical solutions and the power-law logistic model. Ecol. Model. 2007, 204, 531–534. [CrossRef] 16. Gilpin, M.E.; Ayala, F.J. Global models of growth and competition. Proc. Natl. Acad. Sci. USA 1973, 70, 3590–3593. [CrossRef] 17. Wu, H.; Chakraborty, A.; Li, B.-L.; Kenerley, C.M. Formulating variable carrying capacity by exploring a resource dynamics-based feedback mechanism underlying the population growth models. Ecol. Complex. 2009, 6, 403–412. [CrossRef] 18. Al-Moqbali, M.K.A.; Al-Salti, N.S.; Elmojtaba, I.M. Prey-predator models with variable carrying capacity. Mathematics 2018, 6, 102. [CrossRef] 19. Al-Salti, N.; Al-Musalhi, F.; Gandhi, V.; Al-Moqbali, M.; Elmojtaba, I. Dynamical analysis of a prey-predator model incorporating a prey refuge with variable carrying capacity. Ecol. Complex. 2021, 45, 100888. [CrossRef] 20. von Bertalanffy, L. Quantitative laws for metabolism and growth. Q. Rev. Biol. 1957, 32, 217–231. [CrossRef] 21. Richards, F.J. A flexible growth function for empirical use. J. Exp. Bot. 1959, 10, 290–300. [CrossRef] 22. Banks, R.B. Growth and Diffusion Phenomena: Mathematical Frameworks and Applications; Springer: Berlin/Heidelberg, Germany, 1994. 23. Albano, G.; Giorno, V.; Román-Román, P.; Torres-Ruiz, F. Study of a general growth model. Commun. Nonlinear Sci. Numer. Simul. 2022, 107, 106100. [CrossRef] 24. Gompertz, B. On the nature of the function expressing the law of human mortality and on a new mode of determining the value of life contingencies. Philos. Trans. R. Soc. Lond. 1825, 115, 513–585. 25. Rodrigo, M.R. An elementary method for obtaining the general solution of a system of ordinary differential equations. Electron. J. Differ. Equ. 2021, 2021, 1–20. 26. Zulkarnaen, D.; Rodrigo, M.R. Population models with variable carrying capacities that depend on food availability. 2022, under review. 27. Hillen, T.; Painter, K.J. A user ’s guide to PDE models for chemotaxis. J. Math. Biol. 2009, 58, 183–217. [CrossRef] [PubMed] http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png AppliedMath Multidisciplinary Digital Publishing Institute

Mathematical Models for Population Growth with Variable Carrying Capacity: Analytical Solutions

Rodrigo, M.; Zulkarnaen, D.
AppliedMath , Volume 2 (3) – Aug 29, 2022
Download PDF
14 pages

Loading next page...
 
/lp/multidisciplinary-digital-publishing-institute/mathematical-models-for-population-growth-with-variable-carrying-vxc358Mq4H

References (19)

Publisher
Multidisciplinary Digital Publishing Institute
Copyright
© 1996-2022 MDPI (Basel, Switzerland) unless otherwise stated Disclaimer The statements, opinions and data contained in the journals are solely those of the individual authors and contributors and not of the publisher and the editor(s). MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Terms and Conditions Privacy Policy
ISSN
2673-9909
DOI
10.3390/appliedmath2030027
Publisher site
See Article on Publisher Site

Abstract

Article Mathematical Models for Population Growth with Variable Carrying Capacity: Analytical Solutions 1, 2 M. Rodrigo * and D. Zulkarnaen School of Mathematics and Applied Statistics, University of Wollongong, Wollongong, NSW 2522, Australia Department of Mathematics, Universitas Islam Negeri Sunan Gunung Djati, Bandung 40614, Indonesia * Correspondence: marianit@uow.edu.au Abstract: A general population model with variable carrying capacity consisting of a coupled system of nonlinear ordinary differential equations is proposed, and a procedure for obtaining analytical solutions for three broad classes of models is provided. A particular case is when the popula- tion and carrying capacity per capita growth rates are proportional. As an example, a generalised Thornley–France model is given. Further examples are given when the growth rates are not propor- tional. A criterion when inflexion may occur is also provided, and results of numerical simulations are presented. Keywords: population model; variable carrying capacity; exact solution; power-law logistic; logistic; Gompertz MSC: 92D25; 34A12; 34A34 1. Introduction The logistic model has been extensively used to study the cause and effect relationship Citation: Rodrigo, M.; Zulkarnaen, D. between the ‘carrying capacity’ (i.e., the population size that available resources can support) Mathematical Models for Population and the population size. See, for instance, the references in Brauer & Castillo-Chávez [1], Growth with Variable Carrying Capacity: Analytical Solutions. Gotelli [2], Pastor [3], as well the seminal papers by Verhulst [4] and Pearl & Reed [5]. AppliedMath 2022, 2, 466–479. It is typically assumed that the carrying capacity does not vary with time and, therefore, https://doi.org/10.3390/ the logistic model exhibits a sigmoidal shape when the population size is plotted as a appliedmath2030027 function of time. However, many phenomena such as human population growth exhibit more complex behaviours, unlike population species grown in laboratory cultures, for Academic Editor: Leonid Shaikhet example. Thus, it is of theoretical and practical interest to investigate mathematical models Received: 6 July 2022 that incorporate variable carrying capacities. One of the main objectives of studying such Accepted: 18 August 2022 models is to explore how different functional forms of the carrying capacity influence the Published: 29 August 2022 dynamics of the population variable and its long-term behaviour. Meyer [6] and Meyer & Ausubel [7] studied a bi-logistic model derived from a logistic Publisher’s Note: MDPI stays neutral model but with a sigmoidal time-dependent carrying capacity. Cohen [8] proposed a with regard to jurisdictional claims in published maps and institutional affil- human population growth model with a variable carrying capacity that is a function of the iations. population size itself. The aforementioned models make the case that the inclusion of a variable carrying capacity is more reflective of the human condition. Safuan et al. [9] proposed a coupled system of ordinary differential equations (ODEs) to describe the interaction between the population and its carrying capacity. The model they Copyright: © 2022 by the authors. considered does not require prior knowledge of the carrying capacity, nor does it impose Licensee MDPI, Basel, Switzerland. constraints on the initial conditions. Assuming a special form of the carrying capacity in This article is an open access article the logistic equation, the same authors obtained an analytical solution in series form [10]. distributed under the terms and As pointed out by Cohen [8], there is no consensus with regards to appropriate models conditions of the Creative Commons for human carrying capacity. However, most would accept that human carrying capacity Attribution (CC BY) license (https:// is influenced by food availability, amongst other factors. Hopfenberg [11] postulated that creativecommons.org/licenses/by/ food production data are the sole variable that influences human carrying capacity, and 4.0/). AppliedMath 2022, 2, 466–479. https://doi.org/10.3390/appliedmath2030027 https://www.mdpi.com/journal/appliedmath AppliedMath 2022, 2 467 a simple linear relationship between human carrying capacity and the food production index is assumed. More recently, Zulkarnaen & Rodrigo [12] proposed three classes of human population dynamical models of logistic type where the carrying capacity is a function of the food production index. They employed an integration-based parameter estimation technique [13] to derive explicit formulas for the model parameters. Using actual population and food production index data, the results of numerical simulations of their models suggested that an increase in food availability implies an increase in carrying capacity, but the carrying capacity is ‘self-limiting’ and does not increase indefinitely. Thornley & France [14] proposed an ‘open-ended’ form of the logistic equation by considering a system of two ODEs representing the coupled processes of growth and development. Their model is ‘open-ended’ in the sense that dynamic changes in nutrition and environment can influence growth and development, which, in turn, may affect the asymptotic carrying capacity value. Subsequently, Thornley et al. [15] found an analytical solution to the Thornley–France model in the case of constant parameters. The solution of the system of ODEs is expressed in terms of the solution of a single ODE of power-law logistic type, also referred to as the q-logistic model, which frequently arises in ecology and elsewhere [16]. A related article by Wu et al. [17] formulated the variable carrying capacity by exploring a resource dynamic-based feedback mechanism underlying the population growth models. The inclusion of variable carrying capacities in interacting species such as in predator–prey models have been considered in Al-Moqbali et al. [18], Al-Salti et al. [19]. Power-law logistic models have been investigated by von Bertalanffy [20] and Richards [21]. The principal (nonnegative) parameter of these models is denoted by q. The Gompertz model and the logistic model are recovered when q = 0 and q = 1, respectively. Larger values of q behave like a logistic model but with an increasingly sharper cessation of growth as the asymptote (i.e., the constant carrying capacity limit) is approached [15]. As a fraction of the asymptote, inflexion can occur over the range from (Gompertz) through (‘ordinary’ logistic) and then to 1 (for large q). Determining the point where inflexion takes place can be especially important when fitting the model to actual data that exhibit a sigmoidal trend. See Banks [22] for a detailed analysis of the q-logistic model. More recently, Albano et al. [23] considered a general growth curve that includes the Malthus, Richards, Gompertz and other models. Their generalised model is essentially a Bernoulli ODE, which is well known to be analytically tractable. They investigated the analytical and numerical properties of the solution as the parameters varied. Here we propose a general population model with a variable carrying capacity, which includes the Thornley–France [14], Safuan–Jovanoski–Towers–Sidhu [9] and Meyer– Ausubel [7] models as special cases. Moreover, when the carrying capacity is kept constant, the proposed model system reduces to a single ODE population model and recovers the Gompertz, ‘ordinary’ logistic and q-logistic models, amongst others. The idea is to extract the essential properties of such models without getting ‘bogged down’ by particular cases. We provide a procedure for obtaining, when possible, the analytical solution of this general population model. Different models can then be chosen depending on the particular phe- nomenon being studied. An important tractable special case is when the per capita growth rates of the population and carrying capacity are proportional to each other. We give a criterion for when inflexion may occur. Several illustrative examples are also provided. 2. Population Models with Variable Carrying Capacities Consider the initial value problem (IVP) dN dK = N f (N, K), = K g(N, K), N(0) = N , K(0) = K , (1) 0 0 dt dt where N(t) and K(t) are the population and carrying capacity, respectively, at time t. The initial values N and K are positive and given. Let f and g be sufficiently smooth 0 0 bivariate functions on R  R , where R = (0,¥). Denote by D , where j = 1, 2, the + + + partial derivative with respect to the independent variable in the jth position. We assume AppliedMath 2022, 2 468 1 dN that the population per capita growth rate decreases with increasing population N dt (D f (x, y) < 0) and increases with increasing carrying capacity (D f (x, y) > 0). As 1 2 for the signs of D g(x, y) and D g(x, y), it is not obvious what these should be since 1 2 1 dK the behaviour of the carrying capacity per capita growth rate may depend on the K dt particular population species. When g is identically zero, then K(t) = K and (1) reduces to dN = N f (N, K ), N(0) = N . 0 0 dt A well-known example is the q-logistic model [22] r N f (N, K ) = 1 , q K where r > 0 is the intrinsic growth rate, and q > 0 is a parameter related to the point of inflexion of the solution. The case q = 1 gives the ‘ordinary’ logistic model. In the limit as q ! 0 , the q-logistic model reduces to the Gompertz model [24] f (N, K ) = r log . In the case of a variable carrying capacity, the Thornley–France model [14] takes the form N N (2) f (N, K) = a 1 , g(N, K) = b 1 , K K where a, b > 0. Note that D g(x, y) > 0 and D g(x, y) < 0. Safuan et al. [9] proposed 1 2 the model (3) f (N, K) = a 1 , g(N, K) = b cN, where a, b, c > 0. Here we see that D g(x, y) < 0 and D g(x, y) = 0. Meyer [6] and Meyer & Ausubel [7] assumed that f (N, K) = a 1 , g(N, K) = b cK, where a, b, c > 0. This time, D g(x, y) = 0 and D g(x, y) < 0. 1 2 2.1. Construction of Analytically Tractable Population Models with Variable Carrying Capacities Suppose that there exists a positive C -function F such that y = F(x) solves the IVP yg(x, y) dx x f (x, y) dy = 0, y = K when x = N . (4) 0 0 Define G(x) = dz. (5) z f (z, F(z)) We claim that G(N(t)) = t, K(t) = F(N(t)) (6) gives the formal solution of (1) in implicit form. Indeed, (4) and (6) imply that F(x)g(x, F(x)) K(t)g(N(t), K(t)) 0 0 F (x) = , F (N(t)) = . x f (x, F(x)) N(t) f (N(t), K(t)) AppliedMath 2022, 2 469 Implicitly differentiating G(N(t)) = t with respect to t yields 0 0 0 1 = G (N(t))N (t) = N (t) N(t) f (N(t), K(t)) dN or = N f (N, K). Hence dt K(t)g(N(t), K(t)) 0 0 0 0 K (t) = F (N(t))N (t) = N (t) = K(t)g(N(t), K(t)) N(t) f (N(t), K(t)) dK or = K g(N, K). It is clear that G(N(0)) = G(N ) = 0 and K(0) = F(N(0)) = F(N ) = 0 0 dt K . This proves the claim. Therefore, the task of finding an analytical solution of (1) essentially reduces to solving the first-order ODE (4). 2.2. Determination of Inflexion Points for the Population Species Variable Before we consider some examples, let us first investigate where inflexion may occur for the population species variable. As mentioned previously, this is particularly relevant during model fitting. Using (6) and differentiating the first ODE in (1) with respect to t, we have d N dN = [ f (N, F(N)) + N D f (N, F(N)) + N D f (N, F(N))F (N)] . 1 2 dt dt Therefore, inflexion for the function N may occur at some positive root N of the equation H(z) = f (z, F(z)) + zD f (z, F(z)) + zD f (z, F(z))F (z) = 0. (7) 1 2 If this occurs, it will be when t = t , where t = G(N(t )) = G(N ) = dz. (8) z f (z, F(z)) 2.3. First Class of Analytically Tractable Models: f and g Are Proportional Suppose that there exists a 2 Rnf0g such that g(x, y) = a f (x, y). This basically as- 1 dN 1 dK sumes that the per capita growth rates and are proportional. Then (4) and (5) give N dt K dt a a y = F(x) = K N x , G(x) = dz, N z f (z, K N z ) 0 0 respectively. From (6) we obtain the exact solution N(t) a a dz = t, K(t) = K N [N(t)] (9) a 0 N z f (z, K N z ) 0 0 of the system (1). Example 1. Let q q a x b x f (x, y) = 1 , g(x, y) = 1 , (10) q y q y where a, b > 0 and q  0. Then g(x, y) = a f (x, y), where a = < 0. Thus, (1) becomes q q dN a N dK b N = N 1 , = K 1 , N(0) = N , K(0) = K . (11) 0 0 dt q K dt q K AppliedMath 2022, 2 470 The ‘value’ when q = 0 is meant to be understood as the limit when q ! 0 ; this is related to the Gompertz model in the case of a constant carrying capacity. The Thornley–France model (2) is a special case of (10) if we take q = 1. From the first equation of (9), we obtain Z Z a q q(1a) N(t) (K N ) [N(t)] q 1 1 0 1 t = dz = du. a q(1a) q q(1a) q a a(1 a) u(1 u) N z[1 (K N ) z ] (K N ) N 0 0 0 0 0 (12) Recall the formula 1 u I(u) = du = log , 0 < u < 1. u(1 u) 1 u If I(b ) I(a ) = c , where 0 < a < b < 1 and c > 0, then it is not difficult to show that 0 0 0 0 0 0 b = . 1a 0 c 1 + e q(1a) a q 1 q a q q(1a) Taking a = (K N ) N = (K N ) , b = (K N ) [N(t)] and 0 0 0 0 0 0 0 0 0 c = a(1 a)t, we deduce that the exact solution of (11) is " # q(1a) a q (K N ) 0 a a N(t) = , K(t) = K N [N(t)] . (13) 1 0 1(K N ) 0 a(1a)t 1 + e (K N ) As a, b > 0 and a = < 0, we have from (13) that 1 1 a a 1a 1a N = lim N(t) = (K N ) , K = lim K(t) = (K N ) . ¥ 0 ¥ 0 0 0 t!¥ t!¥ This is consistent with the fact that f (N , K ) = g(N , K ) = 0 in (11) yields that any ¥ ¥ ¥ ¥ equilibrium of (11), and there are infinitely many, lies on the curve K = N . ¥ ¥ Let us now look at the inflexion points of N by studying the roots of H(z) = 0 in (7). Straightforward calculations give q q a x a x D f (x, y) = , D f (x, y) = , xF (x) = aF(x), x y y y xD f (x, y) = a + q f (x, y), yD f (x, y) = a q f (x, y). 1 2 If we define H (x, y) = f (x, y) + xD f (x, y) + ayD f (x, y), 0 1 2 then H(z) = 0 if and only if H (z, F(z)) = 0. However, H (x, y) = [1 + q(1 a)] f (x, y) 0 0 a(1 a), so that H (z, F(z)) = 0 if and only if a z a(1 a) 1 = f (z, F(z)) = . q F(z) 1 + q(1 a) This is equivalent to q 1a 1 z z z 1 1a = = = , K = (K N ) . ¥ 0 1 + q(1 a) F(z) K K N z ¥ Thus H(N ) = 0 if and only if q(1a) N = K . (14) 1 + q(1 a) AppliedMath 2022, 2 471 Finally, we see from (8) that q(1a) q(1a) 1 N N t = I I . (15) a(1 a) K K ¥ ¥ As I is an increasing function of u, and assuming that N < N , it follows that t > 0. It was already noted that q = 1 in (10) reduces to the Thornley–France model (2). Then, if q = 1, (13) recovers the analytical solution found in Thornley et al. [15]. Equation (14) indicates where inflexion occurs as a fraction of the asymptotic carrying capacity, while (15) gives the time of inflexion (compare with (17) and (18), respectively, found in Thornley et al. [15] with an appropriate renaming of parameters). On the other hand, when we let q ! ¥, then (14) shows that tends to unity so that, similar to the Thornley–France model when q = 1, exponential growth is sustained for longer and the inflexion value N moves closer to the asymptotic carrying capacity value K . If q > 0 and a = 0 (corresponding to b = 0 in (10)), then K(t) = K and (1) reduces to a single ODE dN a N = N 1 , dt q K which is the q-logistic model. In particular, if q = 1 and a = 0 (i.e., the ‘ordinary’ logistic model), then (14) shows that the inflexion value is one-half of the asymptotic carrying capacity value, which + N 1 b is well known. If q ! 0 , then we deduce from (14) that tends to for any a = < 0. K e a This is similar in behaviour to the case when a = 0 (i.e., b = 0 in (10)), so that K(t) = K and (1) recovers the Gompertz model dN K = aN log . dt N N 1 + Thus tends to as q ! 0 even for the variable carrying capacity system (11). K e Here we include some results of numerical simulations of (11) (equivalently, (13)). Let a = 0.2, b = 0.1, N = 5 and K = 20. Take four different values of q = 0.01, 1, 5, 20 to see and compare the 0 0 population dynamics, as shown in Figure 1, as well as the carrying capacity behaviour in Figure 2. As can be seen in both figures, the larger the value of q, the longer it takes for the population and carrying capacity to reach equilibrium. 20.000 13.333 6.667 inflexion point 0.000 0.000 50.000 100.000 150.000 Figure 1. Profiles of N(t) in (13) for q = 0.01, 1, 5, 20 (with a = 0.2, b = 0.1, N = 5 and K = 20) and 0 0 inflexion points calculated from (14) and (15). AppliedMath 2022, 2 472 20.000 13.333 6.667 0.000 0.000 50.000 100.000 150.000 Figure 2. Profiles of K(t) in (13) for q = 0.01, 1, 5, 20 (with a = 0.2, b = 0.1, N = 5 and K = 20). 0 0 Figure 3 (q = 0.01) and Figure 4 (q = 1.2) depict some trajectories in the NK-plane and illustrate that the carrying capacity experiences a decline as the population grows larger. This is a consequence of the proportionality of the per capita growth rates with a negative proportionality constant. Moreover, when the initial conditions are varied, the trajectories approach the equilibrium line K = N . This behaviour is typical for other values of q as well. ¥ ¥ 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 Figure 3. Representative solution trajectories of (11) in the NK-plane with q = 0.01, a = 0.2 and b = 0.1. AppliedMath 2022, 2 473 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 Figure 4. Representative solution trajectories of (11) in the NK-plane with q = 1.2, a = 0.2 and b = 0.1. 2.4. Second Class of Analytically Tractable Models: f and g Are Homogeneous Functions Suppose that f and g in (1) are homogeneous functions, i.e., f (lx, ly) = l f (x, y), g(lx, ly) = lg(x, y) for every l > 0. Let v = in (4), so that g(x, y) = g(x, xv) = xg(1, v) and f (x, y) = f (x, xv) = x f (1, v). Then (4) is transformed into the separable ODE dv g(1, v) f (1, v) x = v , dx f (1, v) whose general solution is y f (1, v) F(v) = F = log(x) + log(C), F(v) = dv, (16) x v[g(1, v) f (1, v)] where C > 0 is an arbitrary constant of integration. Hence, assuming that F is invertible, we deduce that the solution of (4) is 1 0 y = F(x) = xF (log(Cx)), log(C) = F log(N ). (17) Next, we determine G. From (5), we see that G(x) = dz. (18) z f (z, zF (log(Cz))) Finally, (6) implies that the exact solution of (1) in the case when both f and g are homoge- neous functions is N(t) dz = t, K(t) = N(t)F (log(C N(t))). (19) z f (z, zF (log(Cz))) Example 2. Let x cx f (x, y) = a 1 , g(x, y) = a , y y AppliedMath 2022, 2 474 where 0 < a < c, so that (1) becomes dN N dK cN = aN 1 , = K a , N(0) = N , K(0) = K ; (20) 0 0 dt K dt K compare with model (3) proposed by Safuan et al. [10]. Observe that f and g are homogeneous functions but are not proportional since a < c. Then (16) gives h  i a y a y y F(v) = [v log(v)], F = log = log(Cx) (21) a c x a c x x and, therefore, c K K K 0 0 0 a 0 1 log(C N ) = log or (C N ) = e . (22) 0 0 a N N N 0 0 0 The second equation in (21) can be rewritten as y a c c x a e = exp log(Cx) = (Cx) . (23) x a Let us recall some properties of the Lambert W-function. Suppose that u, z 2 R. Then z = W(u) 1 1 satisfies the equation ze = u if u  . Furthermore, W(u) < 0 for  u < 0 and e e W(u) 1 W (u) = , u 6= , 0. (24) u[1 + W(u)] e Identifying z = , u = (Cx) , 1 1 ac and assuming that 0 < x  e so that u  , we obtain from (23) that C e y = F(x) = xW((Cx) ). We see that F(x) > 0 since the W-function is negative here because  u < 0. Following a similar argument, we deduce from (21) that ac c 1 z 1 1 a a F (z) = W(e ), F (log(Cz)) = W((Cz) ). Therefore, (18) yields Z Z 1 1 a x (Cx) 1 W((Cz) ) 1 W(u) G(x) = dz = du 1 1 a a c a a u[1 + W(u)] N (C N ) 0 z[1 + W((Cz) )] 0 1 c c 1 1 a a = [W((Cx) ) W((C N ) )], c a where we used (24) in the last step. Equation (22) gives c K K 0 0 W((C N ) ) = W e = N N 0 0 z z since z = W(u) = W(ze ) if u = ze . Hence the exact solution of (20) from (19) is c K c 1 1 a a W([C N(t)] ) = + (c a)t, K(t) = N(t)W([C N(t)] ). (25) Example 3. Suppose that y y f (x, y) = a log , g(x, y) = a log + b, x x AppliedMath 2022, 2 475 where a, b > 0. Therefore, (1) becomes dN K dK K = aN log , = K a log + b , N(0) = N , K(0) = K . (26) 0 0 dt N dt N We see that f and g are homogeneous functions but are not proportional as b > 0. Then (16) implies that h  i a y a y F(v) = [log(v)] , F = log = log(Cx) 2b x 2b x and 2 1 a K 1 a K 0 0 log(C) = log log(N ), [log(C N )] = log . (27) 0 0 2b N 2b N 0 0 Since 1 1 h i 2b 2 2b 2 1 1 F (v) = exp v , F (log(Cv)) = exp log(Cv) , a a Equation (17) yields h i 2b y = F(x) = x exp log(Cx) . Furthermore, h i 2b 1 1 2 2 f (z, zF (log(Cz))) = a log(Cz) = (2ab) [log(Cz)] and from (18), we obtain Z Z x log(Cx) 1 1 1 G(x) = dz = u du 1 1 1 N log(C N ) 2 2 2 (2ab) 0 z[log(Cz)] (2ab) 0 1 1 2 1 2 1 2 2 2 2 = [log(Cx)] [log(C N )] . ab ab 1 1 1 ab 2 2 2 Substituting into (19) gives [log(C N(t))] = [log(C N )] + ( ) t, which simplifies to K ab log(C N(t)) = log(C N ) + at log + t N 2 using (27). Thus, at ab 2 N(t) = N e . Moreover, (19) also implies that h i 2b 2 K(t) = N(t) exp log(C N(t)) , so that h i 2b log(K(t)) = log(N(t)) + log(C N(t)) 1 1 1 2b 1 2b ab 2 2 2 = log(N(t)) + [log(C N )] + t a a 2 = log(N(t)) + log + bt. 0 AppliedMath 2022, 2 476 Hence at at K K K ab 2 K ab 2 0 0 0 0 bt bt t t +bt 2 2 K(t) = e N(t) = e N e = K e . 0 0 N N N N 0 0 0 0 Thus, the exact solution of (26) is at at K ab 2 K ab 2 0 0 t t +bt 2 2 N(t) = N e , K(t) = K e . (28) 0 0 N N 0 0 2.5. Third Class of Analytically Tractable Models: f and g Determine an Exact ODE Assume that f and g are such that ¶ ¶ [x f (x, y)] + [yg(x, y)] = f (x, y) + xD f (x, y) + g(x, y) + yD g(x, y) = 0. (29) 1 2 ¶x ¶y This implies that (4) is an exact ODE whose general solution is Y(x, y) = C = Y(N , K ), 0 0 where C is an arbitrary constant of integration and D Y(x, y) = yg(x, y), D Y(x, y) = x f (x, y). (30) 1 2 Then Z Z Y(x, y) = yg(x, y) dx + Q(y) = x f (x, y) dy + P(x), (31) where P = P(x) and Q = Q(y) are arbitrary functions of integration. Differentiating Y in (31) with respect to y and using the second equation in (30) and then (29), we get Z Z Z Q(y) = x f (x, y) f (x, y) dx xD f (x, y) dx dy. (32) An analogous calculation yields Z Z Z P(x) = yg(x, y) f (x, y) dy xD f (x, y) dy dx. (33) Assuming that Y(x, y) = C = Y(N , K ) can be solved for y = F(x), then we have 0 0 from (5) and the second equation in (30) that Z Z x x 1 1 G(x) = dz = dz. (34) z f (z, F(z)) D Y(z, F(z)) N N 2 0 0 Therefore, (6) expresses the exact solution of (1) as N(t) dz = t, K(t) = F(N(t)). (35) D Y(z, F(z)) N 2 Looking at condition (29) more closely, we conclude that Z Z yg(x, y) = f (x, y) dy xD f (x, y) dy + R(x) (36) 1 AppliedMath 2022, 2 477 for some arbitrary function R = R(x). Hence the exactness condition (29) necessarily implies the above form for g. Thus from (31) and (36), we have Z Z Z Y(x, y) = f (x, y) dy + x D f (x, y) dy R(x) dx Z Z Z + x f (x, y) f (x, y) dx xD f (x, y) dx dy Z Z = x f (x, y) dy R(x) dx. This simplified form for Y is then substituted into the first equation in (35). Example 4. Let a > 0. If 1 1 f (x, y) = a 1 = a axy , D f (x, y) = ay , then (36) evaluates to yg(x, y) = ay + 2ax log(y), taking R(x) = 0 for simplicity. Therefore, (1) becomes dN N dK = aN 1 , = aK + 2aN log(K), N(0) = N , K(0) = K . (37) 0 0 dt K dt We have Y(x, y) = x(a axy ) dy = ax[y x log(y)], where C = Y(N , K ) = aN [K N log(K )]. The equation Y(x, y) = C gives 0 0 0 0 0 0 y C y C 2 1 2 log(y) = x , e = x exp x . x a x a Hence, as in Example 2, C C 1 1 2 1 2 y = F(x) = xW x exp x if x exp x  . a a e Then h i 2 1 D Y(x, y) = ax ax y , D Y(z, F(z)) = az 1 2 2 F(z) and 1 C 2 1 + W(z exp( z ) 1 = . 1 2 F(z) W(z exp( z ) From (34) we obtain 1 2 W(z exp( z ) G(x) = dz, 1 2 z[1 + W(z exp( z )] which cannot be further evaluated analytically (unlike in Example 2). The exact solution of (37) using (35) is 1 C 2 N(t) W(z exp( z ) dz = at, 1 2 z[1 + W(z exp( z )] (38) 1 2 K(t) = N(t)W [N(t)] exp [N(t)] . a AppliedMath 2022, 2 478 3. Concluding Remarks In this article, we proposed population models with variable carrying capacities modelled by a coupled system (1) of two nonlinear ODEs and found their analytical solutions. While it was clear that the assumptions D f (x, y) < 0 and D f (x, y) > 0 are 1 2 reasonable since they describe the behaviour of the population per capita growth rate, we showed through several explicit examples that corresponding assumptions for D g(x, y) and D g(x, y) that describe the carrying capacity per capita growth rate are not obvious and may be model dependent. One possible reason for this is that carrying capacity is not directly observable, unlike the population size. If the per capita growth rates are proportional, as in Example 1, then in addition to the analytical solution, we also found a criterion for the occurrence of inflexion in the population profile as a fraction of the asymptotic carrying capacity. This criterion does not apply to the models in Examples 2 and 3 since they do not have a nontrivial equilibrium. However, it should apply for the model in Example 4 with a nontrivial 1 1 2 2 equilibrium (N , K ) = (e , e ). Further classes of analytically tractable models of ¥ ¥ the form (1), not necessarily in the context of population growth with variable carrying capacity, can, of course, be found by considering analytically tractable cases of (4) following the idea in Rodrigo [25]. Future work would involve a more careful investigation of the assumptions for the carrying capacity per capita growth rate. While this article focused on the analytical aspects of population models with variable carrying capacities, fitting the models to actual population data is the next important step. The estimation of parameters in the models can be undertaken by adapting the arguments in Zulkarnaen & Rodrigo [12], Holder & Rodrigo [13]. Parameter estimation using an integration-based technique in population models with variable carrying capacities that depend on food availability is the subject of an article by the authors that is currently under review [26]. Another research direction is the modelling of interacting population species where the carrying capacities are not fixed anymore but may depend on time and/or space. Results for predator–prey models have been obtained in Al-Moqbali et al. [18], Albano et al. [23], although spatio-temporal models (e.g., chemotaxis [27]) can also be considered. Travelling waves are an important class of biologically relevant solutions for spatio-temporal models described by partial differential equations. A travelling wave coordinate transformation leads to a system of higher dimensional ODEs, for which analytically tractable models can potentially be identified using the techniques in the current article. Author Contributions: Conceptualization, M.R.; Formal analysis, M.R.; Investigation, M.R.; Method- ology, M.R.; Project administration, M.R.; Software, M.R. and D.Z.; Supervision, M.R.; Validation, M.R.; Visualization, M.R. and D.Z.; Writing—original draft, M.R. and D.Z.; Writing—review & editing, M.R. and D.Z. All authors have read and agreed to the published version of the manuscript. Funding: This research received no external funding. Data Availability Statement: Not applicable. Conflicts of Interest: The authors declare no conflict of interest. References 1. Brauer, F.; Castillo-Chávez, C. Mathematical Models in Population Biology and Epidemiology, 2nd ed.; Springer: New York, NY, USA, 2012. 2. Gotelli, N. A Primer of Ecology, 2nd ed.; Sinauer Associates: Sunderland, MA, USA, 1998. 3. Pastor, J. Mathematical Ecology of Populations and Ecosystems; Wiley-Blackwell: London, UK, 2008. 4. Verhulst, P.F. Notice sur la loi que la population suit dans son accroissement. Corresp. Mathématique Phys. 1838, 10, 113–121. 5. Pearl, R.; Reed, L.J. A further note on the mathematical theory of population growth. Proc. Natl. Acad. Sci. USA 1922, 8, 365–368. [CrossRef] 6. Meyer, P.S. Bi-logistic growth. Technol. Forecast. Soc. Chang. 1994, 47, 89–102. [CrossRef] 7. Meyer, P.S.; Ausubel, J.H. Carrying capacity: A model with logistically varying limits. Technol. Forecast. Soc. Chang. 1999, 61, 209–214. [CrossRef] AppliedMath 2022, 2 479 8. Cohen, J.E. Population growth and the Earth’s human carrying capacity. Science 1995, 269, 341–346. [CrossRef] 9. Safuan, H.M.; Jovanoski, Z.; Towers, I.N.; Sidhu, H.S. Coupled logistic carrying capacity. ANZIAM J. 2012, 53, 172–184. [CrossRef] 10. Safuan, H.M.; Jovanoski, Z.; Towers, I.N.; Sidhu, H.S. Exact solution of a non-autonomous logistic population model. Ecol. Model. 2013, 251, 99–102. [CrossRef] 11. Hopfenberg, R. Human carrying capacity is determined by food availability. Popul. Environ. 2003, 25, 109–117. [CrossRef] 12. Zulkarnaen, D.; Rodrigo, M.R. Modelling human carrying capacity as a function of food availability. ANZIAM J. 2020, 62, 318–333. [CrossRef] 13. Holder, A.B.; Rodrigo, M.R. An integration-based method for estimating parameters in a system of differential equations. Appl. Math. Comput. 2013, 219, 9700–9708. [CrossRef] 14. Thornley, J.H.M.; France, J. An open-ended logistic-based growth function. Ecol. Model. 2005, 184, 257–261. [CrossRef] 15. Thornley, J.H.M.; Shepherd, J.J.; France, J. An open-ended logistic-based growth function: analytical solutions and the power-law logistic model. Ecol. Model. 2007, 204, 531–534. [CrossRef] 16. Gilpin, M.E.; Ayala, F.J. Global models of growth and competition. Proc. Natl. Acad. Sci. USA 1973, 70, 3590–3593. [CrossRef] 17. Wu, H.; Chakraborty, A.; Li, B.-L.; Kenerley, C.M. Formulating variable carrying capacity by exploring a resource dynamics-based feedback mechanism underlying the population growth models. Ecol. Complex. 2009, 6, 403–412. [CrossRef] 18. Al-Moqbali, M.K.A.; Al-Salti, N.S.; Elmojtaba, I.M. Prey-predator models with variable carrying capacity. Mathematics 2018, 6, 102. [CrossRef] 19. Al-Salti, N.; Al-Musalhi, F.; Gandhi, V.; Al-Moqbali, M.; Elmojtaba, I. Dynamical analysis of a prey-predator model incorporating a prey refuge with variable carrying capacity. Ecol. Complex. 2021, 45, 100888. [CrossRef] 20. von Bertalanffy, L. Quantitative laws for metabolism and growth. Q. Rev. Biol. 1957, 32, 217–231. [CrossRef] 21. Richards, F.J. A flexible growth function for empirical use. J. Exp. Bot. 1959, 10, 290–300. [CrossRef] 22. Banks, R.B. Growth and Diffusion Phenomena: Mathematical Frameworks and Applications; Springer: Berlin/Heidelberg, Germany, 1994. 23. Albano, G.; Giorno, V.; Román-Román, P.; Torres-Ruiz, F. Study of a general growth model. Commun. Nonlinear Sci. Numer. Simul. 2022, 107, 106100. [CrossRef] 24. Gompertz, B. On the nature of the function expressing the law of human mortality and on a new mode of determining the value of life contingencies. Philos. Trans. R. Soc. Lond. 1825, 115, 513–585. 25. Rodrigo, M.R. An elementary method for obtaining the general solution of a system of ordinary differential equations. Electron. J. Differ. Equ. 2021, 2021, 1–20. 26. Zulkarnaen, D.; Rodrigo, M.R. Population models with variable carrying capacities that depend on food availability. 2022, under review. 27. Hillen, T.; Painter, K.J. A user ’s guide to PDE models for chemotaxis. J. Math. Biol. 2009, 58, 183–217. [CrossRef] [PubMed]

Journal

AppliedMath – Multidisciplinary Digital Publishing Institute

Published: Aug 29, 2022

Keywords: population model; variable carrying capacity; exact solution; power-law logistic; logistic; Gompertz

There are no references for this article.