Journal of Applied Mathematics and Physics
Vol.07 No.01(2019), Article ID:90023,18 pages
10.4236/jamp.2019.71013

Application of Optimal Control to the Epidemiology of Dengue Fever Transmission

Okey Oseloka Onyejekwe*, Ayalnesh Tigabie, Biruk Ambachew, Abebe Alemu

Computational Mechanics and Dynamical Systems Group, Computational Science Program, Addis Ababa University, Arat Kilo Campus, Addis Ababa, Ethiopia

Copyright © 2019 by author(s) and Scientific Research Publishing Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

Received: December 3, 2018; Accepted: January 18, 2019; Published: January 21, 2019

ABSTRACT

In this paper, we build an epidemiological model to investigate the dynamics of the spread of dengue fever in human population. We apply optimal control theory via the Pontryagins Minimum Principle together with the Runge-Kutta solution technique to a “simple” SEIRS disease model. Controls representing education and drug therapy treatment are incorporated to reduce the latently infected and actively infected individual populations. The overall thrust is the minimization of the spread of the disease in a population by adopting an optimization technique as a guideline.

Keywords:

Epidemiological Model, SEIRS, Dengue Fever, Optimal Control, Pontryagins Minimum Principle, Fourth-Order Runge-Kutta

1. Introduction

Dengue fever is a painful, debilitating mosquito-borne disease caused by one of four closely related dengue viruses (Noorani [1] ). It is transmitted by the bite of an infected Aedes mosquito. Until now, more than 100 million cases of dengue fever occur worldwide in the Indian subcontinent, Southeast Asia, Southern China, Taiwan area, The Pacific Islands, The Caribbean, Mexico, Africa, Central and South America, Southern United States, and Southern Australia. In Indonesia, dengue cases increase yearly in almost all regions (Rodriguez and Monteiro [2] ). The virus can be spread partly due to an increase in urbanization and also by climate change. Since considerable damage can result from the effects of dengue fever infection, effective control strategies are of vital importance.

A very important aspect of the strategy related to dengue fever spreading is quick and effective action (Rodriguez and Monteiro [2] ). Dengue hemorrhagic fever (DHF) is a more severe form of dengue infection especially if unrecognized and not properly treated in a timely manner. However it has been shown that with good medical management, mortality due to DHF can be less than 1% (Chikaki and Ishikawa [3] ). Patients who are already infected with dengue virus can transmit the infection via the Aedes mosquito just after the first symptoms appear (during 4 - 5 days; maximum 12). Hence, in order to devise effective means of control, it is important to understand the epidemiology of the transmission. Transmission can occur more than once; as a result, if a person has suffered from one virus, there can be a repeat occurrence if a different strain is subsequently involved. Usually people who suffer repeat infections have it worse. They come down with dengue hemorrhagic fever and suffer massive internal bleeding and possible liver damage. The virus causing dengue fever comes in four strains, and immunity to one seems to make infection by a second strain more dangerous (Laurencia et al. [4] ). Experiments for producing and testing those control measures, such as education, antiviral drugs, are costly and time consuming, so any tool, such as a mathematical model(Lasalle [5] ) that will enable us to predict the outcome is highly valuable. Mathematical models avail us with useful predictions about the potential transmission of a disease and the effectiveness of possible control measures. In addition, epidemiology has emerged as an effective tool in disease control. The relationship between mathematics and epidemiology has been increasing. For the mathematician, epidemiology provides new and exciting branches, while for the epidemiologist; mathematical modeling offers an important research tool in the study of the spread of diseases. Bernoulli [6] proposed an epidemiological model which is considered by many authors as the first epidemiological mathematical model. Further work between 1927 and 1933 including those of Kermack and McKendrick [7] largely influenced the development of mathematical epidemiology models (James [8] ). These attempts provided the fundamental framework for the compartmentalization of epidemiological models. Understanding them is vital in gaining important knowledge of the underlying aspects of the dengue fever spread (Thome et al. [9] ). It is also important in assessing the impact of control measures for reducing mortality.

