# On a model of American cutaneous leishmaniasis.

1. Introduction

American Cutaneous Leishmaniasis (ACL) is a resurgent disease that is caused by Leishmania parasites [5], and transmitted by female phlebotomine sandflies ("vectors") from reservoir hosts. This is why a reservoir host is considered to be the main factor in changing the dynamics of the disease. It is a parasitic disease occurring throughout the Americas from Texas to Argentina. Mathematical models of general dynamical systems have proved to be essential for the local and global study of numerous phenomena not only in biology but also in physics, astronomy, economics, etc. ACL unfortunately has not been studied thoroughly from the mathematical point of view when compared to other diseases such as malaria [1].

A study of ACL as a four-dimensional discrete dynamical system is performed in [3]. In [2] threshold conditions for the establishment of the disease are calculated, in the framework of a three-dimensional continuous system, which includes a study of stability of the equilibria.

In this paper we study three and four-dimensional continuous models of ACL, further exploring on the existence of periodic solutions, connecting orbits between equilibria, and other special solutions. We also perform smooth continuation of solutions, that is, we allow some parameters to vary continuously and compute the corresponding new solutions, tracking on the way for possible bifurcations and other special solutions. This naturally includes a study of the stability of such solutions.

2. Setting of the Problem

Let us denote the population of infected humans, donkeys, dogs and vectors as H = H(t), R = R(t), P = P(t), V = V(t) respectively. Here donkeys and dogs are considered to be potential reservoir hosts, and humans are only incidental hosts (they are sink but not source of infection). The main assumption for the model is that the net rate of spread of the infection is proportional to the product of the density of susceptible individuals and the density of infectious individuals.

The parameters of the model (all positive) are [[beta].sub.H], [[beta].sub.R], [[beta].sub.P]: rate of infection per individual in vector-human, vector-donkey and vector-dog encounters respectively; [[gamma].sub.H], [[gamma].sub.R], [[gamma].sub.P]: rate of recovery per unit time of infected humans, donkeys and dogs respectively, and [mu]: rate of mortality of vectors per unit time. Depending whether one considers one or two species as reservoir hosts, the corresponding models to be considered are

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2.1)

and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2.2)

where I in (2.1) can be R or P, depending if one considers donkeys or dogs as the unique reservoir; A, B, C, and D represent the total population of humans, donkeys, vectors and dogs respectively in the ecosystem under consideration. Due to the nature of the equations, we can consider three-dimensional systems as particular cases. By instance, if dogs are not part of the species under study, then we have as in [2]:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2.3)

To study possible existence of periodic solutions of the systems under consideration, we need to introduce some definitions and results. For details, see [8], [13]. For simplicity, the results are mainly given with (2.3) in mind, but they can easily be extended to the four-dimensional systems.

For x, y [member of] [R.sup.n], we write x [less than or equal to] y [??] [x.sub.i] [less than or equal to] [y.sub.i], i = 1, 2,...,n; x < y [??] x [less than or equal to] y and for some 1 [less than or equal to] i [less than or equal to] n, [x.sub.i] < [y.sub.i]. Finally, we denote x [much less than] y [??] [x.sub.i] < [y.sub.i], for all i = 1, 2,...,n.

Let K = [R.sup.n.sub.+], that is, the cone of nonnegative vectors, and consider the dynamical system

[??] = f(x), (2.4)

where f: D [subset] [R.sup.n] [right arrow] [R.sup.n], and D is an open set.

Definition 2.1. The system (2.4) is called K-cooperative (or simply, cooperative) if f is differentiable and

[partial derivative][f.sub.i] (x) [greater than or equal to] 0, for all i [not equal to] j, and it is called monotone if

[x.sub.0] [less than or equal to] [x.sub.1] [??] x([t.sub.0], [x.sub.0]) [less than or equal to] x([t.sub.0], [x.sub.1]).

Definition 2.2. A set D in [R.sup.n] is called p-convex if for every x, y [member of] D satisfying x [less than or equal to] y, we have tx + (1 - t)y [member of] D for all t [member of][0, 1].

Proposition 2.3. Let D be p-convex. If (2.4) is K-cooperative, then it is monotone.

Proposition 2.4. If (2.4) is monotone, then it does not have any attracting periodic orbits.

3. Dynamics of the Model

Periodic solutions. We want to apply the results of the previous section to system (2.3). We start with the following

Theorem 3.1. The ACL dynamical system (2.3) does not have any attracting periodic solution.

Proof. It is clear that the set D = [0,A] x [0,B] x [0,e] is convex and therefore it is also p-convex with the partial order defined before. We also let x = (H, R, V) and f(x) = ([f.sub.1](x), [f.sub.2](x), [f.sub.3](x)) with

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3.1)

so that (2.3) has the general form (2.4). Then, for system (2.3) we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

so that [[partial derivative][f.sub.i]/[partial derivative][x.sub.j]] [greater than or equal to] 0 for i [not equal to] j, and i, j = 1, 2, 3; that is, system (2.3) is cooperative in D. Therefore, an immediate consequence from Propositions 2.3 and 2.4 is that system (2.3) cannot have attracting periodic orbits. Similar conclusions can be arrived at when changing the reservoir hosts, and when dealing with system (2.2).

One approach to study the nonexistence of any periodic solutions of (2.3), or more generally, of (2.2) is to apply a generalization of Bendixon's criterion. The idea is to compute the Jacobian J of say, (2.2) and define the symmetric matrix

S = [1/2] ([J.sup.T] + J).

Theorem 3.3 in [11] says that if S has eigenvalues [[lambda].sub.1] [greater than or equal to] [[lambda].sub.2] [greater than or equal to] ... [greater than or equal to] [[lambda].sub.n] such that [[lambda].sub.1] + [[lambda].sub.2] < 0, then no periodic solution can exist for system (2.2). Numerical calculations indicate that for any

