Correlated dynamics of immune network and sl(3, R) symmetry algebra

1 Introduction

The understanding of coupled multi-component dynamical system (steady state or bifurcation) requires mathematical understanding of the system manifold. In particular, generic bifurcation theory with symmetry, normal forms and unfolding theory all make vital contributions to explaining and predicting behavior in such systems. We consider pathogen dynamics model or immune network model where immune response in target invasion and proliferation in body system is considered to follow very nonlinear/complex path. The growth and interaction pattern follows a nonlinear predator-prey-type interaction. CD8+ T cells are one of the most crucial component of the adaptive immune system that play key role in response to pathogen. Upon antigen stimulation, naive CD8+ T cells get activated and differentiate into effector cell. This mechanism may create small subset of memory cell after antigen clearance. Emerging evidences support that metabolic reprogramming not only provides energy and biomolecules to support pathogen clearance but also is tightly linked to T-cell differentiation [3]. In brief, naive CD8+ T cells are activated in response to the coordination of three signals, including T cell receptor, co stimulation and inflammatory cytokines via multistep strongly connected complex pathway. The majority of CD8+ T cells undergo contraction phase and die by apoptosis [3]. Therefore, we try to undertake nonlinear pathogen dynamics mathematical model for present studies. Since symmetry is a fundamental invariant structure associated with various mechanical/physical systems, it influences the functionality of the dynamical system. In systems, with Hamiltonian dynamics, the equation of motion exhibits symmetry in which total energy conserved (Noether Symmetry). On the other hand, various physical dynamical systems exhibit symmetry feature that conserves action of dynamics and gives rise to Euler-Lagrange equation as equation of motion. Many dynamical systems represented by coupled autonomous equations exhibit presence of attractor (local or global) to sustain stability of the dynamics. The work of Ashwin et al. [2] showed connection between transitive symmetry algebra and periodic behavior of dynamical system, which indicates reflection of pattern behavior through existence of symmetry structure. Symmetric attractor is a signature property of equivariant dynamical system. Continuous group such as compact Lie group is used in many mathematical models to understand connection between symmetry and invariant quantity in the dynamics. Followed by such idea, we try to explore pattern behavior under action of continuous group symmetry. Since Lie group action under one parameter point transformation leaves the manifold invariant (under linear vector field), this can be used to unfold evolution structure to obtain pattern structure at any time. Moreover, maximal Lie algebra for a second-order ordinary differential equation (ODE) can leave a diffeomorphic manifold invariant; this can be implemented to integrate nonlinear ODE. In most autonomous equations, symmetry algebra is transitive in nature (group generators follow simple time translation and population growth). Conn [6] described transitive Lie algebra over a ground field K (real or complex field) as topological Lie algebra whose underlying vector space is linearly compact, which possesses a fundamental sub-algebra with no ideal (opposite to the case of primitive action algebra).

From the standpoint of the geometric analysis of Lie algebra, the generator should take the form

which we view as formal vector field under Lie infinitesimal transformation. Under Lie group of infinitesimal transformation, system follows connected manifold. Given a connected differential manifold M and the action of a compact Lie group g on M , g t represents the isotropy subgroup of g at any time t for a dynamical system represented by autonomous ODE. Normalizer N ( g t ) of the subgroup must exist corresponding to t ˜ . If g ( t ) and g ( t ˜ ) are equivariantly diffeomorphic then g t and g t ˜ are conjugate subgroup of g . Equivalently, g ( t ) and g ( t ˜ ) are equivalently diffeomorphic. If v → is an equivariant vector field (group generator), then v → f = 0 , where f represents differential equation of the dynamical system. Local pseudogroup (local diffeomorphism) under transitive Lie algebra preserves structure of the manifold very much, which is very important in stability of dynamics. Theorem of Guillemin and Sternberg [11] asserts that, given a transitive Lie algebra L and a fundamental sub-algebra L 0 ⊂ L , one can realize L as a transitive sub-algebra of formal vector fields in such a way that L 0 is realized as isotropy sub-algebra of L ; such realization of ( L , L 0 ) is very unique under formal change of coordinates. This means under group action (vector field v → ), flow t → exp ( t v → ) is maintained. In most physical problems dictated by Hamiltonian of the system (Kepler’s law of planetary motion), infinitesimal transformation of the variables under Lie group of continuous transformation showed momentum conservation (Noether symmetry). In system dictated by Lagrangian (action integral) of the system, such group symmetry can manifest in obtaining some of the first integral (under action of proper sub-algebra), which can be related to Lagrangian of the system.

