On the stability and Hopf bifurcation of the non-zero uniform endemic equilibrium of a time-delayed malaria model

This article considers a time-delayed mathematical model of immune response to Plasmodium falciparum (Pf) malaria. Infected red blood cells display a wide variety of surface antigens to which the body in turn responds by mounting specific immune responses as well as cross-reactive immune responses. The model studied here tracks these infected red blood cells as well as the two types of immune responses. It is assumed that the immune responses are timedelayed, and hence a system of nonlinear delay differential equations is considered. The goal of the paper is to provide a vigorous analysis of the stability and Hopf bifurcation of the non-zero uniform endemic equilibrium of the mathematical model.


Introduction
The present paper is a development and refinement of some of the ideas first presented in [15] and [4].For the sake of completeness, we reproduce here some portions of the introduction of [15] -which provides the background of the mathematical model under consideration.We consider a time-delayed modification of the seminal Recker et al. [18] intra-host mathematical model of immune response to Plasmodium falciparum (Pf), a species of parasites that cause malaria in humans.The model incorporates the effects of time-delayed immune response (IR) mounted by the human host and was first introduced in the work of Mitchell et al. [12][13][14] as a natural development of the model proposed by Recker et al. [18].Over the years, there has been considerable work conducted on this type of model in the research literature (see [1][2][3][4][5][6][10][11][12][13][14]17] and [18] for example).In the original work of [18], the authors proposed a mathematical model of immune response to Plasmodium falciparum predicated on the hypothesis that a given antigenic variant experiences two different types of immune response.These are the long-lasting variantspecific immune response and the short-lived cross-reactive immune response mounted c Vilnius University, 2016 I. Ncube against a set of epitopes shared with other antigenic variants.The main achievement of the model is its ability to replicate the sequential appearance of dominant antigenic variants, each of which is the most immunologically distinct from its preceding types [18] -this being a strategy employed by the parasite population to evade the host's immune system [14].
The time-delayed modification of the Recker et al. [18] is expressible in the form [13,14] where Y i denotes the amount of antigenic variant i, Z i and W i denote variant-specific and cross-reactive immune responses, respectively, φ is the intrinsic parasite growth rate, α and α are the removal rates associated with specific and cross-reactive immune responses, respectively, β and β are the proliferation rates of immune responses, µ and µ are the decay rates of variant-specific and cross-reactive immune responses, and T d is the discrete time-delay of the IR.The coefficients c ij of the connectivity matrix characterise cross-reactive inter-variant interactions [2,6,18].Following in the footsteps of [13] and [14], we assume that all of the variants have identical temporal dynamics, namely, Y i (T ) = Y (T ), Z i (T ) = Z(T ), and W i (T ) = W (T ) for all i.
Let us define dimensionless variables y, z, and w as deviations from the endemic steady state and a new dimensionless time [14]: We now recast (1) in the new dimensionless variables to obtain the system where y| τ := y(t − τ ), τ is the rescaled time-delay, and http://www.mii.lt/NAHere, the term nβ denotes the growth rate of the cross-reactive immune response, and n denotes the number of variants that share minor epitopes.In an effort to facilitate the analysis, [14] introduced the new variable x = z + w representing the total IR.Equation (3) thus becomes [14] x System (5) has two different types of equilibria, namely: 1) a disease-free equilibrium E 0 := (x 0 , y 0 , z 0 ) = (−(1 + qb)/(ab), −1, −q/a) and 2) a non-zero uniform (endemic) equilibrium E 1 := (x 1 , y 1 , z 1 ) = (0, 0, 0).We must comment that system (5) describes the dynamics of deviations from the uniform endemic equilibrium.The details of this derivation can be found in [15,18] and [4].Thus, even though the equilibrium E 1 of system ( 5) has all its components identical to zero, it infact corresponds to the non-zero uniform endemic equilibrium of the original system (1) [4].Likewise, the equilibrium E 0 is the genuine disease-free equilibrium of (1) [4].It is worth noting at the onset that the studies presented in [12,13] and [14] primarily focus on the dynamics of the non-zero uniform endemic equilibrium.Employing asymptotics and perturbation analysis techniques, the authors show that a range of interesting dynamics result as a consequence of the time delay.The current paper seeks to investigate the stability and Hopf bifurcation of the non-zero uniform endemic equilibrium from the viewpoint of the vigorous direct analysis of the associated characteristic equation.
In this article, we study a mathematical model of Plasmodium falciparum malaria [12][13][14]17,18], which is given by the system of DDEs given in (5).As established in [15], when τ = 0, the linearisation of ( 5) about the non-zero uniform endemic equilibrium E 1 = (0, 0, 0) is given by Consider solutions of (5) of the form where c 1 , c 2 , c 3 ∈ R and λ ∈ C. Substituting ( 7) into (6) yields the system Nonlinear Anal.Model.Control, 21(6):851-860 which can be re-arranged and simplified to give the matrix equation whose solutions c 1 , c 2 , and c 3 are non-trivial if, and only if, Denoting the above determinant by D(λ), we see that c 1 , c 2 , and c 3 are non-trivial if, and only if, which defines the characteristic equation of ( 6) about the non-zero uniform endemic equilibrium point E 1 = (x 1 , y 1 , z 1 ) = (0, 0, 0), and where λ ∈ C is the characteristic exponent.

