Some Basic Differential Models in Mathematical Biology

It is widely known that many problems in the mathematical biology can be modeled by the differential equations, and the corresponding systems are usually called as biological dynamic systems. In this review, we will introduce some basic biological dynamic models which are widely applied in various fields, such as population ecology models, epidemic models, virus dynamic models, pharmacokinetic models, chemostat models and so on. It is expected that these models could provide some references and enlightenments for the researchers in the related fields, then some more realistic and generalized models could be established based on them.

Keywords: Differential equations; Mathematical biology; Population ecology; Epidemics; Virus dynamics; Pharmacokinetics; Chemostat models

Mathematical biology is a multidisciplinary science of life science, biology, agriculture, medicine and public health and other disciplines. It not only studies and solves some problems in life science, but also studies on the relevant mathematical methods and theories.

Since the late 1940s, with the invention and application of the computers, mathematical biology entered a new period of rapid development. And many branches of mathematical biology gradually emerge, such as numerical taxonomy, bio-cybernetics, biodynamics, bioinformatics etc. By the middle 1970s, with the help of the new theories and methods of differential equations, dynamical systems, as well as the greater power of the mathematical software and the large databases, mathematical modeling by differential equations and corresponding dynamical theories are successfully applied in population ecology, epidemiology, neurobiology, viral kinetics, pharmacokinetics, physiology, environmental pollution problems and other fields, and many phenomena have been well explained in these fields.

Moreover, in the recent years, mathematical biology have played an very important role in the study of the growth of the cancer cells and bacteria, plasma drug concentration, the periodic regulation of cells, gene regulation, molecular biology, medicine and so on.

Since the authors are always doing some rough work in some of these areas, such as population ecology, epidemic dynamic models, chemostat models in the last decade. And it happens to the time that the international journal «Journal of Bioequivalence Studies» will publish the inaugural issue, the authors collect some basic differential models in some areas of the mathematical biology according to limited knowledge and understanding. Moreover, we have referenced some books published in the last decade [1,2]. The aim of this article is to provide a reference and enlightenment for the researchers in the related fields so that they could establish more generalized models and make better work in future.

Any specie in the natural world does not exist solely, but is closely related to other species in biological communities and then constitute a population ecosystem. However, a single population is the basic unit which composed of the entire ecosystem. In order to predict the change law of the population, the famous demographer Malthus proposed the following basic model for single specie,

Where N (t) denotes the population density at the time t, r denotes the intrinsic growth rate, which is the difference between the birth rate and the death rate.

Based on model (1), the famous ecologist Logistic proposed the famous insect population model when studying the growth rate of the insect in the laboratory as follows,

Where r denotes the intrinsic growth rate, K is the carrying capacity of the environment.

This model revised the deficiency of Malthus' model (1) which does not consider the limited resources of the environment, and described the population growth law has a characteristic of S-shaped, and this conforms to the actual ecology much better.

As far as the dynamic models of the ecosystem with multiple species are concerned, it should date back to the famous Lotka- Volterra model of two species, which is proposed by American population scientist Lotka when studying chemical reaction in 1921 and Italian mathematicians when considering the fish competitive system in 1923. The model is,

Where N_i(t) denotes the density of the i-th population at time , a_i>0 denotes the inner competitive coefficient of the i-th population, which is also called density dependence coefficient. If b_i<0, system (3) is a competitive system, and if b_i>0 system (3) is a cooperative system, while if b₁b₂<0, system (3) is a predator-prey system, i = 1,2.

More complicated dynamic models, the interaction of three species for example, can be seen in the reference [1,3]. For these population dynamic models, the global stability of the equilibriums of the system, extinction and persistence for the system, existence and global stability of periodic solution or almost periodic solution, bifurcations, chaos and other properties are usually discussed. And by studying of dynamical behavior of corresponding models, we can provide valuable information and advice to control pest outbreaks, to keep variety of some precious specie, to maintain ecological balance and so on.

As the epidemic models are concerned, it should date back to the famous SIR model, which was proposed by Kermack and Mckendrick in 1927 when studying the propagation law of the Black Death in London from 1655-1666 and the plague in Mumbai in 1906 [4,5]. And the famous SIR model is as follow,

Where S(t) denotes the number of the susceptible individuals which has not yet infected but may be infected by the bacteria at time t, I(t) denotes the number of the infected individuals at time t, and R(t) denotes the number of the removed individuals at time t, β, denotes the infected rate, γ denotes the cure rate and γ^-1 denotes the average cure rate.

A basic assumption of above model is that the birth and the death was not considered in a closed environment and there is no inflow and outflow in the system, then the total number of the individuals keeps constant. Further, it is supposed that the infected individuals have permanent immune capacity after the cure.

