September 24, 2017, Sunday
University of Colorado at Boulder Search A to Z Campus Map University of Colorado at Boulder CU 
Search Links

MBW:Computer Modeling of the Biotic Cycle Formation in a Closed Ecological System

From MathBio

Jump to: navigation, search


This article is a summary, analysis, and expansion of [1]. [1] All information is taken from this source unless otherwise noted.

Executive Summary

The authors of this paper sought to illustrate biotic turnover in a closed ecological system (CES). A producer-consumer-reducer system was modeled using a set of differential equations, each relating to the formation and development of the six species in the system. At the same time the authors tracked the energy flow into the system by tracking the primary production. The term primary production refers to the conversion of carbon dioxide into organic compounds. In this case photosynthesis is the mechanism by which the system creates these compounds which are then the energy source for other species. The system tends to increase its energy uptake as layers of complexity are introduced. When a predator is placed into the system at steady state, or competition is introduced in the primary consumer level, the biotic turnover increases. The authors concluded that for every ecological system it is possible to find a suitable predator which will increase the energy flow into the system.

History and Context

This paper examines a direct application of the laws of thermodynamics to an ecological (specifically a biotic cycle) system. When applying the laws of thermodynamics to the system, there is a misconception that nature contradicts what is defined in the laws. Specifically the second law of thermodynamics seems violated by living systems. These living systems are seemingly very structured which implies a lower entropic state, because entropy is in many ways a measure of disorder. This understanding is not correct. If you were to measure the disorder you need to figure out the total microstates that the system could be in for its given macrostate. Though the system appears to evolve to be more ordered, one must look elsewhere to decide what the entropy is doing. The entropy of a system is also proportional to the energy state of that system, so we track the primary production in order to track the energy being used by the system as consumers or other structures are introduced. The system must burn through more energy in these more "ordered" states because, quite simply, there is more life to support and life needs energy to survive. It should be noted, however, that the physical criteria that is linked to evolution does not correspond to the theoretical biotic cycle that is posed by the authors. According to a paper by Pechurkin, most living systems have a higher use of their free energy as well as a higher dissipation of that same energy. [2]

The authors decided to add complexity to a simple producer-reducer system by adding primary and secondary consumers while tracking the flow of energy into the system. Their hypothesis was that the rate of energy used by the system would increase as the complexity increased.

Monod/Michaelis-Menten Kinetics

The Monod equation is an essential a mathematical model used to represent and structure growth of microorganisms. [3] It is given as:

                        \mu =\mu _{{max}}{\frac  {S}{{K_{S}}+S}}

where the parameters are represented as:

Species Description
\mu Growth Rate
{\mu _{{max}}} Max Growth Rate
S The Concentration of the Limiting Substrate for Growth
{K_{S}} Half Velocity (i.e. when {\frac  {\mu }{\mu _{{max}}}}=0.5)

Much of Monod's equation is based on basic Michaelis-Menten kinetics used to describe the growth rate of certain microorganisms. Monod's derived equation uses the same basic approach that Michaelis-Menten kinetics uses, however Monod built his model through observations he conducted. Michaelis-Menten kinetics are based on theoretical ideas. Monod's equation is extremely useful in certain disciplines in this day and age and is specifically used in the "Activated Sludge Model". This model is used in treatment facilities using biological floc of bacteria. Growth of those certain bacteria is modeled using Monod. The applications of Monod are very important and play an important part in understanding our closed ecological system.

Mathematical Model

The model that the authors use contains several species:

Species of the Model

Species Description
P Producers (mg/L)
{C_{1}} Consumer of First Order(mg/L)
{C_{2}} Consumer of Second Order(mg/L)
D Detritus(mg/L)
R Reducers(mg/L)
N Nitrogen (essential life nutrients) (mg/(L*hour))

Examples of these parameters in the biotic system:
Producers: autotrophs, green algae
Consumers of the first order: water flea, small crustaceans
Consumers of the second order: guppy fishes, small fresh water fish
Reducers: microorganisms
The other two parameters are D (detritus) which essentially marks death (non living organic matter) and N (Nitrogen) which also marks the concentration of essential life nutrients which is generally categorized as nitrogen.

Figure 1: The biotic cycle

Each parameter is modeled with its own differential equation:

                        {\frac  {dP}{dt}}={\mu _{P}}P-{\epsilon _{P}}P-{\frac  {{\mu _{{C_{1}}}}{C_{1}}}{Y_{{C_{1}}}}}

This differential equation models the producer, specifically the birth and death rate of the producer and consumption of the producer by C1. $ Y_{C_1}$ is the yield coefficient of consumer 1 with respect to the producer.

                        {\frac  {d{C_{1}}}{dt}}={\mu _{{C_{1}}}}{C_{1}}-{\epsilon _{{C_{1}}}}{C_{1}}-{\frac  {{\mu _{{C_{2}}}}{C_{2}}}{Y_{{C_{2}}}}}