Equilibria and Hopf bifurcation
Linear stability of the non-zero uniform (endemic) equilibrium in the case of instantaneous IR was studied in [12,13] and [14].In [12,13] and [14], the authors also studied the linear stability of the non-zero uniform (endemic) equilibrium in the case of delayed IR, where the IR delay was taken to be discrete, very small, and fixed for all episodes of infection.In the cited studies, the authors relied heavily on perturbation and asymptotic analysis techniques.Their calculations were somewhat simplified by their assumption that the IR delay τ is very small.In Blyuss et al. [3], the authors did consider the case of an arbitrary immune response time delay and subsequently studied the phenomenon of symmetry breaking.In this article, we systematically investigate the local stability and Hopf bifurcation of the equilibrium E 1 = (0, 0, 0).Our approach differs from previous studies primarily because we rely on the direct analysis of the associated characteristic equation (8).We begin our study by noting that the characteristic equation ( 8) can be expressed in the form Setting τ = 0 in (9) leads to the equation Using the well-known Routh-Hurwitz criterion, we conclude that polynomial (10) is stable if, and only if, http://www.mii.lt/NABecause of the fact that all of the parameters in the Recker et al. model [18] are strictly positive, it follows that conditions (11) are always valid.Let λ = iω (ω > 0) in ( 9), and decompose into real and imaginary components to get We obtain from (12) that Using the identity sin 2 (ωτ ) + cos 2 (ωτ ) = 1, equation ( 13) yields the polynomial where When the Recker et al. model [18] is derived, all parameters are taken to be positive with a very specific biological meaning.The implication of this fact is that the parameters a, b, and q defined in equation ( 4) are positive.As a consequence of this, q 5 defined in equation ( 15) is always negative.Let us introduce the change of variable z = ω 2 in ( 14) to obtain H(z) := z 5 + q 1 z 4 + q 2 z 3 + q 3 z 2 + q 4 z + q 5 = 0.
Using MAPLE, we find that the transversality condition is satisfied if, and only if, where In order to simplify (17), one would typically substitute (13) into the expression.However, in this case, the resulting expression is too lengthy to include here.In the light of the foregoing discussion and the general Hopf bifurcation theorem of functional differential equations [8], we arrive at the following result concerning the stability of the non-zero uniform endemic equilibrium of system (5) when τ > 0.