In this work, we construct evolution structure driven by the presence of symmetry (Lagrangian or Hamiltonian). Since biological evolution does not follow conservative system, one can assume the system is driven by Lagrangian action (followed by Euler-Lagrange equation). The stable structure of the dynamics is intrinsically connected to its inherent symmetry to the system. Once we are able to obtain such symmetry, that is used to obtain analytical structure of the solution, which means that irrespective of the initial condition, the solution will follow the same behavior globally/partially. The work by Ashwin et al. [2] gave necessary condition for a subgroup of a finite group (solvable sub-algebra) to have symmetry of a chaotic/nonchaotic attractor. Their numerical study showed that if any solution of the dynamical system (described by ODE) is obtained from G invariant system of ODE partial differential equation, where G is a compact Lie group of symmetries, then even if solution varies chaotically in time, there exists some invariant quantity under nontrivial subgroup of symmetries in G . The observation of periodic behavior in solution can describe system having symmetry. With such extensive work and results by several groups to understand the existence of symmetric structure of a dynamical system and corresponding regular pattern of solution, we try to investigate immune dynamics model to study relation between existence of symmetry and solution structure. Our goal is to study possible existence of symmetry algebra of continuous group under invariance and use that to develop analytical structure of solution if it exists. In this context, we discuss methods to construct vector fields that act continuously on variables on diffeomorphic manifold (under action of compact Lie group and transitive algebra). Since Lie group with maximal symmetry can have canonical coordinates that converts an ODE into quadrature form, it will be easy to obtain solution in linear sub space. Once solution structure is obtained, one can add nonlinear term into linear part to study the effect of nonlinear term. We adopt group theoretic approach to construct symmetry structure and differential invariant of the manifold under action of G . Our goal is to find existence of Lie symmetry obeyed by ODEs conditionally or globally if it exists and then use it systematically to obtain analytical structure of the solution. Since Lie point symmetry vector field is linear, we can assume under action of the symmetry group, system preserves dynamics of least path of action or action of Lagrangian is preserved. Under Lie point group action g , if the dynamical system exhibits convergence properties, then system can contract toward local fixed point. The characteristics of many physical dynamics are that system possesses intrinsic symmetry characteristics of a certain kind of conservation law under symmetry algebra. The first integral method is very significant as the first integral is associated with the Lagrangian of the dynamics. The first integral may confine the solution to a bounded region of phase space. It is very important to have knowledge about the analytical structure of a solution to understand the dynamical system. Certain mathematical community have devoted their research to the algebraic structure of various point symmetries in various dynamic system.

In order to obtain analytic solution of the ODE through symmetry, method involves reduction of order through construction of canonical (normal) subspace using solvable sub-algebra. Derived algebra of a Lie algebra G is analogous to commutator subgroup of a group. It consists of all linear combinations of commutators and is clearly an ideal. For an r -dimensional Lie algebra G r , relation is given by

in terms of structure constant C γ a b . g r is r -dimensional solvable algebra if there exists a chain of sub-algebras

such that G ( k ) is a k -dimensional Lie algebra and g ( k − 1 ) is called an ideal of g ( k ) . To obtain solution in differential invariant subspace, it must satisfy

for any positive λ . First, method involves construction of normal form of generators in linear subspace of canonical variables and invertible mapping. Under solvable sub-algebra, normal form of generators convert ODE into quadrature form. The second method involves the construction of a differential invariant function in the space of the invariant function, which can then be utilized to construct a solution. This method requires k th extended generator formation of an ODE.

