Modelling of a One-sex Age-structured Population Dynamics with Child Care

The Sharpe-Lotka-Mckendrick-von Foerster one-sex population model and Fredrickson-Hoppensteadt-Staroverov two-sex population one are well known in mathematical biology. But they do not describe dynamics of populations with child care. In recent years some models were proposed to describe dynamics of the wild population with child care. Some of them are based on the notion of the density of offsprings under maternal (or parental) care. However, such models do not ensure the fact that offsprings under maternal (or parental) care move together with their mothers (or both parents). In recent years to solve this problem, some models of a sex-age-structured population, based on the discrete set of newborns, were proposed and examined analytically. Numerical schemes for solving of a one-sex age-structured population model with and without spatial dispersal taking into account a discrete set of offsprings and child care are proposed and results are discussed in this paper. The model consists of partial integro-differential equations subject to conditions of the integral type. Numerical experiments exhibit the stability of the separable solutions to these models.


Introduction
Many species of wild animals care for their offsprings.This phenomenon is natural for many species of mammals and birds and forms the main difference between the behavior of the population taking child care and that without maternal (or parental) duties [1,8].But child care for every species is different.Offsprings of mammals and birds spend some time with their mother or both parents, while young offsprings of some species of fishes, reptilia, and amphibia are left to the own fate.Mammals and birds feed, warm, and defend their young offsprings from enemies.If one of these native duties is not realized, young offsprings die and the population vanishes.For many species of mammals, (e.g., bear (Thalarctos maritimus and Ursus arctos horribilis), whale (Balaenoptera musculus), and panther (Pannthera onca)), only females takes care of their young offsprings.For some species of mammals and birds, (e.g., red fox (Vulpes vulpes), gnawer (Dolichotis patagonium), penquin (Pygoscelis adeliae), heron (Ardea purpurea), falcon (Falco ciolumbarius), and tawny owl (Strix aluco)), both parents take care of their young offsprings.
The Sharpe-Lotka-McKendrick-von Foerster one-sex population model and Fredrickson-Hoppensteadt-Staroverov two-sex population models (see, e.g., [1,8] and references there) are well known in mathematical biology.However, all these models do not treat the child care phenomenon and cannot be used to describe the evolution of the population taking care of its offsprings.Thus, these models can only be applied to describe the dynamics of populations which do not take care of their young offsprings, e.g., some species of fishes, reptilia, and amphibia.Since the behavior of the population taking care of its young offsprings essentially differs from that of the population without child care, the models mentioned above have to be essentially modified.
In [2][3][4][5], four age-structured population dynamics models with child care are proposed, two for one-sex and the other two for two-sex population.The main requirement in these papers is that all offsprings under maternal (or parental) care die if their mother (or any of their parents) dies.In [6], a model for two sex population with strong maternal and weak paternal care is given.The care of this type means that all young offsprings die if their mother dies, but they are protected from the inevitable death in the case of death of their father.All these models are based on the notion of the density of young offsprings.However, such models do not ensure the fact that offsprings under maternal (or parental) care move together with their mothers (or both parents).
To solve this problem, some models of a sex-age-structured population, based on the discrete set of newborns, were proposed.In [1,7,8] two-sex with temporal pairs and one-sex population models, based on a discrete set of offsprings, are proposed.In [9], we discussed numerical results of model [2] for two-sex population with temporal pairs and child care.The present paper is devoted to the numerical investigation of the models studied in [1,8].
The paper is organized as follows.In Sections 3 and 4, the nondispersing population model and the model with spatial diffusion are given.In Section 5, we give the numerical schemes and discuss the numerical solution of models given in Sections 3 and 4. Remarks in Section 6 conclude the paper.

