Population Models

---

Growth or decay is a common occurrence in many systems ranging from biology to economics. Since mathematical language is more concise and less ambiguous than any natural language, the biological models were revised leading to the proposal of a variety of growth curves and the theoretical frameworks that lay behind them. Due to environmental trends and fluctuation, probably no biologist believes that one equation would suit all growth processes. This situation is completely different from those that we observe in physics, where usually only one law is used. For instance, modeling of the fall of a body in a vacuum requires only one of Newton's laws. In many systems, growth patterns are ultimately subjected to some form of limitation. In practice, we usually observe a set of measured values Y_t, at some time sequence { t }, called a time series. Modeling these observations means finding a curve P(t) that is closed to the time series. For example, to model a US population growth, we deal with an equally spaced time series, based on the census, conducted every ten years.

In this section, we discuss some simple models that are expressed solely in terms of differential equations of the first order that ignore random fluctuations in the population. Such models are referred to as deterministic and they are written as

where P(t) is the size of population at time t and r is the rate of change. A number of curve models have been developed to describe and predict such behaviors. For example, Kiviste (1988) described 75 of them in a two-volume monograph. See also the comprehensive work by Boris Zeide. There are two basic trends in population models: either autonomous differential equations or separable equations: In many practical cases, these two approaches lead to the same population curve. Also, with different parametrizations, the same model may be written in distinct (equivalent) ways. Mathematical models to describe growth phenomena can be useful in many fields of interest including biology, ecology, medicine, or economics. Another important example gives a tumor growth, and, in particular, cancer growth. Although cancer as well as many other tumor diseases are uncontrolled cell proliferation, developing their models may contribute to our more insightful description of changing phenomena.

 

With developing mathematical methods, biologists started in the late 18th century to use these techniques in population modeling in order to understand the dynamics of growing and shrinking populations of living organisms. Thomas Malthus (1766--1834) was one of the first to note that populations grow with a geometric pattern (dP/dt = r P) while contemplating the fate of humankind. The Malthus model leads to exponential (unbounded) growth, which is clearly not viable in the long-term run:

In 1960, H. von Foerster, P.M. Mora, and L.W.Amiot proposed a Doomsday model:

which has the solution From the formula above, a very disturbing conclusion follows: there is a finite time t = N^-k k/r at which the population size becomes infinite. Therefore, the doomsday model gives more striking prediction of future catastrophe than even Malthus could come up with.

Von Foerster and his colleagues present many examples of accepted models of physical systems with singularities: pressure as a function of velocity approaching the speed of sound, current approaching breakdown voltage in gaseous conduction, index of refraction in optical absorption bands, and magnetic susceptibility at Curie temperature in ferromagnetism. It is reasonable to expect a singular behavior of a social structure subjected to extreme overpopulation. Nevertheless, the doomsday model provides a remarkably good fit of practical examples for certain time intervals. Example: Exponential growth model for population

Another fundamental growth mode is one subjected to limitations that affect growth rates, making growth slow down over time and asymptotically strive towards a maximum value. Growth may actually continue interminably but the rate of growth approaches zero as time tends to infinity. This is well-known in many biological systems where an organism may grow quickly in its juvenile stage but then growth slows down as it reaches maturity. Bounded growth curves can be used to model renewable energy sources since some kind of upper limit must exist even for renewable energy. Rather in the growth process itself, the limiting factor lies in the form of the increasing monetary, energy, or resource costs required for continued expansion. The upper limit may be high, or virtually non-existent, but the steps in development become more and more challenging to take, thus slowing the growth process. In the biological sciences, this is often seen as proportionality between growth rate and actual size or age.

 

The best example of exponential growth is seen in bacteria. Bacteria are prokaryotes that reproduce by prokaryotic fission. This division takes about an hour for many bacterial species. If we use E. coli bacteria, we could start with just one bacterium and have enough bacteria to cover the Earth with a 30-cm layer in just 36 hours.

Usually, we never observe such exponential growth in long time run because of limitation of resources. Therefore, the Malthus model can be used only for relatively short period of time.    ■