This differential equation models consumer 1, specifically the birth and death rate of the consumer 1 and consumption of consumer 1 by C2. $ Y_{C_2}$ is the yield coefficient of consumer 2 with respect to the consumer 1.

                        {\frac  {d{C_{2}}}{dt}}={\mu _{{C_{2}}}}{C_{2}}-{\epsilon _{{C_{2}}}}{C_{2}}

This differential equation models consumer 2 and portrays the growth and death rate of the population. Notice that this is the highest order consumer so there is no consumption term.

                        {\frac  {dD}{dt}}={\epsilon _{P}}P+{\epsilon _{{C_{1}}}}{C_{1}}+{\epsilon _{{C_{2}}}}{C_{2}}+{\epsilon _{R}}R-{\frac  {{\mu _{R}}R}{Y_{R}}}

This differential equation encapsulates all of the death rates of all the population and accounts for the reducer population's consumption of detritus. The reducer is what helps break down the already dead population and converts it to organic nutrients (nitrogen) for the producer population.

                        {\frac  {dR}{dt}}={\mu _{R}}R-{\epsilon _{R}}R

This equation describes the aforementioned reducer species and simply describes the growth and death rate of the population.

                        {\frac  {dN}{dt}}={k_{R}}{\frac  {{\mu _{R}}R}{Y_{R}}}+{k_{{C_{1}}}}{\frac  {{\mu _{{C_{1}}}}{C_{1}}}{Y_{{C_{1}}}}}+{k_{{C_{2}}}}{\frac  {{\mu _{{C_{2}}}}{C_{2}}}{Y_{{C_{2}}}}}-{\frac  {{\mu _{P}}P}{Y_{P}}}

This equation is the nitrogen production per time. Each coefficient k is in units of nitrogen/food consumed (in the case of consumer 1 it would be nitrogen/producer).

{\mu } relates to specific growth rates within producers, consumers (of both first and second order) and reducers. {\epsilon } corresponds to death rates within the same groups of species. Y are yield coefficients for the producer, consumers and reducer populations.

Initial Results

Figure 2: where 1 - producers, 2-consumer, 3-nutrients,4-detritus (all in mg/liter), and 5-primary production (in mg/(liter*hour))

Complexity of the Biotic Cycle

The authors were able to find an increase in the primary production beginning from the producer, detritus, reducer system and adding the consumers of the first and second order. In the graph at the right, you can see the increase of the primary production after the introduction of the consumers, until it asymptotically reaches a higher steady state. The authors found that there was a higher impact from first order consumers over second order consumer within the whole cycle.

Competition between Producers

When there is competition between certain producers of a biotic cycle there is an increase of primary production.

Competition between Consumers (of order 1)

There is an increase of primary producers if the introduced consumer is not effective in consuming its prey, but if the consumer has a higher growth rate and yield coefficient than the competitor, the flow of energy can potentially decrease after new consumers are introduced.

Formation of Biotic cycle

Throughout the experimentation the authors were able to look into co-evolution. Specifically they found that with cooperative co-evolution there tends to be an increase of energy flow in the form primary producers. Essentially the authors were able to correlate the idea that as consumers tend to increase there is an increase in energy from producers due to the need to support a stable biotic cycle.

Adjustments to Model

As we further explored the authors model we went ahead and made certain adjustments to the model. Within the given model the parameters were undefined. So we used estimates and then used a tool (described below in the analysis) to optimize our parameter values. The ODE's listed above did not include the Michaelis-Menten terms and we added those parameters in for growth rates.

Initial Analysis

No Consumers

We began by running the simulation with no consumers and looking at concentrations of each species:


This steady state is a product of our initial conditions; more nitrogen to begin could lead to an explosion of producer population, but as it is the producer population grows until the amount of nitrogen is too low, and then it dies off. Also, if we had our reducer convert Detritus to Nitrogen more rapidly, the producer population could be steady or explode. As it is, this is the steady state of the system.

One Consumer

We continue and add in a single consumer:


We see that now the population of producers does not go to zero. Instead we see a steady population for a short time and then exponential growth. In real life the system populations would rise to the carrying capacity of the system (which is not described in our model).

Two Consumers

Next, we introduced the second consumer to see how the system progresses: Thurimella2Consumer.jpg

It is important to note that there is not much change in the behavior (when compared to the introduction of consumer 1) but there is a higher overall production within the system. Both consumers are very small but are non-zero and are assumed to reach a carrying capacity within the system and be able to sustain itself. Also note that {\frac  {dP}{dt}} is choppy because our model is doing an analysis at discrete steps (one point per day).


Finally, we introduced varying levels of competition between each of the two consumers, changing the food for each one and comparing how the system progresses in each case:

When there is no competition:
Shows our system starting from a steady state predicted in mathematica.
When the second order consumer's role is changed so that it eats producers and consumers we find the following: ThurimellaCompetition.jpg
The notion of competition changes because now this describes when second order consumers can now eat first order consumers and producers which the populations growing exponentially sooner.
When there is direct competition (i.e. both order consumers, consumer producers): ThurimellaDirectCompetition.jpg
Now we see that the system explodes even sooner, in a real closed ecosystem it would reach some carrying capacity.

These results show that higher competition leads to much larger numbers of producers and agrees with the prediction made by the authors.

Parameter Optimization

We used the following algorithm to optimize our unknown parameters. Initially we used the known numbers that we had from the paper, and then made estimates for the others based on the known parameters. Then using mathematic we entered in our system and found the Jacobian which simply takes the partial derivative of the right hand side of each equation with respect to each species, and puts it into a 6x6 matrix.

An example of a jacobian in a 3 species system where the x1, x2, x3 are the right hand sides of the ode's:

J_{F}(r,\theta ,\phi )={\begin{bmatrix}{\dfrac  {\partial x_{1}}{\partial r}}&{\dfrac  {\partial x_{1}}{\partial \theta }}&{\dfrac  {\partial x_{1}}{\partial \phi }}\\[3pt]{\dfrac  {\partial x_{2}}{\partial r}}&{\dfrac  {\partial x_{2}}{\partial \theta }}&{\dfrac  {\partial x_{2}}{\partial \phi }}\\[3pt]{\dfrac  {\partial x_{3}}{\partial r}}&{\dfrac  {\partial x_{3}}{\partial \theta }}&{\dfrac  {\partial x_{3}}{\partial \phi }}\\\end{bmatrix}}

At this point our jacobian has no parameters defined, so we plug in our estimated values and leave one parameter as a free variable.

Using the eigenvalue equation det(A-I*\lambda )=0 we find the 6 eigenvalues of the system as a function of our free parameter. By examining the roots of the functions we can find where the eigenvalues are negative, and therefore have a stable equilibrium.

Below we see 6 graphs which are plots of the 6 eigenvalues as functions of \mu _{P}. In the first 4 cases we see that the eigenvalues are negative throughout the region, on the 5th we have a zero eigen value and on the 6th we see that below about \mu <0.9 the eigenvalue is positive. So in this case we picked a value of 1.5 to ensure a stable steady state. We then continued on with this process until all our unknown parameters were optimized.

ThurimellaEigenvalue1.jpg ThurimellaEigenvalue2.jpg

We were able to reproduce the results in the paper more directly. After fixing our code so that we could calculate the primary production properly, we were able to produce the following results. In the first graph you can see the blue line for primary production, and the red line for the nitrogen concentration. What we see is that the nitrogen runs out and the primary production drops to a steady state around 1.5. By introducing the consumers into the system we arrive at the second graph where we see that the primary production increases to a steady state at about 2.8. Clearly the introduction of the consumers has increased the energy flow and biotic turnover of the system.



In general, we do in fact see that introducing more consumers increases the overall energy production as predicted by Bril'kov. As shown above when introducing the first consumer increases the number of producers by quite a large number. Also when we introduce in a second consumer we end up getting a similar result where the number of produces increases by quite a bit. Our analysis on competition between the consumers clearly illustrates how competition fuels an increase in the number of producers. However, the conclusions we were able to draw were very general and qualitative. So, in order to get more quantitative results, we performed further analysis on the model as described below.


The original model proposed in the article by Brilkov et al. successfully shows how the energy within a biotic cycle changes when additional predators are introduced. However, for the most part their conclusion is qualitative, as the paper is lacking in specific numbers corresponding to the parameters, and the concentrations of the different populations are unclear. We have expanded upon their model, performing analysis looking at competition between the consumers, the stability of the system with respect to each variable, and attempted to optimize the parameters. Ultimately, we were able to confirm their results by observing the expected changes in the biotic cycle as the numbers of consumers, and the competition between them, changed. Although we cannot speak to the quantitative accuracy of our analysis, as we our model is in terms of numbers of species rather than biotic mass concentrations, the successfully adaptation of the Monod equation to Brilkov et al's model shows the expected behavior of a closed ecological system.

Our Presentation

Our Presentation (Thurimella, Hass and Reisman) in pdf format.


  1. A.V. Brilkov, V.V. Ganusov, E.V. Morozova, and N.S. Pechurkin. "Computer Modeling of the Biotic Cycle Formation In A Closed Ecological System," Adv. Space Res., 27(9): 1587-1592, 2001.
  2. Shirobokova IM and N.S. Pechurkin. "Mathematical modeling of response of ecosystems with different structure to external impact.," Adv. Space Res., 27(9): 1593-1598, 2001.
  3. J. Monod. "The Growth of Bacteria Cultures.," Annu. Rev. Microbiol., 1949.3:371-394.