R [greater than or equal to] 2, P [greater than or equal to] 0, V [greater than or equal to] 3255 and H [member of] [0, 124],

it holds indeed that [[lambda].sub.1] + [[lambda].sub.2] < 0, and therefore there are no periodic solutions on those intervals.

Explicit expressions for the eigenvalues of S in terms of H, R, P and V are too complicated to be useful, especially if we consider the full system (2.2). However if in particular we assume that the entire population of vectors, dogs and donkeys are infected, then we can arrive to a conclusion about periodic solutions.

Theorem 3.2. The system (2.2) has no periodic solutions when entire population of vectors, dogs and donkeys are infected.

Proof. The matrix S takes the form

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

In such a case, the eigenvalues of S are (in order of magnitude)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where a = -[[beta].sub.H]V - [[gamma].sub.H], b = [1/2] [[beta].sub.H](A - H), g = -[[beta].sub.P]P - [[beta].sub.R]R - [mu].

One can readily verify that we always have [[lambda].sub.1] + [[lambda].sub.2] < 0, and therefore there are no periodic solutions to such system.

Invariance. It is possible to explicitly define some sets that are invariant with respect to the flow of the dynamical system. To this end, let <[tau] stand for one of the relations previously defined <, [less than or equal to], [much less than] and let

[P.sub.+] = {x [member of] [OMEGA]: f(x) [>.sub.[tau]] 0}, [P.sub.-] = {x[member of]e [OMEGA]: f(x) [<.sub.[tau]] 0}.

For (2.3) it is readily checked that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The following proposition is a direct consequence of system (2.3) being cooperative. See [13].

Proposition 3.3.

(a) [P.sub.+] and [P.sub.-] are positively invariant.

(b) If we denote [[PHI].sub.t]([x.sub.0]) the solution of (2.3) with [[PHI].sub.0]([x.sub.0]) = [x.sub.0] and [x.sub.0] [member of] [P.sub.+] (resp. [x.sub.0] [member of] [P.sub.-]) then [[PHI].sub.t]([x.sub.0]) is non-decreasing (resp. non-increasing).

Remark. From the definition of invariance, we can say e.g. that the first inequality in the set [P.sub.+] gives an explicit upper bound on the number of infected humans that will result when starting with values of H, R and V that satisfy the given inequalities.

We now move in the direction of the spectral properties of Jacobian matrix associated to (2.3) specifically.

Theorem 3.4. Consider the system (2.3). Then,

(a) the origin is a saddle.

(b) The nontrivial equilibrium is a saddle if [alpha] < [[beta].sup.2.sub.R]BC, and depending on the sign of certain discriminant, it is either a stable node or a stable focus if [alpha] > [[beta].sup.2.sub.R]BC,

where [alpha] = [[beta].sub.R]([[beta].sub.R]BV + [[beta].sub.R]CR + [[gamma].sub.R]R + [mu]V) + [mu][[gamma].sub.R].

Proof. The critical points of (2.3) are, the origin [x.sub.0] = (0, 0, 0) and a non trivial one

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where we require that

[[beta].sup.2.sub.R]BC/[mu] > [[gamma].sub.R], (3.2)

so that the equilibrium [x.sub.1] lies in the positive octant. The general Jacobian of (2.3) is

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3.3)

At the origin, this reduces to

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and its eigenvalues can be computed directly. They are

[[lambda].sub.1] = -[[gamma].sub.H], [[lambda].sub.2,3] = -([[lambda].sub.R] + [mu]) [+ or -] [square root of [([[gamma].sub.R] + [mu]).sup.2] + 4([[beta].sup.2.sub.R]BC - [mu][[gamma].sub.R])/R]

and it is clear that with the condition (3.2), the origin is a saddle point, and therefore unstable. (Also observe that if had [[beta].sup.2.sub.R]BC/[mu] = [[gamma].sub.R] then the origin would be nonhyperbolic and if [[[beta].sup.2.sub.R]BC/[mu]] < [[gamma].sub.R] then it would be stable).

To study the eigenvalues of the Jacobian (3.3) at the nontrivial equilibrium [x.sub.1], we see that one of the eigenvalues is [[lambda].sub.1] = -[[beta].sub.H]V - [[gamma].sub.H] < 0. The other two eigenvalues come from the 2 x 2 lower right block. If we denote with 8 the determinant of that block and with [tau] its trace, then the two eigenvalues are given by

[lambda] = t [+ or -] [square root of [[tau].sup.2] - 4[delta]]/2, (3.4)

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Observe that [delta] can be written as:

[delta] = ([[beta].sup.2.sub.R]BV + [[beta].sup.2.sub.R]CR + [[beta].sub.R][[gamma].sub.R]R + [mu][[beta].sub.R]V + [mu][[gamma].sub.R]) - [[beta].sup.2.sub.R]BC = [alpha] - [[beta].sup.2.sub.R]BC.

Thus, if [alpha]< [[beta].sup.2.sub.R]BC, then [delta] < 0 and therefore there are two eigenvalues of opposite sign and [x.sub.1] is saddle. Similarly, if [alpha] > [[beta].sup.2.sub.R]BC, then [delta] > 0, and since we know that [tau] < 0, we can conclude that [x.sub.1] is a stable node if [[tau].sup.2] - 4[delta] [greater than or equal to] 0 or a stable focus if [[tau].sup.2] - 4[delta]< 0.

Several observations are in order:

(a) Biologically, the inequality (3.2) says that the ratio of the rate of infection of the main reservoir hosts (donkeys) to the rate of mortality of the vectors (scaled by the population product BC) will have to be greater than the rate of recovery of the infected reservoir host, so that there is some infected population. This is in fact a threshold condition for infection [2].

