www.biomathforum.org/biomath/index.php/biomath ORIGINAL ARTICLE Determination of parameters for Cauchy’s problem for systems of ODEs with application to biological modelling A.N. Pete∗, P. Mathye∗, I. Fedotov∗ and M. Shatalov∗† ∗Department of Mathematics and Statistics Tshwane University of Technology, Pretoria, South Africa mathyeph@gmail.com; FedotovI@tut.ac.za; ShatalovM@tut.ac.za †Manufacturing and Materials Council for Scientific and Industrial Research (CSIR), Pretoria, South Africa Received: 13 March 2015, accepted: 23 April 2016, published: 9 May 2016 Abstract—An inverse numerical method that esti- mates parameters of dynamic mathematical models given some information about unknown trajectories at some time is applied to examples taken from Biology and Ecology. The method consists of de- termining an overdetermined system of algebraic equations using the experimental data. The solution of the overdetermined system is then obtained using, for example the least-square method. To illustrate the effectiveness of the method an analysis of ex- amples and a numerical example for the model that monitors the dynamics of HIV is presented. Keywords-Inverse problem; least squares meth- ods; parameter estimation; dymamic systems; predator-prey system; I. INTRODUCTION Function approximation on a fixed interval by means of an initial value problem of an ordinary differential equation with unknown coefficients and unknown initial values as presented by M Shatalov, I. Fedotov and S.V. Joubert in [7], is central to this study. In their paper al method to determine both the unknown coefficients and initial values of a dynamic system by minimizing a certain goal function is presented. In earlier collaborations Shatalov and Fedotov suggested the use of such an approach in identifying dynamic systems’ parameters from experimental data, see [6]. Several other parameter estimation methods are presented in literature, for instance the stochastic models; the Bayesion approach, the Monte Carlo technique, the numerical method with combined Adomain/Alienor approach, the differential evo- lution (DE) and the hybrid Taguchi-differential evolution algorithm. The method used is based on integrating both sides of equations of a dynamic system, and ap- plying regression methods to the overdetermined system of linear algebraic equations with possible constraints. The unknown parameters and initial values can then be obtained using the method of least squares. In this paper, the proposed method gives parameter estimates that have a percentage relative error that is mostly less than 0.4% for artificially generated data and parameters that are in the expected range for real data. As an illustration of the method of identifica- Citation: A.N. Pete, P. Mathye, I. Fedotov, M. Shatalov, Determination of parameters for Cauchy’s problem for systems of ODEs with application to biological modelling, Biomath 5 (2016), 1604231, http://dx.doi.org/10.11145/j.biomath.2016.04.231 Page 1 of 7 http://www.biomathforum.org/biomath/index.php/biomath http://dx.doi.org/10.11145/j.biomath.2016.04.231 A.N. Pete et al., Determination of parameters for Cauchy’s problem ... tion, we present the analysis for problems from Biology and Ecology. The following mathematical models are taken as examples: model of a free pop- ulation (with negligible mortality and with mortal- ity different from zero), population with negligible mortality and unknown initial conditions, inter- species modifications of the Lotka-Voltera model and a model that monitors the dynamics of HIV. Numerical results with artificially generated data and real data is then given for the model that monitors the dynamics of HIV. A. Inverse methods The general approach to the inverse system identification of a n –dimensional parameter a ∈ A ⊂ Rn, where A can coincide with Rn (no constrains between entries of a) and A can be subset of Rn (there are constrains) in a system of ordinary differential equations of the following form ẋ = f(t,x,a), (1) where x is vector-function [0,T] 3 t −→ x (t) ∈ Rm subject to experimental information concern- ing the values x (tj) at the point tj ∈ [0,T] is known: t0 · · · tj · · · tN x0 = x(t0) · · · xj = x(tj) · · · xN = x(tN ) (2) The general approach to find a solution of the formulated problem (1) consists of determining an overdetermined system of algebraic equations using the experimental data (2) Aa = h, (3) with respect to unknown vector a. The solution of the system (3) is then obtained using any method of solution of the overdetermined system, for example the least-square method which minimize the difference Aa − h using Euclidean metric. It is known that in this case (see, for example [5]), the solution of (3) can be obtained by solving the following system A>Aa = A>h, (4) to obtain a. II. MATHEMATICAL MODELS 1) Problem 1: Free population: Consider a single species that grows by sexual reproduction. Assume that the individuals move in the pop- ulation like Brownian motion particles (or that the population is colonial), then the frequency of contact between the individuals is proportional to the squared population density [1]. Further assuming that mortality is different from zero and is independent of the population size, the scalar function f(t,x,a) is given by f(t,x,a) = a1x 2 x + a2 −a3x, (5) with the unknown vector a = (a1,a2,a3)> ∈ R3, n = 3, m = 1 and a1 > a3. The parameters a1 and a2 represent per capita birth rate (fecundity) and the population at which half of the females are able to reproduce, respectively. The mortality rate of the population is represented by a3. The parameters a1, a2 and a3 are positive and x(t) = x(t) ∈ R is a scalar function. We also assume that x0 is specified. A special case arises if we assume that the population is of negligible mortality. That is, if a3 = 0 the scalar function f(t,x,a) is then given as f(t,x,a) = a1x 2 x + a2 . (6) 2) Problem 2: Population with negligible mor- tality and x0 unknown: Here we consider the special case of Problem 1 (that is, if a3 = 0) but with the value x0 considered as an incorrect value which must be corrected. The statement of such a problem naturally appears since the initial value x0 = x(t0) plays an important role, namely; it defines the Cauchy’s problem for equation (1) and, secondly x0 is included widely in computations below. Therefore if x0 given from an experiment has low accuracy, it would be desirable to define this value with more accuracy. 3) Problem 3:Nonlinear predation at small prey population: In this problem we consider an inter- species modification of the Lotka-Voltera model Biomath 5 (2016), 1604231, http://dx.doi.org/10.11145/j.biomath.2016.04.231 Page 2 of 7 http://dx.doi.org/10.11145/j.biomath.2016.04.231 A.N. Pete et al., Determination of parameters for Cauchy’s problem ... where there is a nonlinear predation at small prey population. Consider the system ẋ = −a1 x2y xp + a2 + a3x, ẏ = a4 x2y xp + a2 −a5y,   (7) where p = 1, 2, a1,a2,a4 > 0, a3,a5 ≥ 0 and n = 5, m = 2. The prey and predator densities are represented by x(t) and y(t), respectively. In this model, the parameters a1 represent the rate of the consumption of prey by the predator popula- tion, a2 the prey population density at which the predator’s consumption is half the maximum value [1], a3 the prey’s growth rate, a4 is rate at which the prey contributes to the predator’s growth rate and a5 is the predator’s death rate. 4) Problem 4: A model that monitors the dy- namics of HIV: Consider the following two- dimensional model that monitors the dynamics of HIV. The model considers two sub-populations: HIV susceptible (x), the HIV infected population (y). The total population size is given by N = x + y. The model is described by ẋ = −a1 xy N −a2x + a3, ẏ = a1 xy N −a4y,   (8) where a1, a2, a3 and a4 are all positive con- stants model parameters. The parameter a1 and a2 denotes the average rate of infection by HIV, a2 the natural cessation of sexual activity, a3 the recruitment rate of susceptible and a4 denotes the death rates of the infected population due to HIV. This model is a modified version of the three dimensional one by Gumel (see, [3]). The vector a = (a1,a2,a3,a4) T ∈ R4 is unknown. Note that in this case n = 4 and m = 2. III. CONSTRUCTION OF OVERDETERMINED SYSTEMS Consider equation (1) where f is defined by f(t,x,a) = a1x 2 x + a2 , (9) The resultant equation can be rewritten as ẋx + a2ẋ = a1x 2, (10) Integration of (10) with respect to t from t0 to tj (j = 1, 2, . . . ,N) gives −a1Pj + a2∆xj = −∆hj, (11) where ∆xj = x (tj) −x (t0) , Pj = ∫ tj t0 x2dt, ∆hj = [ 1 2 x2(tj) − 1 2 x2(t0) ] .   (12) The integral Pj can be calculated by using a quadrature rule, for example trapezoidal rule. Thus, System (3) is solved with A =   −P1 ∆x1... ... −PN ∆xN   , a = ( a1 a2 ) and h =   ∆h1... ∆hN   . Now, suppose that (1) is defined as in (10) and the initial condition is unknown. Let us replace in (11) ∆xj and ∆hj by those given in (12): −a1Pj + a2xj − ( a2x0 + 1 2 x20 ) = −hj, (13) where hj = 12x 2 j under the assumption that for Pj we use the old values of x0 defined by (2) for j = 1 and for j ≥ 2 the values of Pj are evaluated by an open quadrature rule. Formally speaking system (13) is nonlinear but setting − ( a2x0 + 1 2 x20 ) = a3, (14) we obtain the linear system of the form (3) with A =   −P1 x1 1... ... ... −PN xN 1   , a =   a1a2 a3   ,h =   −h1... −hN   . Biomath 5 (2016), 1604231, http://dx.doi.org/10.11145/j.biomath.2016.04.231 Page 3 of 7 http://dx.doi.org/10.11145/j.biomath.2016.04.231 A.N. Pete et al., Determination of parameters for Cauchy’s problem ... The value of x0 can then be obtained from Equa- tion (14). Equation (1) with f defined by (5) can be written in the form xẋ + a2ẋ = a1x 2 −a3x2 −a2a3x, (15) In contrast to (10) two supplementary terms −a3x2 and −a2a3x arise. Integrating (15) we obtain (a1 −a3)Pj −a2xj + ( 1 2 x20 + a2x0 ) −a2a3Xj = hj, (16) where Xj = ∫ tj t0 xdt. Introducing the new vector b = (b1,b2,b3,b4) where b1 = a1 −a2, b2 = a2, b3 = 1 2 x20 + a2x0 and b4 = a2a3 the system can be written in the form (3) where A =   P1 −x1 1 −X1... ... ... ... PN −xN 1 −XN   , b =   b1 b2 b3 b4   and h =   h1... hN   In other words, we obtain the following system b1Pj − b2xj + b3 − b4Xj = hj, (j = 1, . . .N) (17) Since the Jacobian is ∂b ∂a = a2 + x0 6= 0, we can find the vector b and hence the unknown vector a. The system (7) of an interspecies modification of the Lotka–Voltera model for nonlinear predation at small prey population can be written as ẋxp + a2ẋ = −a1x2y + a3xp+1 + a2a3x ẏxp + a2ẏ = a4x 2y −a5yxp −a2a5y } (18) where p = 1, 2, a1,a2,a4 > 0 and a3,a5 ≥ 0. Integrating (18) we obtain ∆hj = −a1Zj −a2∆xj + a3Pj + a2a3Xj Sj = −a2∆yj + a4Zj −a5Qj −a2a5Yj } (19) where,using the same notation as in (12): ∆xj = x (tj) −x (t0) , Pj = ∫ tj t0 x2dt, ∆hj = [ 1 2 x2(tj) − 1 2 x2(t0) ] .   (20) and Xj = ∫ tj t0 x(t)dt,Yj = ∫ tj t0 y(t)dt, (21) Zj = ∫ tj t0 y(t)x2(t)dt,Qj = ∫ tj t0 xp(t)y(t)dt, Sj = ∫ tj t0 xp(t)dt, (22) with j = 1 . . .N. Introducing the new vector b = (b1,b2,b3,b4,b5,b6,b7) > where b1 = a1, b2 = a2, b3 = a3, b4 = a4, b5 = a5, b6 = a2a3 and b7 = a2a5, we write the System (19) in the form (3) where A=   −Z1 −∆x1 P1 0 0 X1 0 −Z2 −∆x2 P2 0 0 X2 0 ... ... ... ... ... ... ... −ZN −∆xN PN 0 0 XN 0 0 −∆y1 0 Z1 −Q1 0 −Y1 0 −∆y2 0 Z2 −Q2 0 −Y2 ... ... ... ... ... ... ... 0 −∆yN 0 ZN −QN 0 −YN   and h = (∆h1, . . . , ∆hN,S1, . . . ,SN ). In other words, we obtain the following system of 2N equations with seven unknowns: ∆hj = −b1Zj − b2∆xj + b3Pj + b6Xj, Sj = −b2∆yj + b4Zj − b5Qj − b7Yj, } (23) subject to −b6 + b2b3 = 0, −b7 + b2b5 = 0. } (24) To evaluate numerically the right hand side of Sj can be considered as the Riemann–Stietjes integral Sj = ∫ tj 0 xp(t)dy(t), (see, for example Dragomir and Fedotov [3]). Otherwise this integral can be computed using numerical differentiation. Biomath 5 (2016), 1604231, http://dx.doi.org/10.11145/j.biomath.2016.04.231 Page 4 of 7 http://dx.doi.org/10.11145/j.biomath.2016.04.231 A.N. Pete et al., Determination of parameters for Cauchy’s problem ... Another approach to construct the over deter- mined system is generally appealing since the system obtained is simpler than (19). This method is possible due to the special form of the system (7); namely the system can be written as ẋ = −a1g(x,y,a2) + a3x, ẏ = a4g(x,y,a2) −a5y, } (25) where g(x,y,a2) = x2y xp+a2 . Eliminating g(x,y,a2) we get ẏ + c1y + c2ẋ− c3x = 0, (26) where a5 = c1, a4 a1 = c2 and a4a3 a1 = c3. (27) Using the integration techniques that we discussed earlier, we determine the vector c = (c1,c2,c3) > and hence c3 c2 = a3, which will be considered as unknown in the next step. Multiplying the first equation of (7) by (xp + a2), we obtain xpẋ + a2ẋ = −a1x2y + a3xp+1 + a2a3x. (28) After integration of (28) the overdetermined sys- tem can be written as a1Zj + a2 (∆xj −a3Xj) = −∆hj + a3Pj, (29) where j = 1, . . . ,N, and Zj, Pj, Xj, ∆hj, ∆xj were defined in (20) and (22). The unknown parameters a1 and a2 may then be determined from equation (29) and a4 and a5 from (27). A. A model that monitors the dynamics of HIV The system (8) can be written as ẋ = −a1Pj −a2Xj + a3, ẏ = a1Pj −a4Yj. } (30) where ∆xj = x(tj) −x(t0), ∆yj = y(tj) −y(t0), Pj = tj∫ 0 x(t)y(t) x(t) + y(t) dt, Xj = tj∫ 0 x(t)dt and Yj = tj∫ 0 y(t)dt. The system (30) can be written in the form (3) with A =   −P1 −X1 1 0 −P1 −X2 1 0 ... ... ... ... P1 0 0 −Y1 P1 0 0 −Y2 ... ... ... ...   and h = (∆x1, ∆x2, . . . , ∆y1, ∆y2, . . .) > IV. NUMERICAL EXTRACTION OF MODEL PARAMETERS In order to illustrate the effectiveness of the method the parameter identification of the system (8) is presented. Using parameter values from Gumel [4], we generate points of solutions of the system (8) by the adaptive Runge-Kutta method. A mathematical software Mathcad was used for the adaptive Runge-Kutta method and for the in- tegration by quadratue rules. These solutions are then perturbed by a normal distribution with mean x̄ and standard deviation δ = 0.5 and further taken as “experimental data”. Furthermore, system (8) is studied in the context of the Gauteng Province, South Africa. Data used in the numerical simulation is obtained from the Acturial Society of South Africa (ASSA) (see [2]) HIV prevalence estimates. The data is compiled from the population sensus, antenatal survey and registered deaths [2]. The ASSA 2003 tables gives the population N and the number infected with HIV y. From the relationship N = x + y, (31) we obtain the susceptible population x as x = N −y. (32) The results given in Table 2 below, show a comparison of the parameters used (from Gumel [4]) , the estimated parameters and the percentage error given by ||α− α̃||/||α||× 100. Biomath 5 (2016), 1604231, http://dx.doi.org/10.11145/j.biomath.2016.04.231 Page 5 of 7 http://dx.doi.org/10.11145/j.biomath.2016.04.231 A.N. Pete et al., Determination of parameters for Cauchy’s problem ... Table 2: Problem 4: The parameter estimates and errors. Parameter Actual value Estimated value % Error α α̃ a1 0.24 0.2404 0.375 a2 0.03125 0.0313 0.16 a3 2000 2000 0.000 a4 0.531 0.532 0.16 Applying the parameter estimation method in the case of the transmission dynamics of HIV in the Gauteng Province, South Africa the following estimates are obtained:  a1 a2 a3 a4   =   0.393 0.049 565143 0.149   . Figures 1 and 2 below show the comparison of the estimated solutions with the “experimental data” and the absolute error in the decimal loga- rithmic scale. The initial value problem was solved with the new coefficients and old initial conditions. 0 10 20 2 .10 4 4 .10 4 X i X1i ti 0 10 20 2 1 0 log Xi X1i−( ) ti Fig. 1. Comparison of the estimated solutions, Xi, with the “experimental data” and the absolute error in the decimal logarithmic scale. 0 10 20 100 200 Yi Y1i ti 0 10 20 4 2 log Yi Y1i−( ) Fig. 2. Comparison of the estimated solutions, Yi, with the “experimental data” and the absolute error in the decimal logarithmic scale. From the comparison of the parameters in Table 2, Figures 1 and 2, it can be seen that the estimated parameters are close enough to the actual values. The estimated parameters of dynamic HIV math- ematical models for the Gauteng Province are all nonnegative and also in the expected ranges. V. CONCLUSION AND SUGGESTIONS FOR FURTHER RESEARCH A method for estimating parameters of dy- namic mathematical models given some informa- tion about unknown trajectories at some time, by Shatalov, Fedotov and Joubert [7], was applied to problems from Biology and Ecology. In Problem 2 where the initial value was assumed to be unknown with accuracy, the method was applied to find this value. To illustrate the efficiency of the method, a numerical example for the model that monitors the dynamics of HIV was presented. The estimated parameters for the artificially generated data are close enough to the actual parameter values and for the Gauteng Province the estimated parameters Biomath 5 (2016), 1604231, http://dx.doi.org/10.11145/j.biomath.2016.04.231 Page 6 of 7 http://dx.doi.org/10.11145/j.biomath.2016.04.231 A.N. Pete et al., Determination of parameters for Cauchy’s problem ... are in the expected ranges. The method for estimating parameters values could be improved by incorporating a suitable penalty term that minimizes the error caused by numerical quadrature and high observational noise levels in the real data. An improvement of the method will be investigated in future studies. REFERENCES [1] A. D. Bazykin, Nonlinear Dynamics of Interacting Pop- ulations, Singapore, World Scientific, 1998. [2] R. E. Dorrington, L. F. Johnson, D. Brandshaw & T. Daniel, The demographic impact of HIV/AIDS in South Africa: National and provincial indicators, Cape Town: Centre for Acturial Research, South African Medical Research Council and Actuarial Society of South Africa, 2006. [3] S. S. Dragomir & I. Fedotov, A Gruss type inequality for mapping abounded variation and applications for numerical analysis, Nonlinear Functions. Anal., Appl., 6(3) : 425−433, 2001. [4] A. B. Gumel, A competitive numerical method for a chemotherapy model of two HIV subtypes, Applied Math- ematics and Computations, 131 : 329−337, 2002. [5] C. L. Lawson & R. J. Hanson, Solving Least Square Problems, New Jersey, America: Prentice–Hall, 1974. [6] M. Shatalov & I. Fedotov, On identification of dynamic systems parameters from experimental data, RGMIA, Victoria University, 10(1, 2) : 106−116, 2007. [7] M. Shatalov,I. Fedotov & S. V. Joubert, A novel method of interpolation and extrapolation of functions by a linear initial value problem, Buffelsfontein TIME 2008 Peer- review Conference Procedings, 22−26 September, 2008. Biomath 5 (2016), 1604231, http://dx.doi.org/10.11145/j.biomath.2016.04.231 Page 7 of 7 http://dx.doi.org/10.11145/j.biomath.2016.04.231 Introduction Inverse methods Mathematical models Problem 1: Free population Problem 2: Population with negligible mortality and x0 unknown Problem 3:Nonlinear predation at small prey population Problem 4: A model that monitors the dynamics of HIV Construction of overdetermined systems A model that monitors the dynamics of HIV numerical extraction of model parameters Conclusion and suggestions for further research References