If the individuals suffer from a certain incurable epidemic, then the corresponding model degenerates to the following SI model,

Where d is the death rate due to illness.
If the individuals suffer from a certain curable epidemic, but have no immune capacity after cure, then the corresponding model should be following SIS model,

If the individuals suffer from a certain curable epidemic and have a temporary immune capacity after cure, there are a number of δR(t) removed will become infected again per unit time, then the corresponding model can be described by following SIRS model,

Where δ^-1 is usually called as average immune period.
In addition, if there is an incubation period before a susceptible individual becomes an infected one and suppose that the infected individual in the incubation period has no infectivity. E(t) denotes the number of the infected in the incubation period at time t, and the infected individuals in the incubation period will transform to the real infected ones by a proportion coefficient ω, that is, ω^-1 is the average latency.

For above case, if the infected individuals have permanent immune capacity after the cure, then the corresponding model becomes following SEIR model,

In the same way, if the individuals suffer from a certain curable epidemic and have a temporary immune capacity after cure, then the corresponding model can be described by following SIERS model,

Where the meaning of the coefficients β, λ, ω, δ is the same as above.
By some qualitative and quantitative analysis, such as global stability of the disease free equilibrium and the endemic equilibrium, as well as numerical simulations of the epidemic models, we can reveal the propagation regularity, predict the variation trend, analyze the major cause and the key factors (e.g., the basic reproduction number) of an influenza epidemics. And these can help us understand the population rules of the disease more comprehensively, then help us find an optimal strategy to predict and control all kinds of diseases.

Research shows that the infection of Human immunodeficiency virus (HIV), hepatitis B virus (HBV) and Hepatitis B virus (HCV) is a dynamic process of viral replication and clearance. A large number of virus produce and disappear everyday, the emergence of antiviral drugs promotes the mathematical analysis of the virus dynamics. By the periodic examination the viral load in the blood after administration, we can establish suitable mathematical models and make some numerical simulations and then understand the dynamic process of the infection, the replication and the clearance of the virus in the body.

In this section, we will introduce several classic dynamic models for the HIV.
In 1997, Bonhoeffer and May [6] proposed a simple mathematical model of HIV and CD4+T cells, in which they divided CD4+T cells into two parts: healthy and infected CD4+T cells and established following dynamic model,

Where T(t), T*(t) and V(t) denotes the concentration of the healthy, infected CD4+T cells and the free virus at time , respectively. λ, d, β, N, δ, c are positive constants. λ is the velocity of the CD4+T cells produced in the body, d is the death rate of the healthy CD4+T cells, is the burst size of each infected viral particle after death, c is the death rate of the free cells, β is the infected rate of the CD4+T cells.

For the infected process of the HIV, there are two kinds of anti-HIV drugs, one is the reverse transcriptase inhibitor, the other is the protease inhibitor. The function of the reverse transcriptase inhibitor is to inhibit the activity of the reverse transcriptase and prevent the production of RNA and DNA, and prevent the infection of the new CD4+T cells.

Suppose that the efficiency of the reverse transcriptase inhibitor is η_RA, 0 ≤η_RA ≤ 1 and η_RA=1 means that the efficacy of the reverse transcriptase inhibitor is 100%, then the corresponding dynamic model is as follow [7,8],

Suppose that the efficiency of protease inhibitor is η_PI, 0≤ η_PI≤1,and η_PI=1 means that the efficacy of the protease inhibitor is 100%, and the corresponding dynamic model is,

Where V₁ and V_NI denotes the concentration of the infectious and the non-infectious virus particles, respectively. In addition, if both the reverse transcriptase inhibitor and the protease inhibitor is used, and this kind of treatment is often called as highly active antiretroviral therapy (HAART), and the corresponding dynamic model is,

As the research on the virus dynamic models is concerned, the primary concern is the global stability of the virus free equilibrium, the reproduction number of the virus and so on. Particularly, the dynamic changes of the number of the viral load in the body can be analyzed and observed intuitively. Thus, it can be applied to explore the effectiveness and drug resistance of certain treatment, and then some effective strategies for the treatment strategy can be provided for the virus.

Pharmacokinetics is a science to study the variation law of the pharmaceutical dosage or concentration by applying dynamic theory and mathematical methods. Moreover, it studies the variation law of the processes of the absorption, distribution, biotransformation and the excretion of the drug when it enter the body by all kinds of administration methods, then describes these processes by some mathematical models [2,9].

