Critical behavior of the delay-induced chaos transition in a nonlinear model for the immune response

Critical behavior of the delay-induced chaos transition in a nonlinear model for the immune response

Elder de Souza; Marcelo Lyra; Iram Gleria

Instituto de Física - Universidade Federal de Alagoas 57072970 - Maceió, AL, Brazil

ABSTRACT

In this paper we analyze a model for the dynamics of the immune system interacting with a target population. The model consists in a set of two-dimensional delayed differential equations. The model is effectively infinite dimensional due to the presence of the delay and chaotic regimes can be supported. We show that a delayed response induces sustained oscillations and larger delay times implies in a series of bifurcations leading to chaos. The characteristic exponent of the critical power law relaxation towards the stationary state is obtained as well as the critical exponent governing the vanishing of the order parameter in the vicinity of the chaotic transition.

Keywords: Delayed systems, Immune response, Chaos

1. INTRODUCTION

Mathematical models for the dynamics of the immune system received considerable interest in physics literature [1]. The models are simplified versions that try to mimic the complex behavior of real systems. In spite of this simplicity, such models can provide insights on the real situation and help to understand the intricate dynamics involved.

The more complicated models may display more accurately the complex patterns observed in real data, but there are difficulties concerning numerical estimations of the parameters and the meaning of the many interaction terms that appear in such high dimensional models [2]. Alternatively one may propose low dimensional models that present some of the observed basic phenomena. The model of H. Mayer et al proposed in [3] follows this approach. It models the interaction of the immune system with a target population (e. g virus or bacteria). A two dimensional set of differential equations with nonlinear interactions between specific immune components (the immune competence E measured by the concentration of certain immune cells as, for example, cytotoxic T lymphocytes cells (CTLs), killer T cells, etc) and the target population T (e. g bacteria and viruses) is considered. The model is as follows:

In the above system, r > 0 represents the reproduction rate of the target population and k > 0 the elimination rate of the targets by the immune cells. The first term in the dynamics of E represents the speed of the processes triggered by the targets T leading to their elimination by E. For example, non-specific precursor cells or inactivated T-cells are transformed into specific helper cells, CTL cells or plasma cells in the presence of targets, producing certain antibodies. All the constants (p, u, ν) are positive and u < ν. The second term represents the autocatalytic effect of immune responses. Finally, the term - E reflects the finite lifetime of E.

In the above model, the complex interactions of the immune system with the foreign agents T are represented by few measurable quantities. The nonlinearities in the dynamics of E and T are the only source of complexity in the model. Since it is a two-dimensional set of ODE's, only regular trajectories are observed. In many cases real time series data of the immune state look rather irregular, suggesting chaotic behavior. See for example the data shown in [3] for the number of killer cells versus tumor size during a metastatic process of Fibrosarcoma.

In [5] and [6], the Mayer et. al. model was modified to include time delays and the possibility of delay induced chaos was explored. Here, we also consider variants of (1) with delayed responses. We focus on a simple model given by:

If b 0, we are admitting a competition for the limited resources available to the target population. The sigmoid behavior of means that the immune system ignores the T population if it is below a given threshold . Furthermore, this term is bounded reflecting the fact that the precursor E population is limited. Our motivation is to search in this simple model some patterns observed in real situations. In the next section, we show that, as expected, the delayed response may induce oscillations in T and E.

It is well known that the introduction of the delay makes the model effectively infinite dimensional. Chaos can thus be observed in one dimensional delayed systems [4].

In the next section we show that the delayed response induces sustained oscillations and larger delay times imply in a series of bifurcations leading to chaos. The characteristic exponent of the critical power law relaxation towards the stationary state is obtained as well as the critical exponent governing the vanishing of the order parameter in the vicinity of the chaos transition.

We compute the critical delay-time above which the stationary solution is unstable. The dynamics of the system in the neighborhood of the critical point is characterized by the dynamic critical exponent. We show that the bifurcations evolve to a chaotic regime.

The rest of the paper is organized as follows. Next section presents some numerical results concerning the behavior of (2) and section 3 concludes.

2. SOME RESULTS

Let us consider (r, k, p, s, b) = (3,2,2,1.25,0.3) [3]. The (real) fixed points of system (2) are:

The Jacobian matrix is given by:

It is straightforward to show that the trivial stationary solution is unstable for all values of τ. For null delay, the interior fixed point is stable. When τ 0 the analysis is performed with the help of the following theorem [7]:

Theorem 1: Consider the characteristic equation P(λ)+ Q(λ)e^-λτ = 0 associated to a fixed point, where P(λ) and Q(λ) are analytic functions in a right half-plane Reλ > -δ,δ > 0 which satisfy:

i) P(λ), Q(λ) have no common imaginary zero, ii) where the bar means conjugate; iii) P(0)+ Q(0) 0; iv)There are at most a finite number of roots of P(λ)+ Q(λ)e^-λτ = 0 in the right half plane if τ = 0; v) F(y)= |P(Iy)|² -|Q(Iy)|²with y