(b) Under the threshold condition (3.2), the origin is unstable. This means that a state with no infected individuals can be easily perturbed into a state with some infected individuals.

(c) Although it is not easy to determine the sign of [[tau].sup.2] - 4[delta] in general, for the application in mind, the parameter [mu] is many times larger that the other parameters (see Section 4). This makes the term [[tau].sup.2] much larger that 4[delta], and therefore [[tau].sup.2] > 4[delta].

(d) The inequalities above depend much on the quantity [[beta].sup.2.sub.R]BC, that is, the rate of infection in vector-donkey encounters, the total population of donkeys and the total population of vectors, in other words, the very carriers of the disease.

4. Numerical Dynamics

In this section we perform a numerical study of the system (2.1), where there is only one species (dogs or donkeys) considered as reservoir, and of the extended system (2.2), where both dogs and donkeys at the same time are considered to be reservoirs. Among other things we study sability of equilibria and periodic orbits, location of branch points, and special solutions (such as connecting orbits). We also perform pseudo arc- length numerical continuation, following the theoretical background and numerical algorithms from [4] and [12]. This amounts to an extension of the numerical analysis of the local dynamics sudied in [2] to the global numerical analysis for systems in three and four variables.

[FIGURE 1 OMITTED]

[FIGURE 2 OMITTED]

To begin with, for the values [3] [[beta].sub.H] = 0.000982, [[beta].sub.R] = 0.0000527, [[beta].sub.P] = 0.000085, [[gamma].sub.H] = 0.024, [[gamma].sub.R] = 0.0007, [[gamma].sub.P] = 0.003, [mu] = 0.42, A = 124, B = 28, C = 4807, D = 43,

we solve (2.1) with I = R, and we also solve (2.2). The initial values are H(0) = 0, R(0) = 6, P(0) = 4, V(0) = 0. In Figures 1 and 2 we show the numerical results for the four species as functions of time. As expected, all species are nondecreasing functions of time t, which means that in the absence of prophylaxis and vaccines, the number of infected individuals will increase as time evolves. Of special interest is the evolution of the number of infected humans H(t).

For a closer look at the dynamic behavior of the systems, let us drop for a moment the P component in system (2.1), which leaves us the system

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.1)

We want to perform numerical continuation on the parameter [[beta].sub.H]. That is, starting from a solution (in this case the trivial one) of (4.1), we let [[beta].sub.H] vary and generate corresponding solutions of the system, and at the same time we track for possible special points and solutions. It is clear in this case that varying [[beta].sub.H] only affects H whereas R and V remain fixed. Thus, we plot only H against [[beta].sub.H]. We can make continuation for large intervals of [[beta].sub.H], but of course values larger than say, [[beta].sub.H] = 0.2 would be unrealistic for the given application. No special solutions (periodic orbits, etc.) nor any bifurcation point was detected. See Figure 3.

However, when performing continuation with respect to the parameter [[beta].sub.R], a bifurcation (branch point), was detected at [[beta].sub.R] = 4.673665E - 04 for the small values H = 2.591760E - 04,R = 9.549715E - 05,V = 5.108258E - 05. See Figure 4. For these values of the parameters, the inequality (3.2) becomes an equality. Connecting Orbits. Let us now go back to the general model (2.2). If we fix all parameters except [[beta].sub.H], then besides the trivial solution, there are two more equilibria. One of them has some negative components, so that it is not of practical interest. The third equilibrium is given by

Q = (-0.13817 [[beta].sub.H] + 124 [[beta].sup.2.sub.H]/(0.000766193 + [[beta].sub.H])(-0.00111427 + [[beta].sub.H]), 19.6623, 20.2186, 31.3237), (4.2)

[FIGURE 3 OMITTED]

[FIGURE 4 OMITTED]

The general Jacobian of (2.2) is

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.3)

Due to the block triangular form of the Jacobian, we observe that - [[beta].sub.H]V - [[gamma].sub.H] < 0 is one of the eigenvalues and is the only one that depends on [[beta].sub.H]. For [[beta].sub.H] = 0.000982 and the fixed values of the other parameters, the origin is unstable and Q is stable. In fact, the eigenvalues of the Jacobian at the origin are

[[lambda].sub.1] = -0.02400000, [[lambda].sub.2] = -0.42442592, [[lambda].sub.3] = -0.00141438, [[lambda].sub.4] = 0.00214030.

This also tells us the unstable manifold of the origin is one-dimensional. Similary, the eigenvalues of the Jacobian at Q are

[[mu].sub.1] = -0.05475986, [[mu].sub.2] = -0.42489150, [[mu].sub.3] = - 0.00183740, [[mu].sub.4] - 0.00403915,

which indicates that the stable manifold of P is four-dimensional.

We are interested in computing a connecting orbit from the origin to Q. This is of great importance because such a curve explicitly describes the path from non infected population to the stable state of infected individuals. In general, for a dynamical system

[??] = f(x, [lambda]), x(t) [member of] [R.sup.n], [lambda] [member of] [R.sup.p], (4.4)

where f: [R.sup.n] x [R.sup.p] [right arrow] [R.sup.n], and for [M.sub.[+ or - ]]([lambda]) representing equilibria or periodic orbits, a solution x(t, [lambda]), t [member of] R of (4.4) is called a connecting orbit from [M.sub.-]([lambda]) to [M.sub.+]([lambda]) if

dist (x(t, [lambda]), [M.sub.[+ or -]]([lambda])) [right arrow] 0, as t [right arrow] [+ or -][infinity].

According to [4] and [12], the tangent spaces of the corresponding stable and unstable manifolds should span the whole space and the intersection of such manifolds should be transversal. This guarantees the existence of an isolated connecting orbit between the two equilibria. From these conditions, a fundamental relation or balance between the dimensions of the manifolds and the number of parameters is established as