Any equivariant dynamical system possessing Lagrangian/Hamiltonian structure or any kind, should possess recurrent robust attractor [2]. Since symmetry plays fundamental role in many physical/mathematical problem, our main focus will be to develop condition for symmetry that drives immune network. Many systems in nature possess intrinsic dynamical symmetry. Most biological mechanics in nature (non conservative system) posses intrinsic dynamical symmetry in which evolution dynamics is dictated by Lagrangian action functional. In such systems, Lagrangian function remain invariant associated with symmetry. It can be assumed to have rich interplay between symmetry property and dynamical behavior. The experimental work of Ma et al. [14] indicated periodic behavior of infection phase in rubella infection. Since most clinical data takes average value from blood sample, it cannot reflect the detal dynamics in infection and chronic phase of the disease. The clinical data by Liao et al. [13] and references therein in 18 such pathogen-borne infection cases suggest periodic nature of the infection and related symptoms. This means fever or other external symptom follows up-down path over time till it disappears finally or requires external intervention to annihilate target proliferation. All these clinical results support complexity in interaction path. Keeping complex nature of immune-target network, we consider nonlinear autonomous equation for pathogen dynamics in next section.

Following is our plan of work:

Section 1 is the introduction. Section 2 describes basic immune dynamics model with various features, and Section 3 describes the dynamical analysis of the model in detail. In Sections 4 and 5, we construct method of fundamental symmetry generators under infinitesimal transformation (under transitive algebra). Results of numerical simulation is also shown using group theoretic structure of the solution in the last section. In this work, we try to understand the pattern of growth behavior and corresponding interaction phase space.

2 Immune dynamics model

In case of any target invasion or infection (bacteria/virus/immunologic tumor cell ) in body, two types of immune cells, namely effector and memory cells play key role as immune response in combating such infection in short or long term. The proliferation and interaction of target cell in body are multicomponent/multistep nonlinear phenomena following the idea of predator-prey dynamics. The dynamics can be represented as follows:

where the first term designates self-proliferation and the second term denotes mutual interaction. In case of major two-component pathogen dynamics, system is described [15] as follows:

where x ( t ) and y ( t ) represent target cell and immune cell density in local volume at any time t . The term f ( x ) represents velocity of immune stimulation by target invasion x ( t = 0 ) , which leads to competence against them in the network with u and v as degree of stimulation, respectively. The form of f ( x ) and g ( y ) represent Richardson-type logistic functional growth in model of competence. The term g ( y ) represents autocatalytic enforcement of the network, which is very necessary to acquire adaptive immune memory cell infection, as a manifestation of chiral autocatalytic origin.

The parameter m represents threshold presence of target cell that can be recognized as signal by immune network. Parameter d represents constant death rate of immune cell. Immune competence y ( t ) can be defined as elimination capacity of the immune system with respect to target [15] and can be measured by the concentration of certain cytotoxic T-cell, natural killer cell or by concentration of certain antibodies. Cross talks between antigen presenting and T-cell impacts cell homeostasis amid bacterial infection and tumorigenesis [10] dynamics. The condition u = v is characterized by no target burden factor or equivalently target and immune cell growth are in competence [15]. We chose this nonlinear model of infection disease in our studies of the role of symmetry algebra and how to obtain analytical solution structure. This type of growth pattern of immune cell is noted in tumor dynamics and various infectious diseases [13].

3 Dynamical analysis of the model

Through dynamic analysis, we try to study robustness of this kind model in terms of stable invariant phase space.

Here, we implement no effective immune cell present at t = 0 as a trigger effect of f ( x ) and study immune network growth pattern. Following stability analysis, we obtain one equilibrium point of zero target and positive nonzero immune cell with following relationship of essential parameters, which work as basin of attraction of the two-component dynamics, namely 0 , d ( c d − s ) s 2 with