One of the most basic and milestone models of population growth was the logistic model of population growth formulated by Pierre François Verhulst in 1838. (We don't know why Verhulst called this equation logistic, but this name is universally accepted.) An outstanding exposition of its history is given by G. Evelyn Hutchinson (1903--1991) in his book, published in 1978. The logistic model (\( \dot{P} = r\,P \left( 1 - P/K \right) \) ) takes the shape of a sigmoid curve and describes the growth of a population. He observed that if species are not preyed upon, its population is expected to saturate. Its solution is

Here t is the number of days since the starting day, t₀. The population size can also be written as where b is the fitting coefficient: The logistic model can be also written in the form: The relative growth rate (1/P(dP/dt) declines linearly (no maximum) with increasing population size and reaches its zero minimum at P = K, the carrying capacity. The population at the inflection point (where growth rate is maximum) is exactly half the carrying capacity. Example: Logistic equation for US population

 

Since that time, there has been made substantial progress in developing mathematical population models. Obviously, we cannot discuss and analyze a huge variety of population models. Instead, we present some famous models involving one ordinary differential equation. More complicated cases involving interaction of spices require a usage of system of differential equations, and this topic is covered in the second part of the tutorial.

Example: Logistic equation with predation

The following two population models mimic features of asexual and sexual reproduction. The asexual reproduction model is of the form:

The sexual reproduction model is of the form: We plot direction fields for these two differential equations.
The above direction field plots allow us to determine the behavior of solutions near critical points; some of them are asymptotically stable, some only semi-stable, and others are unstable. In each case, the variable x represents a scaled population size and the parameter 𝑎 is assumed to be positive. Both slope functions have the critical point x=0. However, the function f has an additional equilibrium point x=1-1/ 𝑎, while the function g may have two more equilibrium points: So we see that the slope function f has the bifurcation point 𝑎=1, but the slope function g has another bifurcation point 𝑎=4. Their derivatives are If 𝑎>1, then the derivative of f is positive and the critical point 1-1/𝑎 is a source (unstable). However, for 𝑎<1, the critical point is a sink (stable).

The slope function g(x) has two critical points, and the derivative at point \( x_{-}^{\ast} \) is negative for 𝑎>4; therefore, this point is asymptotically stable. The derivative of g(x) at another critical point \( x_{+}^{\ast} \) is positive for 𝑎>4. On the other hand, these two critical points disappear for 𝑎<4. So the point 𝑎=4 is a saddle-node bifurcation point.    ■

 

Consider the dynamics of a harvested population for the logistic equation

where r, N > 0, and where h(P) is the harvesting rate. We consider several important cases starting with a constant harvesting.

Constant harvesting

---

The two simplest harvesting models are constant harvesting rate and proportional harvesting rate. We start with the former. Suppose \( h(P) = h_0 \) is a constant. Then the slope function becomes

Setting its derivative to zero, we get \( P = \frac{1}{2}\, K . \) The function f(P) is a downward-opening parabola whose maximum value is \( f(K/2) = rK/4 - h_0 . \) Thus, if \( h_0 > rK/4 , \) the harvesting rate is too large and the population always shrinks. A saddle-node bifurcation occurs at this value of h₀, and for larger harvesting rates, there are fixed points at with P_- unstable and P₊ stable. When r K < 4 h₀, there are no real critical points, as the following direction fields support.

StreamPlot[{1, 1*x*(1 - x) - 1}, {t, 0, 3}, {x, -1, 2}, FrameLabel -> {t, P}, MeshStyle -> {Thick}, PlotTheme -> "Business"] StreamPlot[{1, 5*x*(1 - x) - 1}, {t, 0, 3}, {x, -0.4, 1.2}, FrameLabel -> {t, P}, MeshStyle -> {Thick}, PlotTheme -> "Web"]

By rescaling the population u = P/K, time \( \tau = r\,t , \) and harvesting rate \( h = h_0 /rK , \) we arrive at the equation The critical harvesting rate is then \( h_c = 1/4 . \)

The proportional harvesting model is, after non-dimensionalization, of the form:

T = 6.0; (* time interval length *) ER = 8/3.; (* extraction rate *) L = 16; (* carrying capacity *) k = 1/3; (* exponential rate *) b = 2.; (* low threshold *) y1[1] = 15.; (* initial population *) nT = 3; (* number of time cycles *) Do[ soln1[i] = DSolve[{y'[x] == k y[x] (1 - y[x]/L) - ER, y[(2 i - 2) T] == y1[i]}, y[x], x] // Quiet // Flatten; y2[i] = y[x] /. soln1[i] /. x -> T (2 i - 1); soln2[i] = DSolve[{y'[x] == k y[x] (1 - y[x]/L), y[(2 i - 1) T] == y2[i]}, y[x], x] // Quiet // Chop // Flatten; y1[i + 1] = y[x] /. soln2[i] /. x -> T (2 i); , {i, 1, nT}]; mytable = Table[{{y[x] /. soln1[i], (2 i - 2) T <= x < (2 i - 1) T}, {y[x] /. soln2[i], (2 i - 1) T <= x < (2 i) T}}, {i, 1, nT}]; myfun[x_] = Piecewise[Partition[Flatten[mytable], 2]]; Plot[{b, L, myfun[x]}, {x, 0, nT 2 T}, PlotRange -> {{0, nT 2 T}, {0.1, 1.1 L}}, Frame -> True, PlotStyle -> {Dashing[0.02], Dashing[0.02], Automatic}, FrameTicks -> {T Range[0, 2 nT], Automatic, None, None}, AxesOrigin -> {0, 0}]

 

One defect of the constant harvesting rate model is that P = 0 is not a fixed point. To remedy this, consider the following model for the harvesting rate:

where 𝑎 and b are positive constants. We can rescale (N,t) to (n,τ) such that where γ and c are positive constants. To examine the denominator of h(p), we must take P = b n. Dividing both sides of \( \dot{P} = f(P) \) by 𝑎, we obtain from which we glean τ = 𝑎t/b,   γ = rb/𝑎, and   c = N/b. Now we examine the scaled population model

StreamPlot[{1, 2*x*(1 - x) - x*x/(1 + x*x)}, {t, 0, 3}, {x, -0.4, 1}, FrameLabel -> {t, P}, StreamStyle -> Arrowheads[{{0.03, Automatic}}], MeshShading -> {Red, None}]

There is an unstable fixed point at n = 0, where g'(0) = γ > 0. The other fixed points occur when the term in the brackets vanishes. We seek the intersection of the harvesting function \( h(n) = n/(n^2 + 1) \) with a two-parameter family of straight lines, given by \( y(n) = \gamma \left( 1- n/c \right) . \) The n-intercept is c and the y-intercept is γ. Provided c > c^* is large enough, there are two bifurcations as a function of γ, which we call \( \gamma_{\pm} (c) . \)

Both bifurcations are of the saddle-node type. We determine the curves \( \gamma_{\pm} (c) \) by requiring that h(n) is tangent to y(n), which gives two equations:

Together, these give γ(c) parametrically, i.e. as

ParametricPlot[{2*n^3/(n^2 - 1), 2*n^3/(n^2 + 1)^2}, {n, -0.6, 0.6}, Axes -> False, PlotStyle -> {Thickness[0.01]}, PlotLabel -> "Bifurcation Diagram"]

Since h(n) is maximized for n = 1, where \( h(1) = \frac{1}{2} , \) there is no bifurcation occurring at value n < 1. If we plot γ(c) versus c(n) over the allowed range of n, we obtain the phase diagram. The cusp occurs at \( \left( c^{\ast} , \gamma^{\ast} \right) , \) and is determined by the requirement that the two bifurcations coincide. This supplies a third condition, namely that the second derivative is zero: Hence, \( n = \sqrt{3} , \) which gives \( c^{\ast} = 3\sqrt{3} \quad\mbox{and} \quad \gamma^{\ast} = 3\sqrt{3} /8 . \) For c < c^*, here are no bifurcations at any value of γ.

 

Suppose that when the population of species falls below a certain level (called threshold, the species cannot sustain itself, but otherwise, the population follows logistic growth. To describe this situation mathematically, we introduce the differential equation where 0 < A < B. The equilibrium solutions for this differential equation are P = 0, P = A, and P = B. To understand the dynamics of the equation, we plot the phase portrait and pahe line:

StreamPlot[{1, 3*x*(1 - x)*(2 - x)}, {t, 0, 3}, {x, -0.4, 3}, FrameLabel -> {t, P}, StreamStyle -> Arrowheads[{{0.03, Automatic}}], MeshShading -> {Red, None}]

 

Separation of variables in the Gompertz equation yields

From the initial condition, we get the value of C to be where N₀ = N(0). Exponentiation leads to The same formula can be obtained upon substitution (from Thieme's book) The derived formulas exhitit a combination of measurable quantities--namely the tumor size N(t) and its asymptotic limiting size K--this yields a linear function of time if tumor growth. This observation plays an important role in fitting the Gompertz model to actual data. It is clear that N = K is a globally asymptotically stable fixed point, and the Gompertz model exhibits sigmoidal growth.

Next exponentiation yields the explicit formula for the Gompertz's curve:

Another re-parameterisation of the Gompertz model can also be written as the three parameter curve:

where Depending on circumstances, one may come across another version of the Gompertz model Example: Gompertz's COVID-19 models for Brazil

 

 

In 1939, the Czech scientist V. Korf proposed the following model to describe the dependence of a tree's height on its width. The Korf differential equation is easy to solve by separating variables. From the initial condition, we we get This yields

b = 2.3; NN[t_] = 0.2*Exp[0.2/b*(1 - 1/(1 + t)^b)] Plot[NN[t], {t, 0, 3}, PlotStyle -> Thick]

 

Another model was proposed by F.J. Richards in 1959: His model contains four common parameters: a generalised growth rate 𝑟 setting the typical time scale of the epidemic growth process, the final capacity K, which is the asymptotic total number of infections over the whole epidemics, the exponent 𝑎, which is introduced to measure the deviation from the symmetric S-shaped dynamics of the simple logistic curve, and one additional parameter 𝑎 ∈ [0,1].

Since the Richards equation is an example of a more general Bernoulli equation, we seek its solution in the form of the product R(t) = u(t) v(t), where u(t) is a solution of the truncated linear equation

Integration yields Substituting the product R(t) = u(t) v(t) into the Richards equation, we obtain a separable equation for the unknown function v(t): Integrating both sides, we get the general solution where C is a constant of integration. From the initial condition \( R_0 = R(0) = v(0) , \) we derive Multiplying u and v, we obtain the solution Its solution can be simplified to

P[t_] = Exp[2*t]/(-1 + Exp[2*t*0.3] + (1/0.2)^0.3)^(10/3) Plot[P[t], {t, 0, 8}, PlotRange -> {0, 1}, PlotStyle -> Thick]

Example: COVID-19 models for Sweden

 

The Weibull model is Example: COVID-19 models for Italy and South Africa

 

where b is a growth parameter for a fixed c and c is a shape parameter. The above equation can be written as For p < 1, this function exhibits exponential growth; and for p > 1, sigmodal growth occurs.

 

In 1949, the Austrian biologist Karl Ludwig von Bertalanffy (1901--1972) introduced the differential equation as a model for organism growth. He held positions at the University of Vienna (1934-48), the University of Ottawa (1950-54), the Mount Sinai Hospital (Los Angeles) (1955-58), the University of Alberta (1961-68), and the State University of New York (SUNY) (1969-72).

In 1934, he was habilitated by Reininger, Schlick and the zoologist Versluys for the first volume of his Theoretische Biologie. The monography postulated two essential aims of a theoretical biology, firstly to clean up the conceptual terminology of biology, and, secondly, to explain how the phenomena of life can spontaneously emerge from forces existing inside an organism. Here the organismic system represented the main problem as well as the still-to-be-formulated program of a theoretical biology. The second volume developed the research program of a dynamic morphology and applied the mathematical method to biological problems.

As a Rockefeller Fellow at the University of Chicago (1937-38), he worked with the Russian physicist Nicolaus Rashevsky. There he gave his first lecture about the General System Theory as a methodology that is valid for all sciences (1949). In 1939 he was appointed to an extraordinary professor at the University of Vienna. There Bertalanffy concentrated his research on a comparative physiology of growth. He was the first biologist who held lectures in zoology for students of medicine and an integrated course on botany and zoology. During this time he wrote an article about organisms as physical systems (1940), the summary of his biological reasoning: Problems of Life.

In 1949, Ludwig emigrated to Canada where he mainly worked on metabolism, growth, biophysics, and cancer cytology. In his biomedical research on cancer he developed, with his son Felix, the Bertalanffy-method of cancer cytodiagnosis. From the 1950's onwards he shifted his research from the biological sciences to the methodology of science, the General System Theory (GST), and cognitive psychology. Based on his humanistic world view, he developed a holistic epistemology (1966) which sharply criticized the machine metaphor of neobehaviorism (Robots).

In 1960 Ludwig was appointed a Professor of the Department of Zoology and Psychology at the University of Alberta in Edmonton (Canada). There Bertalanffy, the psychologist Royce and the philosopher Tenneysen established the Advanced Center for Theoretical Psychology that became a center for cognitive psychology over the next 30 years. In that time his systematic theoretical approach focused on the modern world of technology and has thrown us humans out of nature and has isolated us from each other. Bertalanffy emphasized in his later works the importance of the symbolic worlds of culture which we ourselves have created during evolution.

After his retirement, he became a Professor of the Faculty of Social Sciences at the State University of New York (SUNY). An international symposium celebrated his 70th birthday in 1971. In June 1972, he suffered a heart stroke and died a few days later, on June 12, shortly after midnight.

To sum up his life-work, Bertalanffy wrote 13 monographies, four anthologies, and over 200 articles, He was the chief editor of the Handbuch der Biologie--among many others. His themes encompassed theoretical biology and experimental physiology (Bertalanffy equations), theoretical psychology--particularly in the last two decades of his life, cancer research (Bertalanffy method of cancer cytodiagnosis), and philosophy and history of science. No doubt, the person Bertalanffy was a very fascinating one, proud of his European background, a connoisseur of architectural drawings, Japanese woodcuts, and stamps, who loved to hear the music of Mozart and Beethoven and to become absorbed in the works of Goethe. Sometimes he puzzled his environment by sarcastic remarks, or his black sense of humor. His worldview was an old-fashioned one, however, never outdated that was deeply rooted in a humanistic ethos:

In the last resort, however, it is always a system of values, of ideas, of ideologies - choose whatever word you like - that is decisive. (1964)

Bertalanffy (1957) derived his equation from the assumptions, which he attributed to Pütter (1920), that the rate of anabolism is proportional to the surface area of an organism (or to its mass raised to the power of ⅔), while catabolism is proportional to the organism's mass. Von Bertalanffy formulation (1957) of the equation is based on the biological principles of anabolism and catabolism, that is, growth and destruction. He derived his equation from the assumptions, which he attributed to Pütter (1920), that the rate of anabolism is proportional to the surface area of an organism (or to is mass raised to the power of ⅔), while catabolism is proportional to the organism's mass. The Bertalanffy model is often used as a population model to study the factors that control and affect its growth. It is used to describe the populations of fish, animals, such as pigs, horses, cattle, sheep, etc. It is also suitable for modelling tumors and infectious diseases because it is similar to the growth of individuals and populations:

Since the Bertalanffy equation is a Bernoully equation, we solve it by representing the solution in the form of the product of two functions: P(t) = u(t) v(t). Here u(t) is a solution of the truncated linear part:

Its solution becomes because we need just one such solution but not the general one. Substitution of the product into the Bertalanffy equation, we obtain for v(t) the following separable differential equation: Integration yields where C is a constant of integration. Multiplication of these functions gives the required general solution: From the initial condition, it follows that Therefore, we finally get the solution for the von Bertalanffy model: where t is time, P(0) is the population size at time t = 0, K is carrying capacity, r is the intrinsic growth rate.

This model attains its inflection value 8 K/27 at

Plot[x^(2/3) - x, {x, 0, 1}, PlotStyle -> Thick] P[t_] = Exp[-2*t/1]*(1*Exp[2*t/3/1] - 1 + 0.1^(1/3))^3 Plot[P[t], {t, 0, 8}, PlotRange -> {0, 1}, PlotStyle -> Thick]

To capture the dynamics of harvesting a population or tumour treatment we can add a harvesting term E(t) P(t). Here, E is also considered to vary with time and represents harvesting effort. Adding this term to the Bertalanffy population model results in a new harvesting IVP,

Here 𝑎(t), b(t), and E(t) are assumed to be known positive functions of time.

U-Bertalaffy:

The von Bertalanffy (1957) growth law can be written as

where η and λ are constants. Once this differential equation is integrated, the von Bertalanffy growth curve becomes where is the limit to growth parameter, β = α/η (α is a constant of integration), λ is a growth rate parameter, 1-k is a shape parameter, and the initial value is If k > 1 and both η and λ are negative, the growth curve is sigmodal. von Bertalanffy determined empirically that k = 2/3 for a wide variety animals.

 

The spruce budworm is an insect which is a serious threat to forests. In 1978, Ludwig, Jones, and Holling published a paper in which they proposed an ingenious model of the interaction of the insects, the trees, and the predatory birds. They proposed that the budworm population P obeys a logistic model modified by a predation rate:

with some positive constants r, α, k, and β. The parameter β represents the predation efficiency of birds and it is proportional to the population of birds, considered a constant in this model. The parameter k is proportional to the average leaf surface area (about half of it), and it is a scaled measure of the average foliage density. So trees with larger leaves attract more budworms. The denominator, \( k^2 + P^2 , \) in the last term of the equation models the laziness of birds that go for places where food density is high allowing them to consume much while expending minimal energy.

The analysis of the above equation for P(t) can be simplified by introducing dimensionless variables \( t=a\tau \) and P = by, where \( a= k/\beta \) and b = k. Then the equation can be written in the following form:

where \( R=rk/\beta \) and \( A=\alpha k. \) The above equation contains two parameters. Its slope function has obviously one null: y = 0, but other critical points are obtained by solving the equation The first step in the analysis is to identify the equilibria (where the derivative of y is zero) and determine their stability. Since the slope function is positive in a neighborhood of y = 0, this equilibrium point is unstable. Upon multiplying both sides of the latter equation by 1 + y², we get the cubic equation that may have either one real root or three real roots, depending on the values of parameters A and R. The best way to do the analysis is to examine what can occur for each value of R as A varies.

The interpretation of the above equation is both simple and important. Its left-hand side is per capita growth rate of the scaled variable y (with respect to a scaled time). The right-hand side is per capita death rate due to predation, also in scaled variables. The points where these curves intersect are nonzero equilibria for y.

R = 0.5; A = 0.1; a = Plot[{R*(1 - A*y), y/(1 + y*y)}, {y, 0, 10}, PlotStyle -> {Thick}, Epilog -> Text[Style["unstable", 10, FontFamily -> Times], {2.2, -0.1}], PlotRange -> {-0.2, 0.5}] b = ListPlot[{{0.7, 0.465}, {2.05, 0.4}, {7.34, 0.132}}, PlotStyle -> {Large, Black}] line1 = Line[{{0.7, 0.46}, {0.7, 0}}] line2 = Line[{{2.05, 0.39}, {2.05, 0}}] line3 = Line[{{7.34, 0.13}, {7.34, 0}}] c = {Graphics[{Dashed, Red, line1}], Graphics[{Dashed, Red, line2}], Graphics[{Dashed, Red, line3}]} Show[a, b, c]

The critical values of R and A are where two roots coalesce and disappear. The bifurcation points are obtained from the following system of equations:

 

References