p = [m.sup.u.sub.+] - [m.sup.u.sub.-] - [m.sup.c.sub.-] + 1, (4.5)

where, for the problem considered here, [m.sup.u.sub.+], [m.sup.u.sub.-], [m.sup.c.sub.-] represent the dimensions of the unstable manifold of P and that of the unstable and center manifolds of the origin. So that for this particular case we have p = 0 - 1 - 0 + 1 = 0, so that the connection between the two equilibria can be established without free parameters.

Remark. In general, if (4.5) is not satisfied, then we should either add parameters to the system or add conditions for the parametrization of a manifold of connecting orbits. More precisely, if p < [m.sup.u.sub.+] - [m.sup.u.sub.-] - [m.sup.c.sub.-] + 1, we should add ([m.sup.u.sub.+] - [m.sup.u.sub.-] - [m.sup.c.sub.-] + 1) - p parameters in order to have a well-posed problem, while if p > [m.sup.u.sub.+] - [m.sup.u.sub.-] - [m.sup.c.sub.-] + 1, then we can add p - ([m.sup.u.sub.+] - [m.sup.u.sub.-] - [m.sup.c.sub.-] + 1) constraints to select a unique connecting orbit.

We do not limit ourselves to compute one connecting orbit between the origin and Q but instead we compute a branch of connecting orbits through smooth continuation, for [[beta].sub.H] in the interval [0.00001, 0.05]. See Figures 5 and 6, where we plot the time evolution of the infected humans. Connecting orbits are very important in the global study of dynamical systems, acting as separatrices of solutions in the phase portrait, giving a better understanding of the behavior of nearby solutions.

For system (2.2) we have also performed continuation on the parameter [[beta].sub.H] on the interval [0.000982, 0.1], starting from the equilibrium point (4.2). Numerical continuation for larger values of [[beta].sub.H] is still possible but it is not of practical interest. For the value [[beta].sub.H] = 0.1, the corresponding value of infected humans is already about H = 123. No special solutions or bifurcations are found. For [[beta].sub.R] we have performed continuation on the interval [0.0000527, 0.0037734]. A critical point at the right endpoint of this interval is found to be

U = (121.0175, 27.9947, 41.5222, 991.16723), (4.6)

which is stable. The other nontrivial equilibrium has some negative components. Just as we computed connecting orbits from the origin to P above, connecting orbits from the origin to U can also be computed.

[FIGURE 5 OMITTED]

[FIGURE 6 OMITTED]

[FIGURE 7 OMITTED]

5. Final Remarks

In this work we have studied models of American cutaneous leishmaniasis in the form of continuous dynamical systems, extending research in the same area in the current literature. We have explicitly computed paths that connect different states and we have performed continuation as a problem parameter is allowed to vary continuously. We have argued on the nonexistence of periodic solutions and have provided a study of stability of solutions. A future avenue of research would be to consider perturbations on the system and study their impact on the overall dynamics; this may introduce other interesting mathematical concepts into play such as bifurcations or strange attractors. A ratio-dependent interaction in the sense of [10] between vectors and humans (or between other species) would generate mathematical models of great interest, though more complicated.

References

[1] Awerbuch T., 1994, Evolution of mathematical models of epidemics. Ann. N. Y. Acad. Sci., 740:232-241.

[2] Chaves L.F. and Hernandez M.J., 2004, Mathematical modelling of American cutaneous leishmaniasis: incidental hosts and threshold conditions for infection persistence. Acta Tropica, 5:1-8.

[3] Chaves L., Hernandez M., and Ramos S., 2008, Simulacion de modelos matematicos como herramienta para el estudio de los reservorios de la leishmaniasis cutanea americana, Divulgaciones Matematicas, 16:125-154.

[4] Dieci L. and Rebaza J., 2004, Point-to-point and point-to-periodic connections. BIT Numerical Mathematics, 44:41-62.

[5] Gratz N.G., 1999, Emerging and resurging vector-borne diseases. Ann. Rev. Entomol, 44:51-75.

[6] Miguel de Guzman, 1980, Ecuaciones Diferenciales Ordinarias. Alhambra, Madrid.

[7] Jack Hale, 1989, Ordinary Differential Equations. Malabar.

[8] Hirsch M.W. and Smith H., 2004, Monotone Dynamical Systems, Electronic version.

[9] Killick-Kendrick R. and Ward R.D., 1982, Ecology of Leishmania, Parasitilogy, 82:143-152.

[10] Leard B., Lewis C., and Rebaza J., 2008, Dynamics of ratio-dependent predator-prey models with nonconstant harvesting, Disc. Cont. Dyn. Syst. S, 1:303-315.

[11] Li Y. and Muldowney J., 1993, On Bendixon's criterion, J. Diff. Equations, 106:27-39.

[12] Rebaza J., 2007, Smooth Schur factorizations in the continuation of separatrices, Linear Algebra and its Applications, 421:138-156.

[13] Hal Smith, 1995, Monotone dynamical systems: An introduction to the theory of competitive and cooperative systems, Math. Surv. Mon. AMS, 41.

Jesus Carreno

Departamento de Matematica, Universidad de Oriente, Nucleo de Sucre, Venezuela

E-mail: carreno@sucre.udo.edu.ve

Teodoro Lara

Departamento de Fisica y Matematicas, Universidad de los Andes, N. U. Rafael Rangel, Trujillo, Venezuela

E-mail: tlara@ula.ve

Jorge Rebaza

Department of Mathematics, Missouri State University

E-mail: jrebaza@missouristate.edu