where m in model represents threshold value of target population that triggers immune network. Corresponding eigenvalues of Jacobian are evaluated as follows:

Dynamic simulation data (Figures 1 and 2) exhibit correlated dynamics between two components in the network through production of so-called proper antigen mechanism path.

Figure 1

Linear growth pattern; d = 1.0 , r = 0.701 , p = 0.642 , s = 1.23 , c = 1.0 , u = v = 0 , k = 0.9255 , m = ‒ 0.14 .

Figure 2

Periodic behavior of growth dynamics; d = 1.0 , r = 0.701 , p = 0.642 , s = 1.56 , c = 1.0 , u = v = 1 , k = 0.12 , m = ‒ 0.12 .

Linear dynamics endowed by u = v = 0 shows the pattern behavior in both components in a closed-phase trajectory (Figure 1) around the point of attraction.

Adding nonlinearity into simulation setting u = v = 1 , we increase s (velocity rate of immune cell interaction) to 1.53; the system exhibits oscillation phase even target interaction rate k is low. This is shown in Figures 1 and 2 where very stable phase trajectory is exhibited. The coexistence of stable pattern of both components is illustrated in Figure 2. The dynamics follows limit cycle characteristics with negative value of m showing equilibrium point.

Our analysis reveals the significant role of m in determining the stability of the dynamics ( − 0.1 < m < 0.3 ). When value of m shifts from this range in either direction, coexistence of stable pattern disappears. The oscillation becomes more prominent and stable for a longer time when we set v = 3 for negative value of m , as shown in Figure 3.

Figure 3

Periodic orbit presence in correlated dynamics; d = 1.0 , r = 0.701 , p = 0.642 , s = 1.56 , c = 1.0 , u = 1 , v = 3 , k = 0.12 , m = ‒ 0.9953 .

This regime can be recognized as interacting phase of various T/B cells via cytokinase production mechanism in order to achieve immunity for longer period marked by m = − 0.9953 . Numerical data exhibit characteristic pattern change after 9–12 h of cell infection. Figures 2 and 3 exhibit characteristic pattern with stable closed trajectory around local attractor. This phase can also be termed as long-term antibody production phase and is significant in rubella, German measles, influenza, smallpox or any other lethal diseases. Oscillation occurs in the presence of threshold target cell ( x c + m ) recognized by immune network. This is the case when the immune system recognizes the small presence of a target cell through the presence of a special kind of T cell in the system. Existence of multi layer component in many lethal diseases such as rubella, HIV are very common. It is evident that ratio of target proliferation p x u x v is very dominant to sustain pattern coexistence. Upon increasing a u and v in the network, we observe the existence of chaotic phase trajectory shown in Figure 4.

Figure 4

Chaotic phase trajectory; m = 0.9995 , d = 1.0 , r = 0.701 , p = 0.642 , s = 1.72 , c = 1.0 , u = 2 , v = 4 , k = 0.265 .

This kind of oscillatory/chaotic behavior can be termed as indeterministic dynamics where solution cannot be obtained following deterministic methods. Stochastic variability of target concentration or strain type within a period of time gives rise to such dynamics.

Moreover, we observe target cell growth becomes zero with very sharp increase of immune cell for m > 0.25 . This regime can be termed as antibody pattern recognition by immune network shown in Figure 5.

Figure 5

Behavior far away from critical attractor; m ≥ 0.3 , d = 1.0 , r = 0.701 , p = 0.642 , s = 1.72 , c = 1.0 , u = v = 1 , k = 0.265 .

The open region in Figure 5 is marked as a therapeutic intervention case, often recognized by poor/failed immune system, and physics can be explained via the uncorrelated/random response of T cells in lymphocytes.

4 Transitive Lie algebra method to obtain analytic solution structure

In order to study symmetry algebra and its role on dynamics, we plan to obtain invariant Lie symmetry generator based on Lie symmetry algebra in manifold. The dynamical system

will be equivariant under Lie group g on m ( m = 2 )-dimensional manifold M with symmetry group of operator g i , then following must be true under periodic temporal symmetry, i.e.,