The discovery of antibiotics and vaccines heralded a new hope in disease control. Despite this, new challenges resulting from factors such as drug-resistance have also emerged. Sometimes this led to the emergence of more virulent forms of previously eradicated diseases. For example resistances to such diseases as malaria, tuberculosis, dengue and yellow fever have emerged and, as a result of climate changes, they have been spreading into new regions (Helena &Teresa [10] ). Efforts to cope with this challenge have given rise an increasing trend in the application of mathematical models and interdisciplinary approaches in disease study. Their uses have contributed immensely in decision making and planning in the health sector. In the work reported herein, an SEIRS compartmentalized model is introduced, followed by an optimal control technique in a typical rural environment in Ethiopia. The Pontryagin Minimum Principle is presented as a groundwork for the application of the optimization principle. The resulting system of equations is solved by a fourth-order Runge-Kutta method.

2. Dengue Fever Transmission Model with Education

Quantitative methods are often applied to achieve optimization of investments in the control of a disease. This is necessary in order to obtain maximum benefits from a fixed amount of financial resources. In this case, our efforts will be directed towards the dynamics of the Aedes mosquito vector as well as some management protocols aimed at controlling or alleviating the spread of the disease. Such management principles involving the termination of the reproduction cycle of mosquitoes by avoiding the accumulation of still water in pot-holes and ditches especially after a heavy downpour, are of vital importance as well as educating the local population on issues related to basic hygiene through the television (TV) and radio.

Model Assumptions and Mathematical Formulation

1) The population is uniform and mixes homogeneously. The total population size, N ( t ) = S ( t ) + E ( t ) + I ( t ) + R ( t ) at any time t > 0, where N stands for the total population, E for exposed I for infected, S for susceptible and R for recovered.

2) The natural birth rate b and death rates μn are assumed to be different.

3) Each individual in the population is considered as having an equal probability of contacting the disease with a contact rate β.

4) An infected individual makes contact and is able to transmit the disease with βN per unit time, that is, the contact rate is proportional to the total population size.

5) The fraction of contacts by an infected with a susceptible is S/N. Therefore the number of new infections in unit time per infective becomes (βN)(S/N). This rate is called an infection rate. This gives the rate of new infections or those leaving the susceptible category as (βN(S/N)I = βSI, which is called an incidence of the disease. This type of incidence is called bilinear incidence i.e., proportional to the product of the number of infective individuals and the number of susceptible individuals.

6) The number of infected individuals move from the exposed compartment per unit time is δE at time t.

7) The exposed E move from their compartment to I-compartment at a constant rate δ, so that 1/δ is the mean latent period.

8) The infectious I move from their compartment to R-compartment at a constant rate γ, so that 1/γ is the mean infectious period.

9) The rate of susceptible, exposed, infected and recovered individual removed from each compartments through natural death and disease induced death are μnS, μnE, μnI, μnR and μdI respectively.

10) The recovered individual R move from their compartment to susceptible(S)-compartment at a constant rate α,

11) The differential equations from these assumptions can be represented by a system of ordinary differential equations:

d S d t = b β I S μ n S u S + α R d E d t = β I S μ n E δ E d I d t = δ E ( μ n + μ d + γ ) I d R d t = γ I ( μ n + α ) R + u S (1)

An optimal control strategy aimed at minimizing the objective (cost) functional J of the cost of education for a susceptible population is described by the following differential equations:

d S d t = b β I S μ n S u S + α R d E d t = β I S μ n E δ E d I d t = δ E ( μ n + μ d + γ ) I d R d t = γ I ( μ n + α ) ) R + u S (2)

S ( 0 ) = S 0 0 , E ( 0 ) = E 0 0 , I ( 0 ) = I 0 0 , R ( 0 ) = R 0 0

Subject to: min J ( u ) = 0 T ( A I + 1 2 B u 2 ) d t

where, A is balancing cost factor due to the infective and B is the weight on the cost of education.

Figure 1 is a compartmentalized representation of the mathematical formulation and optimization strategy for education.

Based on the above assumptions, an optimal control problem is formulated by incorporating one of the intervention strategies into our basic mathematical model (see Equations (1) and (2)).

・ u(t) is the control which represents the education ratio of susceptible individuals being educated per unit of time with bounds between 0 and 1.

Figure 1. SEIRS model with education.

・ The inflow of population to the susceptible class is obtained, by combining assumptions 2, 5, 9, 10 and control (education).

・ A number of individuals leaves S and enter E, at the same time, a fraction of exposed E moves to infectious group I with a latent rate δ. δE represents an individual’s move from exposed to infectious. Some of the exposed group die through natural death rate μn, μnE represents movement from exposed to death.