American Cutaneous Leishmaniasis (ACL) is a resurgent disease that is caused by Leishmania parasites [5], and transmitted by female phlebotomine sandflies ("vectors") from reservoir hosts. This is why a reservoir host is considered to be the main factor in changing the dynamics of the disease. It is a parasitic disease occurring throughout the Americas from Texas to Argentina. Mathematical models of general dynamical systems have proved to be essential for the local and global study of numerous phenomena not only in biology but also in physics, astronomy, economics, etc. ACL unfortunately has not been studied thoroughly from the mathematical point of view when compared to other diseases such as malaria [1].

A study of ACL as a four-dimensional discrete dynamical system is performed in [3]. In [2] threshold conditions for the establishment of the disease are calculated, in the framework of a three-dimensional continuous system, which includes a study of stability of the equilibria.

In this paper we study three and four-dimensional continuous models of ACL, further exploring on the existence of periodic solutions, connecting orbits between equilibria, and other special solutions. We also perform smooth continuation of solutions, that is, we allow some parameters to vary continuously and compute the corresponding new solutions, tracking on the way for possible bifurcations and other special solutions. This naturally includes a study of the stability of such solutions.

2. Setting of the Problem

Let us denote the population of infected humans, donkeys, dogs and vectors as H = H(t), R = R(t), P = P(t), V = V(t) respectively. Here donkeys and dogs are considered to be potential reservoir hosts, and humans are only incidental hosts (they are sink but not source of infection). The main assumption for the model is that the net rate of spread of the infection is proportional to the product of the density of susceptible individuals and the density of infectious individuals.

The parameters of the model (all positive) are [[beta].sub.H], [[beta].sub.R], [[beta].sub.P]: rate of infection per individual in vector-human, vector-donkey and vector-dog encounters respectively; [[gamma].sub.H], [[gamma].sub.R], [[gamma].sub.P]: rate of recovery per unit time of infected humans, donkeys and dogs respectively, and [mu]: rate of mortality of vectors per unit time. Depending whether one considers one or two species as reservoir hosts, the corresponding models to be considered are

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2.1)

and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2.2)

where I in (2.1) can be R or P, depending if one considers donkeys or dogs as the unique reservoir; A, B, C, and D represent the total population of humans, donkeys, vectors and dogs respectively in the ecosystem under consideration. Due to the nature of the equations, we can consider three-dimensional systems as particular cases. By instance, if dogs are not part of the species under study, then we have as in [2]:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2.3)

To study possible existence of periodic solutions of the systems under consideration, we need to introduce some definitions and results. For details, see [8], [13]. For simplicity, the results are mainly given with (2.3) in mind, but they can easily be extended to the four-dimensional systems.

For x, y [member of] [R.sup.n], we write x [less than or equal to] y [??] [x.sub.i] [less than or equal to] [y.sub.i], i = 1, 2,...,n; x < y [??] x [less than or equal to] y and for some 1 [less than or equal to] i [less than or equal to] n, [x.sub.i] < [y.sub.i]. Finally, we denote x [much less than] y [??] [x.sub.i] < [y.sub.i], for all i = 1, 2,...,n.

Let K = [R.sup.n.sub.+], that is, the cone of nonnegative vectors, and consider the dynamical system

[??] = f(x), (2.4)

where f: D [subset] [R.sup.n] [right arrow] [R.sup.n], and D is an open set.

Definition 2.1. The system (2.4) is called K-cooperative (or simply, cooperative) if f is differentiable and

[partial derivative][f.sub.i] (x) [greater than or equal to] 0, for all i [not equal to] j, and it is called monotone if

[x.sub.0] [less than or equal to] [x.sub.1] [??] x([t.sub.0], [x.sub.0]) [less than or equal to] x([t.sub.0], [x.sub.1]).

Definition 2.2. A set D in [R.sup.n] is called p-convex if for every x, y [member of] D satisfying x [less than or equal to] y, we have tx + (1 - t)y [member of] D for all t [member of][0, 1].

Proposition 2.3. Let D be p-convex. If (2.4) is K-cooperative, then it is monotone.

Proposition 2.4. If (2.4) is monotone, then it does not have any attracting periodic orbits.

3. Dynamics of the Model

Periodic solutions. We want to apply the results of the previous section to system (2.3). We start with the following

Theorem 3.1. The ACL dynamical system (2.3) does not have any attracting periodic solution.

Proof. It is clear that the set D = [0,A] x [0,B] x [0,e] is convex and therefore it is also p-convex with the partial order defined before. We also let x = (H, R, V) and f(x) = ([f.sub.1](x), [f.sub.2](x), [f.sub.3](x)) with

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3.1)

so that (2.3) has the general form (2.4). Then, for system (2.3) we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

so that [[partial derivative][f.sub.i]/[partial derivative][x.sub.j]] [greater than or equal to] 0 for i [not equal to] j, and i, j = 1, 2, 3; that is, system (2.3) is cooperative in D. Therefore, an immediate consequence from Propositions 2.3 and 2.4 is that system (2.3) cannot have attracting periodic orbits. Similar conclusions can be arrived at when changing the reservoir hosts, and when dealing with system (2.2).

One approach to study the nonexistence of any periodic solutions of (2.3), or more generally, of (2.2) is to apply a generalization of Bendixon's criterion. The idea is to compute the Jacobian J of say, (2.2) and define the symmetric matrix

S = [1/2] ([J.sup.T] + J).

Theorem 3.3 in [11] says that if S has eigenvalues [[lambda].sub.1] [greater than or equal to] [[lambda].sub.2] [greater than or equal to] ... [greater than or equal to] [[lambda].sub.n] such that [[lambda].sub.1] + [[lambda].sub.2] < 0, then no periodic solution can exist for system (2.2). Numerical calculations indicate that for any

R [greater than or equal to] 2, P [greater than or equal to] 0, V [greater than or equal to] 3255 and H [member of] [0, 124],

it holds indeed that [[lambda].sub.1] + [[lambda].sub.2] < 0, and therefore there are no periodic solutions on those intervals.