for any t ′ > 0 , and the system is called g symmetric.

In this case, g acts locally on m -dimensional manifold M . We are interested in the action of g on p - dimensional sub manifolds N ⊂ M , which we identify as graph of functions in local coordinates. The symmetry and rigidity properties of submanifolds are all governed by their differential invariants. Here, group g acts continuously on differential equation. The model considers target invasion with growth rate r and appropriate immune network response given by functions g ( y ) and f ( x ) . If above dynamical pair of equations (7) and (8) undergo evolution following path of symmetry, it is then necessary to obtain symmetry structure to understand the dynamics. Lie point symmetry and the path of transformation are well known to preserve certain dynamical invariants that can be associated with equations of motion, if they exist. Since Lie compact group of continuous point transformation (continuous group) is special linear group, it can be related to dynamics of holonomic constraints. We assume immune cell responds in a minimal nonlinear way with n = 2 and target stimulation parameter u = v = 1 . With this assumption, we convert coupled autonomous equations into second-order nonlinear ODE in x ( t ) as follows:

which can be expressed as nonlinear ODE

with parameters defined as

with x 0 as initial value of x ( t ) . ODE (equation (18)) is the second-order differential equation, nonlinear in x ˙ and may be linear/nonlinear in x . The idea of linear form of f t ( t , x , x ˙ ) and related to infinitesimal transformation that would entail dynamical extravagance. Once, analytical structure of the solution is obtained, nonlinear term can be added into solution.

Through s l ( 3 , R ) algebra, any two noncommuting, nonproportional symmetry generators should follow

with

For any physical/biological dynamics represented by second-order ODE, differential invariant function of the system is connected to the Lagrangian which describes: A nonsingular Lagrangian admits a symmetry group having dimension ≤ 3 . A nonsingular nth order, n ≥ 2 , admits a symmetry of dimension ≤ n + 3 .

If ℒ ( t , x n ) is corresponding g invariant Lagrangian with nonvanishing Euler-Lagrange expression E ( ℒ ) , then every g -invariant evolution equation should satisfy

where I is an arbitrary differential invariant of the group. If g represents generators of solvable algebra, corresponding differential invariant can be constructed.

Under one-parameter infinitesimal transformation of coordinates (continuous map), we can write

in the neighborhood of identity with vector field defined

for any 0 < ε ≪ 1 . The computation of the flow generated by the vector field is often referred to as exponentiation of the vector field under transitive algebra so that each vector in group g can be integrated through origin in R . It contains an open neighborhood of the origin and flow.

is smooth for any vector field. In case of real field R m , constant vector field defined by g a = ∑ a i ∂ ∂ x i with a = ( a 1 , a 2 , ⋯ a m ) exponentiates to the group of translation

with x ∈ R m . Under continuous symmetry group represented by equations (26)–(30), mono-parametric group generates essential generators for r ≤ 8 . Equation (18) under infinitesimal transformation yields

for all ( t , x , x ˙ ) . Assuming invariance of equation (18), we obtain

Each ξ and η here satisfies the above relation for higher-dimensional algebra. Since, f ( t , x , x ˙ ) is polynomial in x ˙ , it yields set of PDEs in ξ and η . Substituting equations (31) and (32) into equation (35) and separating null coefficients of powers of x ˙ , we obtain;

The general solution of this homogeneous linear system can be formally written as a superposition of linear independent basis solutions ξ a ( t , x ) and η a ( t , x ) with a = 1 , 2 , ⋯ r ≤ 8 following Einstein dummy index

to construct a structure constant. Thus, the symmetry generator takes the form

subject to the Cartan killing condition

where f a b c is the corresponding structure constant.

So, full symmetry group of the differential equation (20) should admit following identities [7]

Commutation relations are regular for regular range of values of a . So, we assume initial data to be regular. Since x = 0 cannot be a regular point, we use initial data at regular point is ( t , x ) = ( t 0 , x 0 ) to evaluate structure constants. This nonsingular initial data will uniquely determine structure constant. In order to represent algebra, we introduce parametrization