・ Some individuals leave E and enter into the infected individuals I with latent rate δ.

・ A part of the population leaves I and enter the recovered group with recovery rate γ. Combination of assumptions 2, 5, 9, 10 in addition to the control u, gives the rate of recovered.

3. Combination of Education and Treatment by Drug therapy

Antiviral drugs are known to be very helpful in decreasing or preventing disease symptoms at the first sign of a dengue outbreak even when there is no evidence of fever. Before we incorporate drug therapy as part of our treatment protocol and control measures, we will deal with how the application of drug therapy affects some of the model compartments.

・ Consider control variables u1, u1E as representing an individual’s move from exposed to recovered. The exposed populations change per unit of time becomes,

d E d t = β I S μ n E δ E u 1 E (3)

・ In addition, a number of individuals leaves the infected group I and enter the recovered group with recovery rate γ. A number of individuals also leaves the susceptible and exposed groups S and E to enter the recovered group with controls u and u1 respectively. This gives rate of recovered as:

d R d t = γ I ( μ n + α ) R + u S + u 1 E (4)

The differential equation of the diagram for t ≥ 0 is given in a system of ordinary differential equation. Introducing the controls representing the education and drug therapy treatment the model of Equation (1) becomes

d S d t = b β I S μ n S u S + α R d E d t = β I S μ n E δ E u 1 E d I d t = δ E ( μ n + μ d + γ ) I d R d t = γ I ( μ n + α ) R + u S + u 1 E (5)

where, S(0), E(0), I(0), R(0) are the initial conditions. The definitions of above model parameters are listed in Table 1. The control functions, u(t) and u1(t) are bounded, Lebesgue integrable functions(van den Driessche and Watmough [11] ). The control, u1(t), represents the effort on drug therapy treatment of latently infected individuals to reduce the number of individuals that may be infectious. While the control u(t) is the effort on education of susceptible individuals to increase the number of recovered individuals.

A is balancing cost factor due to the infective, B and B1 are the weight on the cost of education and drug respectively. Figure 2 is now the overall representation of the model formulation.

The control problem involves a number of individuals with latent and active dengue fever infections. The cost of applying education and drug therapy treatment controls u(t) and u1(t) are minimized subject to the differential equations (6). The performance specification involves the numbers of individuals with latent and susceptible components respectively, as well as the cost for applying education control (u) and drug therapy treatment control (u1). The objective functional is defined as:

J ( u , u 1 ) = min [ u , u 1 ] 0 T ( A I + 1 2 ( B u 2 + B 1 u 1 2 ) ) d t (6)

where T is the final time and the coefficients, A, B, B1 are balancing cost factors reflecting the importance of the three parts of the objective function. We need to

Table 1. Value of variables and parameters.

Figure 2. SEIRS model with control.

find an optimal control pair, u and u1, such that

J ( u , u 1 ) = min J ( u , u 1 ) | u , u 1 U (7)

where, U = ( u ( t ) , u 1 ( t ) ) | ( u ( t ) , u 1 ( t ) ) measurable, a i ( u ( t ) , u 1 ( t ) ) b i , i = 1 , 2 , t [ 0 , T ] is the control set.

4. Analysis of Optimal Control

The necessary conditions that an optimal pair must satisfy come from the Pontryagins Maximum Principle (Helena [12] ). This principle converts (5) and (6) into a problem of minimizing point-wise Hamiltonian H, with respect to (u,u1). First we formulate the Hamiltonian from the cost functional (6) and the governing dynamics (5) to obtain the optimality conditions. Pontryagin introduced the adjoint function to relate the differential equation to the objective functional. The necessary conditions needed to solve this OC problem, can be followed stepwise:

Step 1: Formulate the Hamiltonian for the problem and by applying Pontryagin’s principle to the Hamiltonian and find optimal controls u*, u 1 with the corresponding solution S*, E*, I* and R* of equation (5).

Step 2: Write the adjoint differential equation, the optimality condition and transversality boundary condition (if necessary). Using the Hamiltonian to find the differential equation of the adjoint λ, and obtain the adjoint variables λ1, λ2, λ3 and λ4 that satisfy adjoint condition.

λ i ( t ) = H x i , where i = 1 , 2 , 3 , 4