The plasma drug concentration is usually used in the practical applications since the concentration of the action site of the drug is difficult to measure. Therefore, the plasma drug concentration is often used as the variable in the mathematical modeling. Moreover, there are two key factors for the compartment model of the pharmacokinetics, one is administration routes, the other is the elimination process of the drug in the body. There are three kinds of common used administration routes, i.e., oral or extravascular administration, intravenous injection and intravenous drip with constant rate.

Drugs enter the body by different routes, and lead to the change of the amount in different positions and at different times, which involves the elimination process of the drug in the body. In the following, we will introduce three pharmacokinetics models for above three different administration routes.

For most drugs, it is necessary to dose repeatedly so that the plasma drug concentration can reach and keep in the efficient range, and then achieve the desired efficacy. Thus, the periodically repeated intravenous administration is very important for the legitimate use of drugs and preparations. As multiple dosing by a regular time interval and equal amounts repeatedly, the baseline of the plasma drug concentration will change every time. Finally, it will reach a steady state as the drugs accumulate in the body gradually. At the moment, it can be described by the following impulsive differential equation,

Where C(t) denotes the plasma drug concentration at time t, D denotes the drug dosage of the intravenous injection each time, V denotes the apparent volume of distribution, K is the elimination rate constant.

Intravenous drip with constant rate is an administration method by continuous medication into the vascular with a constant rate. If the drug with dosage D is intravenously dripped into the body with constant rate r=D/TV in time T, then the velocity of the drug into the body is independent on the concentration. For this case, the corresponding dynamic model is,

It is easy to solve above initial values problems (15) as,

Making a simple computation that which is usually called as steady-state concentration or plateau concentration.

Extra-vascular administration includes oral, intramuscular or subcutaneous injection, etc. For a single injection, since there is an absorption process for the drugs entering from the medication position to the blood circulation. Thus, the drug can't enter the blood circulation immediately as the intravenous injection. According to this feature, we only need to append an additional absorption compartment before the intravenous injection model, and the drug will enter the absorption compartment at first, then enter the center compartment gradually. That is to say, drug with dosage D will be absorbed after entering the absorption compartment and transmitting it to the center room, then it is easy to obtain the single-compartment model of extra-vascular administration. If C_α denotes the concentration of the drug in the absorption compartment, V_α denotes the apparent volume of distribution, then the variation law of the drug in the absorption compartment satisfy the following dynamic model,

Where K_α>0 is the elimination rate constant in the absorption room.
It is easy to solve initial value problem (17) and obtain C (t). Assume the drug in the absorption will transmit into the center room totally, and K denotes the elimination rate constant in the center room, then the plasma drug concentration C(t) in the center room can be described by following dynamic model,

And it can be also solved as,

Which reflects the relationship of the plasma drug concentration in the body to the time t.

In the same method, we can establish the impulsive differential model of the periodic extra-vascular administration, and we omit here, more details can also be seen in [2,9].

In addition, according to the pharmacokinetic properties in vivo, the compartmental model can be divided into single-compartment model, two-compartment model and multi- compartment model. For a single-compartment model, the transport rate of the drugs in all of the organ or tissue is nearly the same or similar, and the entire body can be regarded as a compartment. But for many drugs, the transport rate between organ and organ, organ and tissue, muscle and bone are usually different. Thus, it can not reflect the difference of the plasma drug concentration in the different parts of the body by the pharmacokinetic analysis of the singlecompartment model. Instead, it is more reasonable to describe the variation law of the plasma drug concentration, the relationship between dosage and concentration or time by the multi-compartment model.

In a multi-compartment model, it is supposed that the organism is divided into a central compartment and several compartments. For example, three-compartment model is divided by the central compartment and two compartments. From the viewpoint of mathematical modeling, modeling methods of multi-compartment model are the same as two-compartment model. Thus, we only introduce the two-compartment model in this review.

For a two-compartment model, the organism is divided into a central compartment (whose apparent volume is smaller) and a peripheral compartment (whose apparent volume is bigger). And the central compartment is composed of some organizations, such as heart, liver, lung, kidney, where the blood is richer, the membrane permeability is better and the transport is easier. Drugs often enter this compartment at first. While the peripheral compartment is composed of some organizations, such as fact and muscles, where the blood is not rich and the transport rate is slower, and it will take a period of time to achieve the dynamic balance of the plasma drug concentration. Therefore, the variation of the dosage or the concentration, the time to achieve the dynamic balance in the central and the peripheral compartment will directly affect the treatment effect. By this means, the modeling process of the extra-vascular administration is similar as this two-compartment model.