has at most a finite number of real zeros.

Then the following statements are true:

a) Suppose F(y) = 0 has no positive real roots. Then if the associated fixed point is stable (unstable) for τ = 0 it remains stable (unstable) for all τ > 0.

b) Suppose F(y)= 0 has at least one positive root and that each positive root is simple. As τ increases, stability switches may occur. There exists τ* > 0 such that P(λ)+ Q(λ)e^-λτ = 0 is unstable for all τ > τ* . As τ varies from 0 to τ* at most a finite number of stability switches may occur.

The characteristic equation for the interior fixed point is given by:

with B = 2.534071424, C = 1.190945352, Q = 0.961886711. Constructing F(y) as above, we find the following roots for F(y)= 0: y =[±0.8837551400,±1.086095210I], which means we are in case b) of the above theorem. In fact, the critical value of the delay, say τ*, above which the point loses its stability, can be obtained as follows. Let λ = u + Iω, the critical delay then occurs when u = 0. With this substitution in (5) we obtain:

which can be solved numerically giving τ* 0.823.

From the discussion above, we conclude that if τ < τ* the system presents damped oscillations toward the stable equilibrium. A delay τ > τ* makes this fixed point unstable and, after an initial outbreak in the T population, we expect to observe sustained oscillations in the concentrations of T and E. Solving (2) with a modified Runge-Kutta method, we show in figure 1 the behavior of T for τ > τ* . Note that the minimum T population is smaller than the equilibrium value of the stationary solution. As pointed out in [4], this may bring advantages to the immune system because other strategies for the elimination of T may be more efficient when the T concentration is lower than the equilibrium value.

We characterized the critical behavior at τ = τ*, by computing the amplitude of the E oscillations as a function of the delay-time. In Figure 2, the amplitudes of oscillations (after the transient) for the E population versus τ are shown. At the critical delay τ = τ*, the amplitude decays with a power-law. In figure 3 we consider the amplitudes in the T population, say A(T ). We have A(T) ∝ (τ - τ*)^β as τ → τ*₊, with β = 0.5 as seen in figure 3. The damped oscillations also become critical. Let A(t; T ) denote the time evolution of the damped oscillations. We observed a power-law behavior A(t; T ) ∝ t^-φ, see figure 4. We estimate the critical exponent as being φ 0.5 (within the error bar).

In figure 5, we show the local maximum values of E as a function of the delay-time. Increasing the delay time, we observe bifurcations in the delay induced oscillations. The bifurcations reflect the emergence of new oscillatory patterns. Large delays induce a series of bifurcations to chaotic behavior. The route to chaos can be seen in figure 5 where E_max is plotted against τ. Windows of regular trajectories are present within the chaotic regions.

Figure 6 shows some typical phase portraits, where the emergence of chaos can be observed. Figures 6 a) and b) show typical oscillatory behavior for τ > τ*. We illustrate the series of bifurcations that appears as we increase the delay. Figure 6 c) and d) shows typical chaotic trajectories.

3. CONCLUDING REMARKS

We analyzed a two-dimensional system with delay that models the dynamics of the interaction of the immune competence E with a target population T . Several different behaviors are present in this model due to the presence of nonlinear functions in the dynamics of E. The dynamics of the delayed model is, in some aspects, richer than that of the original model. For example, it allows the presence of chaos.

We showed that the delay induces bifurcations and the stationary solution becomes oscillatory above a critical delay. The characteristic exponents of this bifurcation and the critical dynamics were obtained. Larger delay times induces bifurcations which ultimately lead to chaos, which is forbidden in the original non-delayed model.

Low dimensional models have the advantage that several phenomena may be viewed as the result of few basic mechanisms [3]. High dimensional models with many variables may furnish better computational results, but our work reinforces the fact that simple models, based on a few number of easily interpreted parameters, are able to simulate some phenomena observed in real situations.

4. ACKNOWLEDGMENTS

The authors acknowledge financial support from CNPq, CAPES and FAPEAL.

(Received on 6 January, 2009)

• [1] Wang K, Wang W, Pang H, Liu X. Physica D, 2007; 226: 197.
• [2] Perelson AS and Weisbuch G. Review of Modern Physics, 1997; 69-4: 1219.
• [3] Mayer H, Zaenker KS and An der Heiden U. Chaos, 1995; 5-1: 155.
• [4] Canabarro A, Iram Gleria and Lyra ML. Physica A, 2004; 342: 234.
• [5] Buric´ N , Mudrinic M and Vasovic´ N. Chaos, Solitons and Fractals, 2001; 12: 483.
• [6] Elder de Souza, Marcelo Lyra and Iram Gleria. pre-print.
• [7] Cooke KL and Van den Driessche P. Funkcial. Ekvac., 1986; 29: 77.

Publication Dates

• Publication in this collection
  10 Sept 2009
• Date of issue
  Aug 2009

History

• Accepted
  06 Jan 2009
• Received
  0000