Adjoint Functions

λ i = H x λ = ( F x + λ g x ) , adjoint condition (8)

λ 1 = H S λ 1 = λ 1 ( β I + μ n + u ) λ 2 β I λ 4 u

λ 2 = H E λ 2 = λ 2 ( μ n + δ ) λ 3 δ

λ 3 = H I λ 3 = A + λ 1 ( β S ) λ 2 β S + λ 3 ( μ n + m u d + γ ) λ 4 γ

λ 4 = H E λ 4 = λ 1 α + λ 4 ( μ n + α )

with transversality conditions λ i ( T ) = 0 , i = 1 , 2 , 3 , 4 .

The optimality condition is given by,

δ H δ u = 0 at u = u * F u + λ g u = 0 δ H δ u 1 = 0 at u 1 = u 1 F u 1 + λ g u 1 = 0

Step 3: Solve for u* and u 1 in terms of S*, E*, I*, R* and λ

δ H δ u = B u + S ( λ 4 λ 1 ) = 0 (9)

In this way we obtain an expression for the OC:

δ H δ u = 0 at u = u * F u + λ g u = 0

u * = S ( λ 1 λ 4 ) B

δ H δ u 1 = B 1 u 1 + E ( λ 4 λ 2 ) = 0 (10)

u 1 = E ( λ 2 λ 4 ) B 1

Step 4: Solve the four differential equations for S*, E*, I*, R* and λ with boundary conditions, substituting u* and u 1 in the differential equations with the expression for the optimal control from the previous step.

Step 5: After finding the optimal state and adjoint, solve for the optimal control.

We solve that system of differential equations for the optimal state and adjoint. The solution of the optimal control in problem terms of S*, E*, I*, R* and λ, represents the characterization of the optimal control (u*). The state equations and the adjoint equations together with the characterization of the optimal control and the boundary conditions constitute the optimality system.

Remark 1: If the Hamiltonian is linear in the control variable u, it can be difficult to calculate u* from the optimality equation, since δ H δ u would not contain u. Specific ways of solving these kind of problems can be found in Lenhart and John [14] .

Backward-forward Sweep Method

From the model the optimal control problem becomes:

min J ( u ) = 0 T ( A I + 1 2 B u 2 ) d t

Subject to:

d S d t = b β I S μ n S u S + α R d E d t = β I S μ n E δ E d I d t = δ E ( μ n + μ d + γ ) I d R d t = γ I ( μ n + α ) R + u S (11)

With initial value,

S ( 0 ) = S 0 0 , E ( 0 ) = E 0 0 , I ( 0 ) = I 0 0 , R ( 0 ) = R 0 0 and

min J ( u ) = 0 T ( A I + 1 2 ( B u 2 + B 1 u 1 2 ) ) d t

Subject to:

d S d t = b β I S μ n S u S + α R d E d t = β I S μ n E δ E u 1 E d I d t = δ E ( μ n + μ d + γ ) I d R d t = γ I ( μ n + α ) R + u S + u 1 E (12)

As previously indicated, any solution to the above optimal control problem must also satisfy

λ i ( t ) = H x i (13)

where, i = 1, 2, 3, 4, x1 = S, x2 = E, x3 = I, x4 = R

δ H δ u = 0 at u * (14)

δ H δ u 1 = 0 at u 1 (15)

The optimal controls are,