Explicit expressions for the eigenvalues of S in terms of H, R, P and V are too complicated to be useful, especially if we consider the full system (2.2). However if in particular we assume that the entire population of vectors, dogs and donkeys are infected, then we can arrive to a conclusion about periodic solutions.

Theorem 3.2. The system (2.2) has no periodic solutions when entire population of vectors, dogs and donkeys are infected.

Proof. The matrix S takes the form

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

In such a case, the eigenvalues of S are (in order of magnitude)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where a = -[[beta].sub.H]V - [[gamma].sub.H], b = [1/2] [[beta].sub.H](A - H), g = -[[beta].sub.P]P - [[beta].sub.R]R - [mu].

One can readily verify that we always have [[lambda].sub.1] + [[lambda].sub.2] < 0, and therefore there are no periodic solutions to such system.

Invariance. It is possible to explicitly define some sets that are invariant with respect to the flow of the dynamical system. To this end, let <[tau] stand for one of the relations previously defined <, [less than or equal to], [much less than] and let

[P.sub.+] = {x [member of] [OMEGA]: f(x) [>.sub.[tau]] 0}, [P.sub.-] = {x[member of]e [OMEGA]: f(x) [<.sub.[tau]] 0}.

For (2.3) it is readily checked that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The following proposition is a direct consequence of system (2.3) being cooperative. See [13].

Proposition 3.3.

(a) [P.sub.+] and [P.sub.-] are positively invariant.

(b) If we denote [[PHI].sub.t]([x.sub.0]) the solution of (2.3) with [[PHI].sub.0]([x.sub.0]) = [x.sub.0] and [x.sub.0] [member of] [P.sub.+] (resp. [x.sub.0] [member of] [P.sub.-]) then [[PHI].sub.t]([x.sub.0]) is non-decreasing (resp. non-increasing).

Remark. From the definition of invariance, we can say e.g. that the first inequality in the set [P.sub.+] gives an explicit upper bound on the number of infected humans that will result when starting with values of H, R and V that satisfy the given inequalities.

We now move in the direction of the spectral properties of Jacobian matrix associated to (2.3) specifically.

Theorem 3.4. Consider the system (2.3). Then,

(a) the origin is a saddle.

(b) The nontrivial equilibrium is a saddle if [alpha] < [[beta].sup.2.sub.R]BC, and depending on the sign of certain discriminant, it is either a stable node or a stable focus if [alpha] > [[beta].sup.2.sub.R]BC,

where [alpha] = [[beta].sub.R]([[beta].sub.R]BV + [[beta].sub.R]CR + [[gamma].sub.R]R + [mu]V) + [mu][[gamma].sub.R].

Proof. The critical points of (2.3) are, the origin [x.sub.0] = (0, 0, 0) and a non trivial one

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where we require that

[[beta].sup.2.sub.R]BC/[mu] > [[gamma].sub.R], (3.2)

so that the equilibrium [x.sub.1] lies in the positive octant. The general Jacobian of (2.3) is

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3.3)

At the origin, this reduces to

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and its eigenvalues can be computed directly. They are

[[lambda].sub.1] = -[[gamma].sub.H], [[lambda].sub.2,3] = -([[lambda].sub.R] + [mu]) [+ or -] [square root of [([[gamma].sub.R] + [mu]).sup.2] + 4([[beta].sup.2.sub.R]BC - [mu][[gamma].sub.R])/R]

and it is clear that with the condition (3.2), the origin is a saddle point, and therefore unstable. (Also observe that if had [[beta].sup.2.sub.R]BC/[mu] = [[gamma].sub.R] then the origin would be nonhyperbolic and if [[[beta].sup.2.sub.R]BC/[mu]] < [[gamma].sub.R] then it would be stable).

To study the eigenvalues of the Jacobian (3.3) at the nontrivial equilibrium [x.sub.1], we see that one of the eigenvalues is [[lambda].sub.1] = -[[beta].sub.H]V - [[gamma].sub.H] < 0. The other two eigenvalues come from the 2 x 2 lower right block. If we denote with 8 the determinant of that block and with [tau] its trace, then the two eigenvalues are given by

[lambda] = t [+ or -] [square root of [[tau].sup.2] - 4[delta]]/2, (3.4)

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Observe that [delta] can be written as:

[delta] = ([[beta].sup.2.sub.R]BV + [[beta].sup.2.sub.R]CR + [[beta].sub.R][[gamma].sub.R]R + [mu][[beta].sub.R]V + [mu][[gamma].sub.R]) - [[beta].sup.2.sub.R]BC = [alpha] - [[beta].sup.2.sub.R]BC.

Thus, if [alpha]< [[beta].sup.2.sub.R]BC, then [delta] < 0 and therefore there are two eigenvalues of opposite sign and [x.sub.1] is saddle. Similarly, if [alpha] > [[beta].sup.2.sub.R]BC, then [delta] > 0, and since we know that [tau] < 0, we can conclude that [x.sub.1] is a stable node if [[tau].sup.2] - 4[delta] [greater than or equal to] 0 or a stable focus if [[tau].sup.2] - 4[delta]< 0.

Several observations are in order:

(a) Biologically, the inequality (3.2) says that the ratio of the rate of infection of the main reservoir hosts (donkeys) to the rate of mortality of the vectors (scaled by the population product BC) will have to be greater than the rate of recovery of the infected reservoir host, so that there is some infected population. This is in fact a threshold condition for infection [2].

(b) Under the threshold condition (3.2), the origin is unstable. This means that a state with no infected individuals can be easily perturbed into a state with some infected individuals.

(c) Although it is not easy to determine the sign of [[tau].sup.2] - 4[delta] in general, for the application in mind, the parameter [mu] is many times larger that the other parameters (see Section 4). This makes the term [[tau].sup.2] much larger that 4[delta], and therefore [[tau].sup.2] > 4[delta].