Notation
We use the notation of papers [1,8].R m the Euclidean space of dimension m with x = (x 1 , . . ., x m ), κ the diffusion modulus, (0, T ), (T 1 , T 3 ), the child care and reproductive age intervals, respectively, u(t, τ 1 , x) the age-space-density of individuals aged τ 1 at time t at the position x who are of juvenile (τ 1 ∈ (T, T 1 )), single (τ 1 ∈ (T 1 , T the jump discontinuity of u at the point the minimal age of an individual finishing care of offsprings of the first generation, T 4 = T 3 + T the maximal age of an individual finishing care of offsprings of the last generation, In what follows, κ, T, T 1 , and T 3 are assumed to be positive constants.In the case of nondispersing populations, all functions u, u k , ν, ν k , ν ks , α k , u 0 , and u k0 do not depend on the spatial position x.

The nondispersing population model
In this section, we give the nondispersing age-structured population model [1,8] with stationary vital rates.This model involves the environmental pressure, ρ(N ), by the death rates of juvenile an adult individuals depending on the sum, N, of their spatial densities.It is assumed that young offsprings are subject to natural mortality and are protected from density related increases of mortality dependent on N directly.At age τ 1 = T all young offsprings go to the juvenile group and at age τ 1 = T 1 all juveniles become adult individuals.The model consists of the equations subject to the conditions Here ∂ t and ∂ τ k denote partial derivatives, while n is the biologically possible maximal number of newborns of the same generation produced by an individual.The first term on the right-hand side in equation ( 1) means the part of individuals who produce offsprings, the second and third terms describe the part of individuals whose all young offsprings die and who finish child care, respectively.The transition term k−1 s=0 ν ks u k on the left-hand side in equation ( 2) describes the part of individuals aged τ 1 at time t who take child care of k young offsprings and whose at least one young offspring dies.Similarly, the term on the right-hand side in this equation describes a part of individuals aged τ 1 at time t who take care of more than k, 1 ≤ k ≤ n − 1, young offsprings aged τ 2 whose number after the death of the other offsprings becomes equal to k.The condition [u| τ1=τ ] = 0, τ = T 1 , T 2 , T 3 , and T 4 means that the function u must be continuous at the point, τ 1 = τ, of the discontinuity of the right-hand side of equation (1).
Given functions ν, ν k , ν ks , α k , u 0 , and u k0 and the unknown ones u and u k are to be positive.The positive constants T and T s are also to be given.
We use the following compatibility conditions: Using the following transformations, model ( 1)-( 5) can be decomposed into the problem for U and U k , and the equations for f and N, The existence and uniqueness theorem for problems ( 7)-( 9), (5), and ( 10), ( 11) is proved and, in the case of stationary vital rates, the large time behavior of solution ( 6) is given in [1,8].
In the case of time-independent vital rates ν, ν k , ν ks , and α k , the existence of separable solutions of the form with unique real λ is also proved in [1,8].Here U > 0 is an arbitrary constant, while the constant λ and positive functions v λ and v λ k are to be determined.Functions v λ and v λ k satisfy the problems where the prime denotes the differentiation of v λ with respect to age τ 1 , and while the constant λ is a root of the characteristic equation Set where v 0 and v 0 k satisfy equations ( 14) and (15) for λ = 0. From equations ( 16) and (17) we get the characteristic equation for λ, Under minimal restrictions on the time-independent vital rates ν, ν k , ν ks , and α k and initial functions u 0 and u 0 k , it is proved in [1,8] that the solution of problem ( 7)-( 9) and ( 5) for u 0 and u 0 k of the general type tends, as t → ∞, to separable solution (13) with unique U and exponent λ satisfying equation (18).Therefore, λ can also be determined by the formula Equation ( 10) for the large time can be written in the form where β 0 is a positive constant and ρ has a positive derivative.It is well known that, for Results of numerical investigation of problem ( 5) and ( 7)-( 11) are discussed in Section 5.

A population model with spatial diffusion
In this section, we consider the population dynamics including the spatial diffusion in the interval (0, 1) with the extremely inhospitable points x = 0 and x = 1.We use the model analyzed in [1,8], where the constant κ is the diffusion modulus.We assume that ν, ν k , α k , and ν sk do not depend on t and x and use the following compatibility conditions: In the case of product initial functions, with positive U 0 and U k0 , and f 0 (0) = f 0 (1) = 0, f 0 (x) > 0 in (0; 1), by the transform we split (20)-(23) into problem ( 7)-( 9) and ( 5) for U and U k and the problem for f, where β(t) and λ are defined by equations ( 11) and (19), respectively.Numerical results are discussed in Section 5.
The general type of initial functions does not allow to split problem (20)-( 23) into two problems.In this case we solved system (20)-( 23) directly.
In all calculations, excluding Fig. 16, we use the vital rates and initial functions for nondispersing population model and for the population with the spatial diffusion.Here Using compatibility conditions (23) 1 , we determine the function Note that formula (30) with p 1 = 0 determines the constant β 2 in equations (28).The positive constants , α k2 , τ 0 < T, and ρ 0 remain free.
It is easy to see that (29) represent the sum of two product functions, and is the particular case of functions (4.2.1) considered in [1,8].Now we describe numerical schemes for both (1)-( 5) and ( 20)-( 23) models.
Brief description of the procedure for solving system (1)-( 5).We solve this system using the Chiu [10] (see, also, [11]) scheme applied to the Gurtin-MacCamy model.To do this, we write model ( 1)-( 5) on the characteristic lines and then replace it by a discrete system using the same step for ages and time.At each time level t s+1 we use the following iteration procedure: 1. Determine N (t 0 ), t 0 = 0; 2. Determine N (0) (t s ) by the formula N (0) (t s ) = N (t s−1 ), s = 1, 2, . . .using known u and u k for t s−1 ;

Determine u
(n) k for t s except u (n) k τ2=0 , k = 1, 2, 3; 4. Determine u (n) for t s by the explicit scheme.To do this, we use the rectangle formula to calculate the second term of the right-hand side of equation (1);

Determine u
(n) k τ2=0 , k = 1, 2, 3, for t s ; 6. Determine N (n) using found u (n) and u n k for t s , and repeat steps 3, 4, 5, 6 until inequality max Brief description of the procedure for solving system (7)-( 12) and ( 5) is similar to that used to solve model ( 1)-( 5).Here, we do not use iteration procedure.We apply steps 2, 3, and 4 from the procedure above and use the explicit scheme to determine f from equation (10).Then, by equation ( 6), we determine u and u k .
To solve system (20)-(23), we use the Crank-Nicholson scheme and the iteration process: 1. Determine N (t 0 , x), t 0 = 0; 2. Determine N (0) (t s , x) by the formula N (0) (t s , x) = N (t s−1 , x), s = 1, 2, . . .using known u and u k for t s−1 ; 3. Determine u 4. Determine u (n) for t s by the explicit scheme.To do this, we use the rectangle formula to calculate the second term of the right-hand side of equation ( 1); 6. Determine N (n) using found u (n) and u n k for t s , and repeat steps 3, 4, 5, 6 until inequality max Procedure for solving of system (7)-( 9), (5), and (24)-( 26).We determine U and U k , k = 1, 2, 3, in the same way as in the procedure for solving of system ( 7)-( 12), (5).Then determine f from problem (26) using the iteration procedure with the Crank-Nicholson scheme.At last, by equation ( 6), we determine u and u k .
Brief description of the procedure for solving equation (18).
3. Determine v 0 (τ 1 ), T 1 < τ 1 < T 4 .To do this, we apply the rectangle formula to calculate the second term on the right-hand side of equation ( 14) and then use the explicit scheme.
and the constants given in Table 1.Here

Define
where Here v(t, τ 2 , x), V (t, x), and z k (t, x) mean the number of offsprings aged τ 2 at time t at the positions x, the total number of offsprings at time t at the positions x, and the number of individuals taking care of k offsprings at time t at the location x.
Figs. 1-7 represent the solution to problem (1)-( 5).The influence of death rates, environmental pressure, and the rate α 3 on the behavior of densities u, u k , v, and N + V is investigated.
Fig. 1 describes the behavior of u 1 , u 2 , and u 3 for t=10 and the vital parameters from Set (31) and Table 1.We observe the inequality u 1 < u 2 < u 3 for α 11 < α 21 < α 31 and fixed the other parameters.
Fig. 5 illustrates the total population N + V versus t for µ k = µ k1 = µ k2 = 0.0001, 0.001, 0.01 with k = 1, 2, 3, and the other parameters from (31) and Table 1.We see that N + V decreases as deaths rates increase.It also possesses a local maximum in time.
Fig. 6 describes the behavior of N + V versus t for α 31 = 0.14, 0.12, 0.1, and the other parameters from (31) and Table 1.Function N + V decreases with decreasing α 31 and possesses a local maximum in time.
Fig. 7 exhibits the behavior of V with respect to t for T 1 = 1.9, 2, 2.1, and the other parameters from (31) and Table 1.Function V decreases as T 1 increases.
Fig. 11 illustrates the influence of ρ 0 on the behavior of N (t, x).The other parameters are used from (31) and Table 1.
Figs. 12 and 13 describe the influence of κ on the behavior of N (t, x) for ρ 0 = 0.001 and ρ 0 = 0, respectively.The other parameters are used from (31) and Table 1.N decreases as as ρ or κ increase.
Fig. 16 illustrates the influence of initial functions (29) and (32) on the behavior of z 1 for the parameters from (31) and Table 1.
We also solved model ( 20)-( 23) with initial functions of the form where Note that (32) do not belong to the class of the product initial functions.It is shown, in [1,8], that solution of model ( 20)-( 23) with ρ = 0 and initial functions of type ( 24) and (29) tend to a separable solution with exponent Λ = λ − κπ 2 , where λ is defined by formula (19).Numerical calculations confirm this result and also show that the solutions of model ( 20)-( 23) with ρ = 0 and initial functions (32), for large time, possess the exponent Λ = Λ.This result can also be shown analytically using the solution expansion by the eigenfunctions of the Dirichlet problem to the Laplace operator.and parameters from (31) and Table 1.  1.

Concluding remarks
A one-sex age-structured population dynamics model [1] (see also [8]) is examined numerically herein, taking into account an environmental pressure, a discrete set of offsprings, and child care.Both the nondispersing population model and that describing the spatial diffusion are considered.Numerical schemes are based on the method of the characteristic lines.Numerical results are exhibited in the graphs and illustrate the behavior of the solution to these models depending on the various values of the parameters determining vital functions.Nondispersing population model ( 1)-( 5) is solved directly and by using decomposed problem ( 7)-( 12) and (5).
Numerical experiments show that, in the case of zeroth environmental pressure and initial functions of the general type (at least for functions (32)), the solution of model (20)-(23) tends to a separable solution with an exponent Λ = λ − κπ 2 as time tends to infinity, where λ is defined by formula (19).This result can also be proved analytically.For initial functions of type ( 24) and (29) this result is proved analytically in [1,8].
3 )), or postreproductive (τ 1 > T 3 ) age, u k (t, τ 1 , τ 2 , x) the age-space-density of individuals aged τ 1 at time t at the position x who take care of their k offsprings aged τ 2 at the same time,ν(t, τ 1 , x)the natural death rate of individuals aged τ 1 at time t at the position x who are of juvenile or adult age, ν k (t, τ 1 , τ 2 , x) the natural death rate of individuals aged τ 1 at time t at the position x who take care of their k offsprings aged τ 2 , ν ks (t, τ 1 , τ 2 , x) the natural death rate of k-s young offsprings aged τ 2 at time t at the position x whose mother is aged τ 1 at the same time, α k (t, τ 1 , x)dt the probability to produce k offsprings in the time interval [t, t + dt] at the location x for an individual aged τ 1 , N sum of spatial densities of juvenile and adult individuals, ρ(N )the death rate conditioned by ecological causes (overcrowding of the population), ρ(0