Original article Biomath 1 (2012), 1209032, 1–6 B f Volume ░, Number ░, 20░░ BIOMATH ISSN 1314-684X Editor–in–Chief: Roumen Anguelov B f BIOMATH h t t p : / / w w w . b i o m a t h f o r u m . o r g / b i o m a t h / i n d e x . p h p / b i o m a t h / Biomath Forum Model and Simulations of a Wood Frog Population Nofe Al-Asuoad∗, Roumen Anguelov†, Keith Berven ‡ and Meir Shillor∗ ∗ Department of Math. and Statistics, Oakland University, Rochester, MI, USA Emails: nalasuoa@oakland.edu, shillor@oakland.edu † Department of Math. and Applied Mathematics, University of Pretoria, Pretoria, South Africa Email: roumen.anguelov@up.ac.za ‡ Department of Biological Sciences, Oakland University, Rochester, MI, USA Email: berven@oakland.edu Received: 15 July 2012, accepted: 3 September 2012, published: 11 October 2012 Abstract—This work presents and simulates a mathe- matical model for the dynamics of a population of Wood Frogs. The model consists of a system of five coupled impulsive differential equations for the larvae, juveniles (early, middle, and late) and the mature adult populations. A simulation result depicts possible dynamics of the frogs’ population when during one year the larvae population dies out. This provides a tool for the study of the resilience of the population and the conditions that may lead to its survival and flourishing or extinction. Keywords-population dynamics; compartmental model; Wood Frog, impulsive ODEs; simulations I. INTRODUCTION We present a model for a Wood Frog population and a preliminary simulation of its solutions. This research is motivated by more than two decades of field observations of a population of Wood Frogs, which has been recently reported in Berven [3]. The aim is, once the model is val- idated by comparison with experimental data collected in [3], to study the conditions that allow for the survival, and possible flourishing of the population. The model is of the compartmental type (see, e.g., [1], [5], [6] and the references therein) and consists of a system of five nonlinear ordinary differential equations (ODEs), and includes impulses that describe the tran- sitions from one population to the next. The equations describe the dynamics of the larval aquatic stage, and juvenile and adult stages, which are terrestrial, in the development of the frogs’ population. When the larvae metamorphose and become juveniles, they leave the pond over a period of two weeks, and it is found in [3] that there is considerable merit in dividing the juvenile population into three groups, those who leave the pond early, late, and in the middle. In this manner we obtain the five compartments with the associated impulsive ODEs. Whereas the applied interest in the model, once it has been validated by comparison with the data from the field, is to study the conditions for the survival of the population, the mathematical interest lies in the facts that the aquatic larval stage is separate from the other stages, and the interactions are via transfer conditions at prescribed times, and the resulting impulses. The long-time interest in the model lies in its ability to provide for qualitative and quantitative predictions on the overall populations growth that will allow to better understanding and management of the populations. The model is described in the following section, then, we present the results of a typical computer simulation of the model which shows the dynamics and the recovery from a year without any larvae. In the conclusions section we also mention some unresolved questions that we plan to address in the future. II. THE MODEL We construct a model for the dynamics of a Wood Frog population. The stages of the life cycle of the Citation: N. Al-Asuoad , R. Anguelov, K. Berven, M. Shillor, Model and Simulations of a Wood Frog Population, Biomath 1 (2012), 1209032, http://dx.doi.org/10.11145/j.biomath.2012.09.032 Page 1 of 6 http://www.biomathforum.org/biomath/index.php/biomath http://dx.doi.org/10.11145/j.biomath.2012.09.032 N. Al-Asuoad et al., Model and Simulations of a Wood Frog Population frogs that we concentrate on are the larvae, juveniles, and mature adults. The first one is aquatic, and the other two are terrestrial, a fact that leads to a nonstan- dard model with impulses and possibly time dependent periodic coefficients. The field data obtained from the population study in Berven [3] allows us to use ODEs. The model deals with the total populations because the spatial distribution of the populations is not taken into account. The frogs’ life cycle we model is as follows. The eggs are produced in very large numbers in the spring, over a period of two weeks, and those that survive become larvae. The larvae that survive undergo metamorphosis in the summer and become juveniles. These become mature adults over the next 1-4 years. We note that since the hatching rate of the eggs is constant, at about 90%, for the sake of simplicity we omit the eggs compartment. An interesting observation made in [3] leads us to split the juvenile population into three groups, those who leave the pond at the beginning of the first week, and the beginning and end of the second week, since these juveniles’ rates of growth and maturation are different. We denote by L and M the total numbers of larvae, and mature frogs, respectively, and by Je, Jm and Jl the three subpopulations, early, middle, and late, of the juveniles, all functions of time t (measured in days). We let [0, T ] denote the time interval over which the populations grow, or have been under observation. We describe the rate of change of each population per day. The model is of the compartmental type, and is depicted schematically in Fig. 1. - L R µL Q Q QQs -� � ��3 Je Jl Jm �µe R µm Rµl H H H Hj � � � �* - M R µM Fig. 1. Compartmental structure and flow chart The eggs hatch within two weeks after fertilization. The larval period lasts about eight to ten weeks. So counting from the laying and fertilization of the eggs, the early larvae metamorphose into juveniles within 10 weeks, and the later ones one and two weeks later. Since these periods are relatively short, compared to the rest of the dynamics of the population, e.g., the life span or growing to maturity, we model these changes using impulsive differential equations (see, e.g., [1] and also [4], [2] and the references therein). We start the time count t = 0 on the day of the eggs’ fertilization and assume that initially the number of larvae that successfully hatch is L0, and so the larvae population undergoes a discontinuous change on that day by jumping from no larvae to L0. Then, at the same day at the year k (for k = 1, 2, ..., T ), that is at the times tk = 365k, all the eggs hatch and L jumps from zero larvae (before hatching) to the number that hatched, σM (tk), which is proportional to the mature population M at that time. The proportionality rate constant σ is the fertility rate of the mature female frogs (which are a third of the matures). The splitting of the juveniles into the three groups of early, middle, and late ones is based on the observation that when the larvae population is large, those who metamorphose and leave earlier develop and mature faster. We denote by τe, τm, and τl the respective days in the year on which the first second, and third group of larvae become early Je, middle Jm, and late Jl juveniles, respectively. The data in [3] indicates that we may set approximately τe = 75, τm = 80, and τl = 85 days, but these choices are somewhat arbitrary, however, in the model we keep the general notation. We assume that at the times tk + τe the fraction δeL becomes ‘early’ juveniles, at the times tk + τm the fraction δmL becomes ‘middle’ juveniles, and at tk + τl the rest of the larvae become ‘late’ juveniles. At the exceptional times the larvae population jumps discontinuously, that is impulses take place. At the beginning of each year, at time tk, for k = 1, 2, 3, . . . , the larvae population is L(tk + 0) = σM (tk). At times tk + τe, tk + τm, and tk + τl we have, L(tk + τe + 0) = (1 − δe)L(tk + τe − 0), L(tk + τm + 0) = (1 − δm)L(tk + τm − 0), L(tk + τl + 0) = 0. The equation for the larvae population is dL dt = −µLL, t 6= tk, tk + τe, tk + τm, tk + τl, for k = 1, 2, 3, . . . . Here, the mortality rate is µL = µL(t, L) = µ1L + µ2L(t)L(t), where, following [7], we let µ1L represent the density independent part and µ2L(t)L(t) is the density dependent Biomath 1 (2012), 1209032, http://dx.doi.org/10.11145/j.biomath.2012.09.032 Page 2 of 6 http://dx.doi.org/10.11145/j.biomath.2012.09.032 N. Al-Asuoad et al., Model and Simulations of a Wood Frog Population part that depends on the available resources and, hence, may be time dependent. We turn to the describe the rates of growth of the ju- venile populations. In normal circumstances, e.g., when food is sufficient or the weather is mild, a juvenile is ready for reproduction in the next mating and egg laying season, that is the next spring. However, some, especially the late ones, may become fertile only in the second year. We assume that the different juvenile populations have different mortality rates, set as µr(t) = µ1r + µ2r(t)Jr(t), for r = e, m, l, where µ1r represent the density indepen- dent rates, and µ2rJr are the density dependent parts. We denote by αe, αm, and αl the rates at which the juveniles mature and move to the M population, which may be time dependent, to take into account possible changes in the environmental conditions and, also, are likely to be periodic functions reflecting the availability of food and growth rates of the juveniles. We denote by βrs (r, s = e, m, l, r 6= s) the influence of population Ls on Lr, which may describe the competition for food. However, the experimental data is not clear about it, so we assume that these rate coefficients are small. Thus, the population growth equations for Je are dJe dt = −µe(t)Je(t) − αe(t)Je(t) − βemJm(t)Je(t) −βel Jl(t)Je(t), t 6= tk + τe, and similar equations hold for Jm and Jl. Here, k = 0, 1, 2, . . . is the kth year. Moreover, we assume that the competition for food between the juveniles and the mature frogs is negligible, since the mature frogs feed on larger insects. Otherwise, a term of the form −γsM M , for r, s = e, m, l, has to be added to the equations that describe the dynamics of Je, Jm, and Jl. The compartment of mature frogs is assumed to con- tain a homogeneous population the growth of which is governed by the equation dM dt = αeJe + αmJm + αlJl − µM M, where µM = µM (t, M ) = µ1M + µ2M (t)M (t) is the mortality rate consisting of density independent term µ1M , and density dependent term µ2M . Collecting the equations and conditions above yields the following model consisting of five impulsive dif- ferential equations for the larvae, juvenile and mature populations. The model for the dynamics of the Wood Frog popu- lation is: Find five functions: (L(t), Je(t), Jm(t), Jl(t), M (t)), for 0 ≤ t ≤ T , such that, for tk = 365k, k = 0, 1, 2, . . . , dL(t) dt = −(µ1L + µ2L(t)L(t))L(t), t 6= tk, tk + τe, tk + τm, tk + τl, (1) L(0) = L0, (2) L(tk) = σM (tk), k 6= 0, (3) L(tk + τe + 0) = (1 − δe)L(tk + τe − 0), (4) L(tk + τm + 0) = (1 − δm)L(tk + τm − 0), (5) L(tk + τl + 0) = 0, (6) dJe(t) dt = −(µ1e + µ2e(t)Je(t))Je(t) − αeJe(t) −βemJm(t)Je(t) − β e l Jl(t)Je(t), 0 < t 6= tk + τe, (7) dJm(t) dt = −(µ1m + µ2m(t)Jm(t))Jm(t) − αmJm(t) −βml Jl(t)Jm(t) − β m e Je(t)Jm(t), 0 < t 6= tk + τm, (8) dJl(t) dt = −(µ1l + µ2l(t)Jl(t))Jl(t) − αlJl(t) −βleJe(t)Jl(t) − β l mJm(t)Jl(t), 0 < t 6= tk + τl, (9) Je(tk + τe + 0) = Je(tk + τe − 0) +δeL(tk + τe − 0), (10) Jm(tk + τm + 0) = Jm(tk + τm − 0) +δmL(tk + τm − 0), (11) Jl(tk + τl + 0) = Jl(tk + τl − 0) +L(tkτl − 0), (12) dM (t) dt = αeJe(t) + αmJm(t) + αlJl(t) −(µ1M + µ2M M (t))M (t), (13) M (0) = M0, (14) Je(0) = Jm(0) = JL(0) = 0. (15) Here, L0 is the number of larvae and M0 is the number of adult frogs at t = 0, i.e., at the beginning of the first year (k = 0). At that time there are no juveniles, (15). Biomath 1 (2012), 1209032, http://dx.doi.org/10.11145/j.biomath.2012.09.032 Page 3 of 6 http://dx.doi.org/10.11145/j.biomath.2012.09.032 N. Al-Asuoad et al., Model and Simulations of a Wood Frog Population Fig. 2. Larvae vs. t for 41 years III. SIMULATIONS An algorithm for the numerical solutions of the model was constructed and implemented in Maple, using the numerical solver dsolve. The main issue in designing the algorithm was the need to solve the equations be- tween the various times of impulse or transfer, and over different time intervals in each year different equations were solved. In particular, the equation for the larvae, (1), was solved in year k only in the intervals 356k < t < 365k + 75, 356k + 75 < t < 365k + 82, and 356k + 82 < t < 365k + 90, and then L(t) = 0 for 356k + 90 < t < 365(k + 1), that is until the first day of the following year. The various input data were either taken or estimated from [3], or chosen reasonably, and taken as follows. µL1 = 2.8 10 −3, µL2 = 1 10 −7, µe1 = 6 10 −4, µe2 = 2 10 −9, µm1 = 6 10 −3, µm2 = 3.33 10 −9, µl1 = 6 10 −3, µl2 = 1 10 −8, µM 1 = 3.5 10 −3, µM 2 = 1.6 10 −8; αe = 6 10 −4, αm = 5 10 −4, αl = 4 10 −4, βe = 6 10 −6, βm = 5 10 −6, βl = 4 10 −6, δe = δm = δl = 0.5, τl = 75, τm = 82, τl = 90. The figures depict a typical run of 41 years, starting with L0 = 540, 000 larvae, no juveniles, and M0 = 3000 mature frogs. The number of eggs per mature female was Fig. 3. Matures vs. t for 41 years Fig. 4. Early juveniles vs. t for 41 years σ = 600 (a characteristic of Wood Frogs ([3])), and a third of the mature population was females. Under these initial conditions, there is a large drop in year k = 1, and then the populations grow steadily, and in longer simulations (not presented here) they level off to what seems to be steady oscillations. The yearly oscillations are of interest since they cannot be observes directly. To study the effects of a year with harsh conditions, the larvae population was set to be zero in the year k = 21. The larvae population in Fig. 2 is set to zero at day 90 of each year, since all the larvae leave the pond by then. Then, on the first day of the next year a batch of Biomath 1 (2012), 1209032, http://dx.doi.org/10.11145/j.biomath.2012.09.032 Page 4 of 6 http://dx.doi.org/10.11145/j.biomath.2012.09.032 N. Al-Asuoad et al., Model and Simulations of a Wood Frog Population Fig. 5. Middle juveniles vs. t for 41 years Fig. 6. Late juveniles vs. t for 41 years larvae that is proportional to the mature population, (3), appears in the pond from the fertilized eggs. Similarly, on days 75, 82 and 90 portions of the early, middle and late juveniles, respectively, move from the pond to join the juveniles from the previous year. These impulses cause the solutions in Figs. 2, 4–6 to be discontinuous. It is seen that there is leveling or stabilization of the populations due to the density dependent mortality rates in the equations. Since in the year 21 there were no larvae, the other four populations clearly had considerable drops. Never- theless, the whole population recovered over the follow- ing 10 years, reaching a similar behavior as before the drop. Unfortunately, the weather conditions in Michigan this year were such that the larvae in the pond were wiped out completely, and this result allows us to hope that next year the population will begin to recover. The trends in the behavior of the system seem to be similar to what was observed in [3], however the details were not observed, and a number of the coefficients were chosen ‘reasonably.’ Finally, there is strong indication that periodic solutions are possible. IV. CONCLUSION The paper presents a new compartmental model, using impulsive ODEs, for the dynamics of a population of Wood Frogs, based on the field observations of Berven [3]. Then, it depicts computational results for the devel- opment of the larvae, juveniles (early, middle, and late) and the mature frogs populations. Under the given choice of the parameters the populations grow and stabilize in about 20 years. The introduction of environmental adverse conditions that wiped out the larvae population in year 21 show that it takes about 10 years for the population to fully recover. The next stages in this research will be to validate the model by comparing the relevant predictions to the observations in [3]. Once we have confidence in the model, we plan to use it to assess the possible behavior of the population when the environmental conditions are adverse as a result of bad weather. We plan to use the model and the numerical simulations to study various possible future scenarios for the population. The well-posedness of the model, its analysis and stability will be described elsewhere. Moreover, is seems from the numerical results, that the model has periodic solutions, and it is of interest to establish this mathema- tically. ACKNOWLEDGMENT The authors would like to thank the School of Natural Resources and Environment at the University of Michi- gan for access to the research site, and the anonymous referees for some useful suggestions. REFERENCES [1] L.J.S. Allen, An Introduction to Mathematical Biology, Pearson Prentice-Hall, 2007. [2] D. Bainov, P.S. Simeonov, Impulsive Differential Equations: Periodic Solutions and Applications, CRC Press, 1993. [3] K.A. Berven, Density Dependence in the Terrestrial Stage of Wood Frogs: Evidence from a 21-Year Population Study, Copeia, 2009 (2), 328–338 (2009). http://dx.doi.org/10.1643/CH-08-052 Biomath 1 (2012), 1209032, http://dx.doi.org/10.11145/j.biomath.2012.09.032 Page 5 of 6 http://dx.doi.org/10.1643/CH-08-052 http://dx.doi.org/10.11145/j.biomath.2012.09.032 N. Al-Asuoad et al., Model and Simulations of a Wood Frog Population [4] V. Lakshmikantham, D. Bainov, P.S. Simeonov, Theory of Im- pulsive Differential Equations, World Scientific, 1989. http://dx.doi.org/10.1142/0906 [5] H.W. Heathcote, The mathematics of infectious disease, SIAM Review 42 (4), 599–653 (2000). http://dx.doi.org/10.1137/S0036144500371907 [6] H.R. Thieme, Mathematics in Population Biology, Princeton University Press, Princeton, 2003. [7] J.R. Vonesh and O. De la Cruz, Complex life cycles and den- sity dependence: assessing the contribution of egg mortality to amphibian declines, Oecologia 2002 133, 325–333 (2002). Biomath 1 (2012), 1209032, http://dx.doi.org/10.11145/j.biomath.2012.09.032 Page 6 of 6 http://dx.doi.org/10.1142/0906 http://dx.doi.org/10.1137/S0036144500371907 http://dx.doi.org/10.11145/j.biomath.2012.09.032 Introduction The Model Simulations Conclusion References