(d) The inequalities above depend much on the quantity [[beta].sup.2.sub.R]BC, that is, the rate of infection in vector-donkey encounters, the total population of donkeys and the total population of vectors, in other words, the very carriers of the disease.

4. Numerical Dynamics

In this section we perform a numerical study of the system (2.1), where there is only one species (dogs or donkeys) considered as reservoir, and of the extended system (2.2), where both dogs and donkeys at the same time are considered to be reservoirs. Among other things we study sability of equilibria and periodic orbits, location of branch points, and special solutions (such as connecting orbits). We also perform pseudo arc- length numerical continuation, following the theoretical background and numerical algorithms from [4] and [12]. This amounts to an extension of the numerical analysis of the local dynamics sudied in [2] to the global numerical analysis for systems in three and four variables.

[FIGURE 1 OMITTED]

[FIGURE 2 OMITTED]

To begin with, for the values [3] [[beta].sub.H] = 0.000982, [[beta].sub.R] = 0.0000527, [[beta].sub.P] = 0.000085, [[gamma].sub.H] = 0.024, [[gamma].sub.R] = 0.0007, [[gamma].sub.P] = 0.003, [mu] = 0.42, A = 124, B = 28, C = 4807, D = 43,

we solve (2.1) with I = R, and we also solve (2.2). The initial values are H(0) = 0, R(0) = 6, P(0) = 4, V(0) = 0. In Figures 1 and 2 we show the numerical results for the four species as functions of time. As expected, all species are nondecreasing functions of time t, which means that in the absence of prophylaxis and vaccines, the number of infected individuals will increase as time evolves. Of special interest is the evolution of the number of infected humans H(t).

For a closer look at the dynamic behavior of the systems, let us drop for a moment the P component in system (2.1), which leaves us the system

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.1)

We want to perform numerical continuation on the parameter [[beta].sub.H]. That is, starting from a solution (in this case the trivial one) of (4.1), we let [[beta].sub.H] vary and generate corresponding solutions of the system, and at the same time we track for possible special points and solutions. It is clear in this case that varying [[beta].sub.H] only affects H whereas R and V remain fixed. Thus, we plot only H against [[beta].sub.H]. We can make continuation for large intervals of [[beta].sub.H], but of course values larger than say, [[beta].sub.H] = 0.2 would be unrealistic for the given application. No special solutions (periodic orbits, etc.) nor any bifurcation point was detected. See Figure 3.

However, when performing continuation with respect to the parameter [[beta].sub.R], a bifurcation (branch point), was detected at [[beta].sub.R] = 4.673665E - 04 for the small values H = 2.591760E - 04,R = 9.549715E - 05,V = 5.108258E - 05. See Figure 4. For these values of the parameters, the inequality (3.2) becomes an equality. Connecting Orbits. Let us now go back to the general model (2.2). If we fix all parameters except [[beta].sub.H], then besides the trivial solution, there are two more equilibria. One of them has some negative components, so that it is not of practical interest. The third equilibrium is given by

Q = (-0.13817 [[beta].sub.H] + 124 [[beta].sup.2.sub.H]/(0.000766193 + [[beta].sub.H])(-0.00111427 + [[beta].sub.H]), 19.6623, 20.2186, 31.3237), (4.2)

[FIGURE 3 OMITTED]

[FIGURE 4 OMITTED]

The general Jacobian of (2.2) is

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.3)

Due to the block triangular form of the Jacobian, we observe that - [[beta].sub.H]V - [[gamma].sub.H] < 0 is one of the eigenvalues and is the only one that depends on [[beta].sub.H]. For [[beta].sub.H] = 0.000982 and the fixed values of the other parameters, the origin is unstable and Q is stable. In fact, the eigenvalues of the Jacobian at the origin are

[[lambda].sub.1] = -0.02400000, [[lambda].sub.2] = -0.42442592, [[lambda].sub.3] = -0.00141438, [[lambda].sub.4] = 0.00214030.

This also tells us the unstable manifold of the origin is one-dimensional. Similary, the eigenvalues of the Jacobian at Q are

[[mu].sub.1] = -0.05475986, [[mu].sub.2] = -0.42489150, [[mu].sub.3] = - 0.00183740, [[mu].sub.4] - 0.00403915,

which indicates that the stable manifold of P is four-dimensional.

We are interested in computing a connecting orbit from the origin to Q. This is of great importance because such a curve explicitly describes the path from non infected population to the stable state of infected individuals. In general, for a dynamical system

[??] = f(x, [lambda]), x(t) [member of] [R.sup.n], [lambda] [member of] [R.sup.p], (4.4)

where f: [R.sup.n] x [R.sup.p] [right arrow] [R.sup.n], and for [M.sub.[+ or - ]]([lambda]) representing equilibria or periodic orbits, a solution x(t, [lambda]), t [member of] R of (4.4) is called a connecting orbit from [M.sub.-]([lambda]) to [M.sub.+]([lambda]) if

dist (x(t, [lambda]), [M.sub.[+ or -]]([lambda])) [right arrow] 0, as t [right arrow] [+ or -][infinity].

According to [4] and [12], the tangent spaces of the corresponding stable and unstable manifolds should span the whole space and the intersection of such manifolds should be transversal. This guarantees the existence of an isolated connecting orbit between the two equilibria. From these conditions, a fundamental relation or balance between the dimensions of the manifolds and the number of parameters is established as

p = [m.sup.u.sub.+] - [m.sup.u.sub.-] - [m.sup.c.sub.-] + 1, (4.5)