As mentioned above, for a two-compartment model, drug will firstly enter the central compartment and reach a dynamical balance quickly after intravenous injection, then eliminate by the process of the first order rate in the central compartment and transport to the peripheral compartment. Moreover, the drug transport between the central and the peripheral compartment is reversible. So it can be described by the dynamical process of the first order. If we denote C₁(t),C₂(t) the plasma drug concentration in the central and the peripheral compartment at time t, and D₁(t),D₂(t) the dosage of the drug in the central and the peripheral compartment at time t, respectively. Then the variation of the dosage after an intravenous injection can be described by following differential equations,

Where k₁₂,k₂₁,k_e is positive constant.
Further, according to the relationship between the dosage and the concentration, we can obtain following models of the dosage,

Where K₁₂ = k_{12V1^-1, K21 = k21V^-12} denotes the first order rate transport constant of the drug from the central compartment to the peripheral compartment and the opposite direction. K_e = k_eV^-11denotes the first order rate elimination constant in the central compartment. V₁,V₂ denotes the apparent volume of distribution. Since the equations (21) are linear homogeneous differential equations with constant coefficients, then it is not difficult to solve it and obtain the solutions, which describe the relationship between the dosage and the time.

In the study of pharmacokinetic models, the analytic solutions of the compartment models in different kinds of administration routes lead a very important role in the pharmacokinetics

And in the practical application, the analytical formulas can be directly used for the parameters estimation of the drug and the model determination, then provide necessary theoretical basis for the choice of the administration routes. However, if the transport and the elimination process which the drug from the central compartment to the peripheral compartment satisfies the Michaelis- Menten's relationship, then the compartment model is not a linear model any longer. And for this kind of nonlinear compartment models, it is impossible to obtain the analytic solutions. Thus, more complex theories on ordinary differential equations and even impulsive differential equations are necessary for the study on these problems, so one can make some further study in this field and maybe it will promote further development of pharmacokinetics.

The chemostat is an experimental device used for continuous cultivation of microorganism, and it is composed of three connected containers. The first container provides nutrients sufficient enough for the microbial growth, whose main nutrition component is called the nutrient medium. It is limited, and the concentration of the nutrient remains constant in the first container. Then, the nutrient medium is pumped to a second container (which is the so-called culture container) at a constant rate. At the same rate, the material (the mixture of the nutrient and the microorganisms) in the container is pumped into the third container to keep its capacity unchanged.
The classic chemostat model can be described by the following differential equations [10,11],

Where s(t) denotes the concentration of the nutrient at time t and x(t) denotes the concentration of the microorganism at time t, s₀ denotes the initial concentration of the nutrient in the culture container, D is the dilution rate, the stoichiometric coefficient δ denotes the proportion of the amount of the produced microorganism to the consumed nutrient, p(s) denotes the function of the growth rate and one of the most popular functions is a so-called Monod function [12] with the following form,

Where μ_m is the maximal growth rate, K_m is the half-saturation coefficient.
For the case of multiple microorganisms competing nutrient in a chemostat, Taylor and Williams [13] establish following dynamic model,

Where s(t) denotes the concentration of the nutrient at time t and x_it denotes the concentration of the i-th microorganism at time t, y_i,μ_i,K_i denotes the growth constant, the maximal growth rate, the half-saturation coefficient of the i-th microorgnism, respectively. D_i = D + ε_i, D is the dilution rate, iε is the special death rate of the i-th microorganism.

If multiple chemostat are linked together, it is often called as gradostat. The corresponding model and dynamical behavior is more complicated. References for the interaction between multiple chemostat can be seen in the literature [14].

Chemostat model can be used to describe the microbial growth process in the laboratory, in the natural rivers, lakes, and so on, and it has very important significance in the field of microbiology and biochemistry. Therefore, it has been long-term concerned by the biological and mathematical researchers.

Besides models in the field mentioned above, there are many differential models in the mathematical biology, such as biochemical reaction model [15], neural dynamics model [16], cell cycle regulation model [17], the gene regulation model [18] and so on. However, we are not going to list all of these in the review by the constraints of knowledge and time.

Generally speaking, in this review, we introduce the basic method of modeling by the differential equations and provided several basic models for the readers. All of these models can be applied in the corresponding fields, and it can be used to explain some phenomena effectively in the related fields, provide a theoretical basis to some control strategies and so on. However, considering the theme of the journal, the introduction of the pharmacokinetic model is longer than the other models, and the main purpose of this review is providing a communication platform for the ecological and biological researchers, especially the medical researchers with the mathematical researchers.

This work is supported by the National Natural Science Foundation of China (51349011) and Scientific Research Fund of Sichuan Provincial Education Department (11ZB192, 14ZB0115).