u * = { 0 if δ H δ u < 0 S ( λ 1 λ 4 ) B if δ H δ u = 0 0.9 if δ H δ u > 0 (16)

u 1 = { 0 if δ H δ u 1 < 0 E ( λ 2 λ 4 ) B 1 if δ H δ u 1 = 0 0.9 if δ H δ u 1 > 0 (17)

The optimality condition can usually be manipulated to find a representation of u in terms of t, state variables and λ. If this representation is substituted back into the ODEs for the state variables and λ then the Equations (11) and (12) form a two-point boundary value problem. The Runge-Kutta method is then applied to solve initial value problems, and resolve the optimality system of the optimal control problem. This approach is generally referred to as the Forward-Backward Sweep method. Information about convergence and stability of this method can be found in (Lenhart & John [14] ). The process begins with an initial guess on the control variable. Then, the state equations are simultaneously solved forward in time and adjoint equations are solved backward in time. The control is updated by inserting the new values of states and adjoint into its characterization, and the process is repeated until convergence occurs.

5. Numerical Illustrations and Conclusions

Numerical solutions to the optimality system comprising the state Equation (5)and adjoint equations are carried out using MATLAB and using parameters in Table 1 and the following weight factors and initial conditions: A = 100, B = 0.04, B1 = 0.06, S(0) = 86.46%, E(0) =4.5%, I(0) = 9.042%, R(0) = 0%. The algorithm is the forward-backward scheme; starting with an initial guess for the optimal controls u and u1, the state variables are then solved forward in time from the dynamics (5) using a Runge-Kutta method of the fourth order. Then those state variables and initial guess u and u1 are used to solve the adjoint equations backward in time with given final conditions (16) and (17) by employing a fourth order Runge-Kutta method. The controls u and u1 are updated and used to solve the state and then the adjoint system. This iterative process terminates when current state, adjoint, and control values converge sufficiently (Helena &Teresa [10] ).

5.1. Results for Optimal Education Only

With this strategy, education (u) is utilized in the disease control while the control on drug therapy treatment (u1) is set to zero, with weight factors B1 = 0, A = 100, B = 0.04. For this strategy, we observed that the number of susceptible individuals is higher when education and drug therapy treatment are absent (Figure 3). For the latently exposed (E) individuals in Figure 4, it can be seen that with the presence of education the percentage rate of the exposed is lower than when there is no education. The same trend is followed in Figure 5, where the percentage of the infected group (I) is lower when exposed to education. However the percentage of the recovered individuals (R) with education is higher than when there is no exposure to education. Figure 5 and Figure 6 are respectively less and greater respectively than the percentage of infected individuals and recovered individuals in the absence of education and drug therapy treatment.

5.2. Optimal Drug Therapy Treatment Only

The control (u1) on drug therapy treatment is utilized while the control on education(u) is set to zero, with weight factors A = 100, B = 0.04, B1 = 0.06. For this strategy, it can be observed in Figure 7, that controls with education and drug therapy treatment lowers the percentage of susceptible individuals (hardly perceptible in the diagram) than with education alone. This is because the recovered individuals go back to susceptible group and increase the susceptible group at higher rate. For the latently infected individuals in Figure 8, it can be seen that in the absence of education, and with an initially exposed population of 4.5%; there is hardly any change in the percentage of the individuals exposed both with

Figure 3. Susceptible with education Vs without education.

Figure 4. Exposed with education Vs without education.

Figure 5. Infectives with education Vs without education.

Figure 6. Recovered with education Vs without education.

Figure 7. Susceptibles with education Vs with education anddrug therapy treatment.

Figure 8. Exposed with education Vs with education and drug therapy treatment.

education and with education and treatment during the first 5 weeks. It is obvious that the impact of education takes time to be felt or manifested in the dynamics. However there is a dramatic change in the dynamics after this period as the percentage of the exposed with education and treatment becomes significantly lower than for those with education alone. For the infected individuals in Figure 9, with an initially infected population of 9.04%, it can be seen that using both intervention mechanisms is better than using education as only control mechanism. As earlier observed, there is a time lag of about ten weeks for the impact of education to be reflected in the dynamics. The same trend is observed in Figure 10 for the percentage of the recovered where the time lag for education is about five weeks before the influence of education with treatment shows a higher percentage than with education alone.

5.3. Optimal Education and Drug Therapy Treatment

With this strategy, the controls on education (u) and drug therapy treatment (u1) are utilized, with weight factors A = 100, B = 0.04, B1 = 0.06. Figure 11 shows that the percentage of susceptible individuals with education and treatment is lower than the susceptible population in the absence of education and drug therapy treatment. Figure 12, shows that without control the percentage of exposed individuals is higher than would be the case with education and treatment options. The positive effect of treatment and education is further confirmed in Figure 13 where there is a higher percentage of individuals recorded without any control measures. Figure 14 shows that as more people get

Figure 9. Infectives with education Vs with education and drug therapy treatment.

Figure 10. Recovered with education Vs with education and drug therapy treatment.

Figure 11. Susceptibles without education and drug therapy treatment Vs with education.

exposed to treatment and education the more they are likely not to get infected.

Figure 12. Exposed without education and drug therapy treatment Vs with education.

Figure 13. Infectives without education and drug therapy treatment Vs with education and drug therapy treatment.

Figure 14. Recovered without education and drug therapy treatment Vs with education and drug therapy treatment.

6. Concluding Remarks

The results displayed herein not only confirm the validity of the mathematical formulation derived but also illustrate how to optimally apply control measures involving treatment and education for the control of dengue fever. Utilizing education and drug therapy treatment lead to better disease control in the population than utilizing drug therapy treatment only. In addition, the application of only one form of control measure though it results in a delayed peak in the percentage of exposed and infected, is not as effective as using both controls. Thus control programs that specialize in an optimal application of multi-control measures can effectively reduce or alleviate the effects of dengue fever spread.

Further work should include other control variables like the effect of bio-immunology on the spread of dengue fever, the use of medicated mosquito nets, development and application of vaccines, creation of sterile mosquito males for the control of mosquito population etc.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

Cite this paper

Onyejekwe, O.O., Tigabie, A., Ambachew, B. and Alemu, A. (2019) Application of Optimal Control to the Epidemiology of Dengue Fever Transmission. Journal of Applied Mathematics and Physics, 7, 148-165. https://doi.org/10.4236/jamp.2019.71013

References

  1. 1. Noorani, M.S.M. (2102) SEIR Model for Transmission of Dengue Fever in Selangor Malaysia. Department of Mathematics, University of Kebangsaan Malaysia, Bangi.

  2. 2. Rodrigues, H.S., Monteiro, M.T.T., and Torres, D.F.M. (2013) Sensitivity Analysis in a Dengue Epidemiological Model. Conference Papers in Mathematics.

  3. 3. Chikaki, E. and Ishikawa, H. (2009) A Dengue Transmission Model in Thailand Considering Sequential Infections with All Four Serotypes. Journal Inf. Dev. Ctries, 3, 711-722.

  4. 4. Laurencia, N.M., Estomih, S.M. and Makinde, O.D. (2015) Temporal Model for Dengue Fever with Treatment. Advances in Infectious Diseases, 5, 21-36.

  5. 5. LaSalle, J.P. (1976) The Stability of Dynamical Systems. Regional Conference Series in Applied Mathematics, SIAM Philadelphia.

  6. 6. Bernoulli, D. (1776) Essai d’une Nouvelle Analyse de la Moralite’ Causee Par La Petite Verole & des Avantage de L’ inoculation Pour La Prevenir. Mem. Math. Phys. Acad. Roy. Sci. Paris, 145.

  7. 7. Kermack, W. and McEndrick, A. (1991) Contributions to the Mathematical Theory of Epidemics, II the Problem of Endemicity. Bulletin of Mathematical Biology, 53, 57-87.

  8. 8. James, D.M. (2002) Mathematical Biology. Springer-Verlag.

  9. 9. Thome, R.C., Yang, H.M. and Esteva, L. (2010) Optimal Control of Aedes, Aegypti Mosquitoes by the Sterile Insect Technique and Insecticide. Mathematical Biosciences, 223, 12-23.

  10. 10. Helena, S.R. and Teresa, M. (2012) Sensivity Analysis in a Dengue Epidemiological Model, School of Business, Polytechnic Institute of Vianna do Castello, Valencia, Portugal, R and D Center Algoritmi, Department of Production and Systems, University of Minho, Braga, Portugal, Rand D Center CIDMA, Department of Mathematics, University of Aveiro, Aveiro.

  11. 11. van den Driessche, P. and Watmough, J. (2002) Reproduction Numbers and Sub-Threshold Endemic Equilibria for Comparmental Models of Disease Transmission. Mathematical BioSciences, 180, 29-48. https://doi.org/10.1016/S0025-5564(02)00108-6

  12. 12. Helena, S.F.R. (2012) Optimal Control and Numerical Optimization Applied to Epidemiological Models. PhD Thesis University of Aveiro.

  13. 13. Esayas, Z.K. (2015) Epidemiological Modelling of Measles Disease with Optimal Control of Vaccination Strategy. MSc Thesis Computational Science Program, Addis Ababa University, Addis Ababa.

  14. 14. Lenhart, S. and John, T.W. (2007) Optimal Control Applied to Biological Models, Mathematical and Computational Biology Series. Chapman and Hall/CRC, Boca Raton.