According to the above relation, we can adopt

or

Because of Lie-Cartan integrability conditions, killing equations (50)–(59) satisfy Lie algebra. Solving the above equations simultaneously, we find that equation (18) admits sl(3, R ) algebra with the following conditions:

where n presents positive integer. Under transitive algebra, we consider all elements of g = ∑ i c i ∂ ∂ x i and higher-order terms of g for any x i .

So, we propose symmetry algebra admitted by nonlinear ODE (equation (18)) as

for any integer value n . We implement most of the two significant methods of symmetric sub-algebra to integrate the ODE and obtain the analytical structure of the solution. The first method is called method to obtain normal form of generators in the space of variables or quadrature [5]. The second method involves normal form of generators in the space of differential invariant function or first integral function [16]. This method can be used when the number of symmetries is higher or equal to the order of equation. It is necessary to use the generators for the integration procedure in a specific order. This depends on the properties of the algebra. In case of existence of solvable sub algebra of dimension equal to the order of the equation, solution will be obtained in quadrature form if integration should be performed in correct order.

The derived algebra of a Lie algebra ( g , [ ⋅ , ⋅ ] ) is the sub-algebra g 1 of g , defined by

while the derived series is the sequence of Lie sub-algebra defined by g 0 = g and

for any k ∈ N . Such a sequence satisfies g ( k + 1 ) ⊂ g k and the Lie algebra g is said to be solvable if the derived series eventually arrives at the zero sub-algebra. And n -level solvable algebra admits a series of invariant sub-algebras defined by

Derived algebra can be constructed by some linearly independent subset of elements of the commutator of the algebra described by

Because of above relation, the sub-algebra g ( 1 ) is an invariant sub-algebra of g ( 0 ) . Since generators have cardinality l = 8 , they must follow projective group algebra under s l ( 3 , ℛ ) group.

5 Normal form of generators in the space of variables

The first integration method is to reduce the order of the ODE into quadrature form, or equivalently, to find out normal form of generators and use those generators to reduce the order of the equation. Under maximal algebra calculation, three generators follow A 3 , 3 solvable sub-algebra given by

which belongs to Weyl group, semidirect product of time dilation and translations D ⊗ x 2 . Corresponding sub-algebra can be identified as

with

Following [17], we compose sub-algebras as semidirect sums of a one-dimensional sub-algebra and an Abelian ideal with e 1 , e 2 , e 3 are the bases. For a Lie algebra g with corresponding Lie group G = ⟨ exp g ⟩ , the sub-algebra can be constructed as semidirect sum of two algebras. Equations (63)–(66) are two-level solvable with g ( 1 ) = { g 5 , g 7 } and g ( 2 ) = { 0 } . A two-level solvable algebra is designated as follows:

Following Theorem 3 of [16], the first integration method can be repeated n times in the following chain of cosets:

where coset is defined between two derived algebras g ( i ) and g ( j ) as

Here,

forms Abelian algebra. Following Theorem 1 of [16], reduced generators form an algebra with structure constants a subset of the original ones such that

The two generators of the coset act transitive. Using this, transform generators into normal form. By introducing new coordinates ( T , X ) such that T is the independent variable and X is corresponding dependent variable in normal form. As T = f ( t , x ) and X = g ( t , x ) , we employ

two-dimensional sub-algebra such that

holds true. Using twisted Goursat algorithm for decomposable algebra [17], we obtain a = 1 and b = 3 , which yields

In terms of normal variable (canonical), we obtain ODE

which is in quadrature form with linear form of solution. Inverse mapping of variables yields solution structure of target component as

with C 1 and C 2 to be determined from initial conditions.

6 Group invariant solution structure

This method uses differential invariant of the subgroup to obtain solution. We use Abelian solvable sub-algebra

under prolonged group to evaluate invariant differential function in the evolution dynamics. We use g ≡ g 7 here.