Stability region in an appropriate two-parameter space
In this section, we discuss the asymptotic stability boundary of the equilibrium E 1 in the (b, q)-parameter space.We note at the onset that the parameters b and q appear directly in system (5).To achieve this goal, we attempt to go via an intermediate step in which we may conduct our analysis in an intermediary two-parameter space [7].This approach is along the lines of the work of [7], who did a similar analysis for a nonlinear scalar delay differential equation arising in the work of [16].Unfortunately, the generalisation of their technique to systems of delay differential equations remains an open problemthis article is a modest attempt to address this problem in the particular case of a system of three nonlinear delay differential equations.
To begin, we now prepare system (5) for the analysis to come by rescaling some of its parameters.In particular, let t := t/τ , q := qτ , a := aτ , and b := bτ .As a consequence of this, we obtain the following rescaled system: For notational convenience, we now drop the tilde's from system (18), leading to its less cluttered counterpart The linearisation of (19) about the equilibrium E 1 = (0, 0, 0) is given by and the corresponding characteristic equation is Instead of directly analysing the characteristic equation ( 21), we shall now digress and investigate whether we can, instead, work with an intermediary characteristic equation, as suggested and discussed in [7].To begin, let us express (20) using matrix notation to obtain where Equation ( 22) admits non-trivial solutions of the form χ(t) = ce λt , λ ∈ C, and c = (c 1 , c 2 c 3 ) T = 0, with c 1 , c 2 , c 3 ∈ R, if, and only if, λ satisfies the characteristic equation where I is the 3 × 3 identity matrix.Equation ( 10) in [7] is a scalar analogue of (22).
Comparison of these two equations makes it abundantly clear that the two-parameter technique championed in [7] may not be possible in dealing with characteristic equations associated with systems of dimension greater than one.We must remark that equation ( 23) is equivalent to (21).It is evident that (21) is complicated, and any attempt aimed at studying the distribution of its zeroes directly is bound to be laborious.The idea behind the derivation of an equivalent and 'simpler' intermediate characteristic equation such as (23) (or equation ( 2) in [7]) is so that we may be able to obtain some insight about (21) by studying (23) in some appropriate 'intermediate' two-parameter space.In the one dimensional case, the intermediate parameters used in the two-parameter analysis would typically be selected to be the coefficients of χ(t) and χ(t − 1) in ( 22).This idea is well-articulated for the one-dimensional case studied in [7], and where the intermediate characteristic equation takes a very simple and prototypical form.The current effort investigates whether idea of a two-parameter analysis promulgated in [7] is extensible to systems of dimension higher than one.It is now clear that analysis of (23) in an appropriate intermediary two-parameter space as posited in [7] is not straightforward in the three-dimensional case, primarily because our two intermediate 'parameters' A and B are matrices (not scalars as in equations ( 2) and (10) derived in [7]), and computing the determinant in (23) will not preserve the matrix structure of these parameters -they will simply disappear during the ensuing calculation.At this point, we turn to looking at (21) directly.Let λ = iω, ω > 0, in (21) and decompose the resulting expression into its real and imaginary components, thus: We then solve (24) for the parameters (b, q) as functions of the frequency ω.We mention that two sets of values of the 2-tuple (b, q) are obtained, as clearly b and q occur nonlinearly in (24).Let us denote the pair of solutions so obtained as (b 1 (ω), q 1 (ω)) and (b 2 (ω), q 2 (ω)).Furthermore, let us denote http://www.mii.lt/NAwhere ξ 1 = 4ω 4 cos 2 ω + 4aω 3 sin(2ω) + 4a 2 ω 2 sin 2 ω, Then, we have the following solutions (b(ω), q(ω)) of ( 24 When b = 0 = q, the characteristic equation (21) simplifies to (λ + a) λ 2 + τ 2 e −λ = 0, whose only root is λ = −a < 0. This implies that any region containing (0, 0) in the (b, q)-parameter space is necessarily the asymptotic stability region.In such a region, the quasi-polynomial (21) has no roots whose real part is located in the right half-plane.As one crosses the boundary of the stability region, two roots whose real part is in the right half-plane appear via a Hopf bifurcation of the equilibrium E 1 = (0, 0, 0).

Conclusion
This article has focussed on the non-zero uniform endemic equilibrium E 1 of the celebrated Recker et al. model [18] endowed with an arbitrary immune response time delay.In the process, we were able to establish concrete conditions under which a Hopf bifurcation occurs at the equilibrium E 1 .In Section 3 of the article, we attempted to give a description of the asymptotic stability boundary of E 1 in an appropriate two-parameter space.In carrying out this analysis, we also showed that the technique advertised vigorously in the work of [7] may not be extensible to higher dimensional systems, unfortunately.