where, for the problem considered here, [m.sup.u.sub.+], [m.sup.u.sub.-], [m.sup.c.sub.-] represent the dimensions of the unstable manifold of P and that of the unstable and center manifolds of the origin. So that for this particular case we have p = 0 - 1 - 0 + 1 = 0, so that the connection between the two equilibria can be established without free parameters.

Remark. In general, if (4.5) is not satisfied, then we should either add parameters to the system or add conditions for the parametrization of a manifold of connecting orbits. More precisely, if p < [m.sup.u.sub.+] - [m.sup.u.sub.-] - [m.sup.c.sub.-] + 1, we should add ([m.sup.u.sub.+] - [m.sup.u.sub.-] - [m.sup.c.sub.-] + 1) - p parameters in order to have a well-posed problem, while if p > [m.sup.u.sub.+] - [m.sup.u.sub.-] - [m.sup.c.sub.-] + 1, then we can add p - ([m.sup.u.sub.+] - [m.sup.u.sub.-] - [m.sup.c.sub.-] + 1) constraints to select a unique connecting orbit.

We do not limit ourselves to compute one connecting orbit between the origin and Q but instead we compute a branch of connecting orbits through smooth continuation, for [[beta].sub.H] in the interval [0.00001, 0.05]. See Figures 5 and 6, where we plot the time evolution of the infected humans. Connecting orbits are very important in the global study of dynamical systems, acting as separatrices of solutions in the phase portrait, giving a better understanding of the behavior of nearby solutions.

For system (2.2) we have also performed continuation on the parameter [[beta].sub.H] on the interval [0.000982, 0.1], starting from the equilibrium point (4.2). Numerical continuation for larger values of [[beta].sub.H] is still possible but it is not of practical interest. For the value [[beta].sub.H] = 0.1, the corresponding value of infected humans is already about H = 123. No special solutions or bifurcations are found. For [[beta].sub.R] we have performed continuation on the interval [0.0000527, 0.0037734]. A critical point at the right endpoint of this interval is found to be

U = (121.0175, 27.9947, 41.5222, 991.16723), (4.6)

which is stable. The other nontrivial equilibrium has some negative components. Just as we computed connecting orbits from the origin to P above, connecting orbits from the origin to U can also be computed.

[FIGURE 5 OMITTED]

[FIGURE 6 OMITTED]

[FIGURE 7 OMITTED]

5. Final Remarks

In this work we have studied models of American cutaneous leishmaniasis in the form of continuous dynamical systems, extending research in the same area in the current literature. We have explicitly computed paths that connect different states and we have performed continuation as a problem parameter is allowed to vary continuously. We have argued on the nonexistence of periodic solutions and have provided a study of stability of solutions. A future avenue of research would be to consider perturbations on the system and study their impact on the overall dynamics; this may introduce other interesting mathematical concepts into play such as bifurcations or strange attractors. A ratio-dependent interaction in the sense of [10] between vectors and humans (or between other species) would generate mathematical models of great interest, though more complicated.

References

[1] Awerbuch T., 1994, Evolution of mathematical models of epidemics. Ann. N. Y. Acad. Sci., 740:232-241.

[2] Chaves L.F. and Hernandez M.J., 2004, Mathematical modelling of American cutaneous leishmaniasis: incidental hosts and threshold conditions for infection persistence. Acta Tropica, 5:1-8.

[3] Chaves L., Hernandez M., and Ramos S., 2008, Simulacion de modelos matematicos como herramienta para el estudio de los reservorios de la leishmaniasis cutanea americana, Divulgaciones Matematicas, 16:125-154.

[4] Dieci L. and Rebaza J., 2004, Point-to-point and point-to-periodic connections. BIT Numerical Mathematics, 44:41-62.

[5] Gratz N.G., 1999, Emerging and resurging vector-borne diseases. Ann. Rev. Entomol, 44:51-75.

[6] Miguel de Guzman, 1980, Ecuaciones Diferenciales Ordinarias. Alhambra, Madrid.

[7] Jack Hale, 1989, Ordinary Differential Equations. Malabar.

[8] Hirsch M.W. and Smith H., 2004, Monotone Dynamical Systems, Electronic version.

[9] Killick-Kendrick R. and Ward R.D., 1982, Ecology of Leishmania, Parasitilogy, 82:143-152.

[10] Leard B., Lewis C., and Rebaza J., 2008, Dynamics of ratio-dependent predator-prey models with nonconstant harvesting, Disc. Cont. Dyn. Syst. S, 1:303-315.

[11] Li Y. and Muldowney J., 1993, On Bendixon's criterion, J. Diff. Equations, 106:27-39.

[12] Rebaza J., 2007, Smooth Schur factorizations in the continuation of separatrices, Linear Algebra and its Applications, 421:138-156.

[13] Hal Smith, 1995, Monotone dynamical systems: An introduction to the theory of competitive and cooperative systems, Math. Surv. Mon. AMS, 41.

Jesus Carreno

Departamento de Matematica, Universidad de Oriente, Nucleo de Sucre, Venezuela

E-mail: carreno@sucre.udo.edu.ve

Teodoro Lara

Departamento de Fisica y Matematicas, Universidad de los Andes, N. U. Rafael Rangel, Trujillo, Venezuela

E-mail: tlara@ula.ve

Jorge Rebaza

Department of Mathematics, Missouri State University

E-mail: jrebaza@missouristate.edu

Printer friendly Cite/link Email Feedback | |

Author: | Carreno, Jesus; Lara, Teodoro; Rebaza, Jorge |
---|---|

Publication: | Global Journal of Pure and Applied Mathematics |

Article Type: | Report |

Geographic Code: | 3VENE |

Date: | Dec 1, 2009 |

Words: | 4027 |

Previous Article: | Some of properties of intuitionistic fuzzy normal subrings. |

Next Article: | On Abel-Grassmann's groupoids defined by vector spaces. |

Topics: |