with

with characteristic equation

which provides two invariant functions, namely ϕ ( t , x ) and ψ ( t , x , x ˙ ) with

We concentrate on the system where number of equations in the system is same as number of dependent variables. That means

Two invariant functions give relation

Implementing equations (80)–(89), we obtain

where C 1 and C 2 will be determined from initial conditions.

7 Numerical simulation and results

In the model, parameters u and v play a critical role in the transition of the dynamics. So, we expand solution around linear (analytical) value by adding nonlinear part in a perturbative way, i.e.,

Parameter u act as immune stimulation and v as correlation parameter in the network. We keep parameter values same as obtained in dynamical simulation. Following is the result of simulation. Linear behavior is exhibited in the simulation using analytical solution without perturbation (Figure 6). Simulation data from linear structure of the solution exhibit the presence of local critical attractor in stable manifold.

Figure 6

Linear solution behavior based on analytic structure; d = 1.0 , r = 0.701 , p = 0.642 , s = 1.23 , m = ‒ 0.14 , c = 1.0 , k = 0.93 , u = v = 0 .

Setting u = 1 , v = 1 and second order perturbation, growth pattern exhibits the presence of limit cycle determined by parameter values d = 1.0 , r = 0.7 , p = 0.642 , c = 1.0 , k = 0.25 (Figure 7). The negative of ( m = − 0.142 ) designates attractive phase of dynamics with stable oscillation in closed manifold. The target growth shows consistent increase with numerous small oscillations, whereas immune response indicates strong oscillation between 35 and 80 h. In the simualtion, time can be hours or days depending on the intensity of infection. In the first phase, negative oscillatory values of immune density indicate not responding to annihilate target or failure to respond in target invasion. The way immune responds in any case of target invasion in the body is through antibody protein production via CD8+ or CD4+ cytokinase protein. The process of the production of these proteins involves multistep interaction via several meta stages.

Figure 7

Dynamical behavior with weak correlation u = v = 1 ; d = 1.0 , r = 0.701 , p = 0.642 , s = 1.23 , m = ‒ 0.142 , c = 1.0 , k = 0.93 , u = v = 1 .

Values u = 1 and v = 3 indicate strong correlation between immune and target dynamics within local domain (Figure 8). Our simulation reflects chaotic correlated dynamics between two components using second-order perturbation. This shows the role of parameter v as correlation in nonlinear growth dynamics.

Figure 8

Chaotic behavior of strongly correlated dynamics; d = 1.0 , r = 0.701 , p = 0.642 , s = 1.634 , c = 1.0 , k = 0.15 , m = ‒ 0.13 .

Under second-order perturbation, pattern behavior shifts away from critical attractor designated by positive value of m ( m > 0.25 ). As a result, phase trajectory becomes open manifold as shown in Figure 9 (therapeutic intervention phase). Existence of open manifold case (Figure 9) with low immune cell growth indicates that disease persists and requires therapeutic intervention. This region is also viable for adaptive immunity in case of lethal infection. The parameter value m ≥ 0.25 plays significant role in the development of adaptive immunity, which is noted by certain T-cell proliferation through antigen production. This feature is prominent in both approaches. The parameter m here plays significant role via antigen production in case of adaptive immunity. The existence of continuous oscillation or unstable dynamics is noted far away from stable equilibrium ( k = 3.4189 ) when we increase death rate of immune cell d = 1.18 . When we set parameter values v = 2.0 , m = 0.3 , second-order perturbation simulation of target growth exhibits steady positive for longer time period with maximum value around 40 h. Positive value of m drives the system away from closed manifold. Immune growth pattern within 30 h of infection cannot compete with target and then starts target annihilation process between 40 and 80 h after first infection.

Figure 9

Pattern behavior far away from critical attractor; m ≥ 0.25 , d = 1 , r = 0.701 , p = 0.652 , s = 1.625 , c = 1.0 , u = 1.0 , v = 3.0 , k = 0.254 .