BISTABILITY ANALYSIS OF VIRUS INFECTION MODELS WITH TIME DELAYS

. Mathematical models with time delays are widely used to analyze the mechanisms of the immune response to virus infections and predict vari- ous therapeutic eﬀects. Using the lymphocytic choriomeningitis virus infection model as an example, this work describes an original computational technol- ogy for searching the bistable regimes of such models. This technology includes numerical methods for ﬁnding all possible steady states at ﬁxed values of pa- rameters, for tracing these states along the parameters and for analyzing their stability.


1.
Introduction. Mathematical modelling of infectious diseases is a rapidly expanding field at the interface of applied mathematics, immunology and virology [22,23,3]. The formulation and analysis of the mathematical models of virus infections aims at a mechanistic understanding of the pathogenesis mechanisms of disease dynamics and the ways to correct the offset dynamics. The mechanisms regulating the development and outcome of human infections are studied in clinics and using experimental animal model systems. One of the broadly used experimental infection models is the virus infection of mice with lymphocytic choriomeningitis virus (LCMV).
The experimental LCMV infection of mice is a fundamental model system of modern immunology that enabled the discovery of key mechanisms of the development of immunopathological processes and chronic infections. The underlying processes were characterized in terms of the systems dynamics of the virus and immune responses (e.g., anergy, exhaustion of T cells) [2]. The phenotypes of the LCMV infection dynamics were quantitatively examined and it was shown that the development and outcome of disease depend on the race between the virus infection spreading and the strength of the immune response. A thorough understanding of the processes underlying the interactions of the virus and the host organism for LCMV infection enabled the development of biologically relevant mathematical models formulated with delay differential equations [4,2,3]. The models can be used to predict the sensitivity of the infection dynamics to antiviral and immunomodulatory treatments. In turn, this understanding can be utilized for investigation of optimal control of the dynamics of chronic virus infections in humans such as hepatitis B or human immunodeficiency virus type 1 (HIV) infections which share a number of common pathogenesis mechanisms with the LCMV infection.
The dynamics of virus infections can take various phenotypes, including high-and low-viral load infections [24,13]. The high viral load infection represents a chronic disease with a low level immune response. A low viral load infection corresponds to the state of a recovered healthy organism with immunological memory maintained by an ongoing low level antigenic stimulation, or to a less severe form of chronic disease resulting from a stronger immune response which is however not efficient to completely eliminate the infection. These two extreme variants of disease dynamics can be represented by the bistability feature of the mathematical model of infection suggesting that there co-exist at least two stable steady states for given values of model parameters (see, e.g. [23,12]).
Bistability is a fundamental property of complex biological systems [1,16] to which the virus infections belong to. It determines the switches between different modes of disease dynamics, hysteresis phenomena and the emergence of oscillations. The examination of bistability properties of mathematical models represents a challenging computational problem. It relies heavily on the prediction of the existence of bistability domains in the model parameter space. Understanding it is most important for elaboration of an optimal treatment strategy of the offset system dynamics. Indeed, a bistable dynamical system can be transferred to a favourable steady state using the optimal disturbance approach as outlined in [5,6,7]. In the situation of monostability, the transfer of the system away from unwanted steady state requires the implementation of control taking the system to a neighborhood of a given trajectory (e.g., feedback stabilization, extremal shift or the optimal programme (open loop) control [8]). An alternative solution could be a mixed control consisting of a parametric shift of the system to a bistable domain of the parameter space and further correction of its dynamics by optimal disturbances of the steady state. These issues remain poorly investigated for mathematical models based on delay differential equations.
An efficient treatment of the above control problems for virus infection requires the following tools of analysis of the mathematical models: • methods for computing all steady states of the model and tracing them as functions of model parameters, • numerical methods for analyzing the stability of the steady states, • numerical methods for constructing optimal disturbances bringing the dynamical system from one stable steady state to another.
In this paper, we propose an original technology for analyzing the bistability of virus infection models with delayed arguments. The technology includes methods which allow one to compute all steady states for fixed parameter values, to trace them along parameters, and to analyze their stability. Methods for constructing the optimal disturbances were proposed and justified in [5,6,7,25].
The models of virus infections we refer to have the following structural elements: the logistic description of virus growth, Lotka-Volterra type of virus-immune system interaction with a bell-shaped functional response, confined exponential equation for immune cell homeostasis, and time delay of immune reaction. The proposed technology is particularly designed for virus infection models with such properties, features guaranteed computation of all steady states in the considered region of parameters and applies an algebraic approach for the analysis of their stability. When studying stability, the leading eigenvalues of the corresponding nonlinear eigenproblem are not traced along the model parameters but are computed for each new parameter value by solving the complete eigenvalue problem for a rational approximation of the original nonlinear matrix pencil and refining the computed approximate leading eigenvalues by a local method. This ensures that the eigenvalue with the maximum real part is always correctly computed.
It should be noted that currently there are numerical software packages DDE-BIFTOOL [15,27], knut (formerly PDDE-CONT) [26] and TRACE-DDE [11], which implement various methods for steady states computation and their stability analysis for time delay systems of more general forms than models of virus infections (e.g., systems with state-dependent delays). The first package is designed for bifurcation analysis of time delay systems. It allows one the continuation of steady states along parameters starting from specified initial steady state solution. However, the problem of guaranteed computation of all steady states is not considered. The eigenproblems are solved by the subspace iteration based on solving the initial value problems for the corresponding linearized equations. The second package is designed for analysis of periodic solutions. Steady states are treated as constant periodic solutions. The third package is not designed for searching or tracing of steady states. It allows one to find the stability regions of a given steady state by two variable parameters. To solve the eigenvalue problems, an approximation of the infinitesimal generator for the linearized equations by the collocation method on the Chebyshev grid is used. We do not exclude that on the basis of the methods cited above, it is possible to develop another technology that also automatically (without handwork) computes all steady states of viral infection models with delay, as functions of model parameters, and investigate their stability, but today the only technology that solves the above problem is the technology which we propose here.
The paper is organized as follows. Section 2 briefly describes the LCMV infection model under consideration. In Sections 3 and 4, the proposed methods for computing steady states and analyzing their stability are described and justified. Section 5 demonstrates the application of the proposed methods. Section 6 describes the proposed bistability analysis technology. Section 7 summarizes this work and discusses a wider application of the developed methods.
2. Mathematical model of LCMV infection. The basic mathematical model of LCMV infection in mice proposed and analysed in [4] is formulated as the following system of non-linear delay differential equations: where g p (W ) = 1/(1 + W/θ p ) 2 , g e (W ) = 1/(1 + W/θ E ) 2 . This system describes the dynamics of the following time-dependent variables: concentration of viruses V (t), population densities of two LCMV-specific cytotoxic lymphocytes (CTLs)precursors E p (t) and effectors E e (t), and the cumulative viral load W (t). The biological meaning of parameters is explained in Table 1. For the biological adequacy of the model (1), the values of all these parameters must be non-negative. Moreover, of interest is only the case when the values of all the parameters are positive, and 0 < τ ≤ τ A . Further we will assume that these conditions are fulfilled.
To determine the solution of system (1) for t > 0, it is necessary and sufficient to define the following initial functions V (t) for −τ A ≤ t ≤ 0, E p (t) for −τ ≤ t ≤ 0, and the values E e (0) and W (0). However, we will assume for the sake of generality that the initial conditions for all variables are specified on −τ A ≤ t ≤ 0.
The initial value problem for system (1) with non-negative initial conditions and non-negative parameters has a unique non-negative solution in any finite time interval [0, T ]. This can be proved using the technique described in [22], which is based on the Bellman's method of steps and makes use of a linear ordinary differential equations system majorizing the right hand side of the DDE system.
Let us denote the vector of system (1) variables as to express this system in the following compact form: As mentioned above, we assume that U (t) is defined on −τ A ≤ t ≤ 0.
3. Computation of steady states. For an arbitrary system (2), algorithms that find all the steady states are not currently known. However, for the system (1)  (1) with fixed values of the parameters are non-negative solutions of the following non-linear system: From the first equation it follows that the necessary condition for the solution to be non-negative is that V belongs to the closed interval 0 ≤ V ≤ V mvc . It is easy to verify that the system (3) has a solution and there are no other solutions with V = 0. If V > 0, then the corresponding value of E e is uniquely found from the first equation, E p from the third equation, W from the fourth one. Thus, for V > 0 any solution of (3) has the following form:

Moreover, (4) is a solution of (3) if and only if it satisfies the second equation in (3). That is
Therefore, the computation of the steady states of (1) reduces to the computation of zeros of the continuous scalar function and for any such V 0 the equation (5) has at least one root in the closed interval It can be shown that (6) is fulfilled for For all the parameter values considered below this formula gives V 0 ≥ 4.5 and An algorithm proposed in [10] allows one to find all zeros of a given continuous function f (x) in a given closed interval [x 1 , x 2 ] with a given accuracy in the case when the number of zeros is finite. This algorithm uses the standard fzero and fmin procedures [17,21] included in many numerical software packages. The first of these procedures, based on a combination of the bisection method with the secant method, allows one to compute an approximate zero of the function f (x) in a given closed interval [x 1 , x 2 ] with a given accuracy ρ in the case when the function takes nonzero values of different signs at the ends of the interval, that is, the procedure finds a point The second procedure, based on a combination of the golden section search method with the parabolic method, allows one to find the minimum of a given continuous function in a given interval.
We will use the algorithm proposed in [10] to find all solutions of (5) with a given accuracy ρ in the interval [V 0 , V mvc ], with V 0 given in (7). These solutions will be a complete set of positive steady state values of V . The corresponding values of the other variables can be computed using the formulas given in (4).

3.2.
The dependence of steady states on a parameter. A numerical study of the dependence of steady states of (2) on the parameters is conveniently performed by varying one of the parameters and fixing the values of the others. Denote the variable parameter by s, and the interval of its variation by σ (this interval is assumed to be open). We denote by Y (s, V ) the right side of equation (5) as a function of the selected parameter s and the variable V . In this case, as before, we assume that the other variables depend on s and V according to the formulas (4).
Since Y (s, V ) is a rational function of the variables s and V , having no poles in the domain σ × (0, V mvc ), the set is a part of a plane real algebraic curve [20]. For definiteness, we assume that this curve is simple, that is, it does not have multiple irreducible components. Then equality ∂Y ∂V (s, V ) = 0 (8) holds only at a finite number of points.
The task is to break the curve A into elementary components of the form where V (s) is a single-valued real analytic function of the variable s defined in a certain interval σ C ⊂ σ, and points at which the equality (8)  . . , m k } with the method described in Section 3.1. We assume that the grid step size along s is chosen so small that each of the smooth elementary components of the form (9) that need to be constructed contains at least one point of the form (s k , V k j ). In addition, we assume that at all such points the equality (8) does not hold.
To construct smooth elementary components we need the operation for refining a given row Z = ( Z 1 , . . . , Z m ), whose elements are real and pairwise different, by a given real set Z = {Z 1 , . . . , Z m } whose elements are also pairwise different, where m ≥ m. This operation consists in replacing the row Z with the row (Z i1 , . . . , Z i m ) which consists of m pairwise different elements of Z and minimizes In addition, we need the standard ode45 numerical integration procedure, based on explicit 4th and 5th order Runge-Kutta methods. We will use this procedure for tracing steady states along a given parameter by numerically solving the initial value problem of the form We will stop the integration if the variable s reaches the specified final value, or if during the integration process it is found that the function f (s) = ∂Y /∂V (s, V (s)) vanishes on the integration trajectory. Note that the ode45 procedure chooses the integration grid adaptively depending on the accuracy requirements. To test the second stopping criterion, stepping from the node s i to the node s i+1 of this grid, the ode45 procedure checks whether the function f (s) accepts values of different signs in these nodes and, if so, finds zero s = s * between these nodes using the fzero procedure mentioned in the previous section. For each value of s, the corresponding value of V (s), which is necessary for computing f (s), is computed by one step of numerical integration from the node s i . After this zero is found, the ode45 procedure stops, and the point (s * , V (s * )) is considered as the last point of the integration trajectory.
At the first step of the algorithm for constructing elementary components, we consider the interval s 1 ≤ s ≤ s 2 . Using the ode45 procedure, we solve the initial value problem (10) with V 0 = V 1 j , s 0 = s 1 and the final value s = s 2 for all j = 1, . . . , m 1 , and form two-point rows where (s 1 j , V 1 j ) means the last point of the integration trajectory of the j-th initial value problem. Note thats 1 j < s 2 if the second stopping criterion was satisfied when solving the j-th initial value problem, ands 1 j = s 2 otherwise. Without loss of generality, we assume that the second stopping criterion was satisfied for the first m 1 initial value problems and was not satisfied for the remaining m 1 − m 1 . Refine the row ( V 1 m 1 +1 , . . . , V 1 m 1 ) by the set V(s 2 ), and replace in the second elements of the rows B m 1 +1 , . . . , B m 1 the numbers V 1 j by their refined values. As a result, we obtain the following two-point rows of two types: where s 1 <s 1 j < s 2 , and If m 1 − m 1 < m 2 then m 2 = m 2 − m 1 + m 1 elements of the set V(s 2 ) were not used for refinement. Without loss of generality, we assume that these elements are V 2 1 , . . . , V 2 m 2 . We solve numerically the initial value problem (10) with V 0 = V 2 j , s 0 = s 2 and the final value s = s 1 for all j = 1, . . . , m 2 , and form two-point rows where (s 2 j , V 2 j ) means the last point of the integration trajectory. Note that when solving each of these initial value problems, the second stopping criterion will be satisfied, that is s 2 >s 2 j > s 1 . This is where the first step of the proposed algorithm for constructing elementary components, consisting in the formation of two-point rows (11), (12) and (13), ends. Each of the obtained rows consists of the initial and final points of an elementary component in the interval s 1 ≤ s ≤ s 2 . In addition, each point of s 1 × V(s 1 ) is the starting point of one of these rows, and each point of s 2 × V(s 2 ) is the end point of one of these rows. Figure 1 illustrates the first step of constructing elementary components in the case when m 1 = m 2 = 4, m 1 = m 2 = 2. The circles are the initial points of numerical integration, the squares are the final ones. Thick lines connect the end points of numerical integration with their refinements. As a result, we obtained the following two two-point sets of form (11):  Further, in exactly the same way, we construct two-point rows for the intervals s k ≤ s ≤ s k+1 , k = 2, . . . , N − 1. This completes the construction of two-point rows and the following sewing procedure is applied to the resulting rows: rows (α 1 , . . . , α k ) and (β 1 , . . . , β l ) are replaced with one row (α 1 , . . . , α k , β 2 , . . . , β l ) if α k = β 1 and at this point the equality (8) does not hold. The procedure is applied as long as there are rows that can be sewed together. As a result, we obtain a set of rows of the form C = ((s 1 , V 1 ), . . . , (s r , V r )) , where 2 ≤ r ≤ N , s 1 ≤ s 1 < . . . < s r ≤ s N . Each of these rows is a row of points with high accuracy located on a separate elementary component of the form (9), that is, V j ≈ V (s j ). Constructing a spline through these points we obtain an approximation of this component.
Recall that the algorithm described above enables one to find all the smooth elementary components of the form (9) of the curve A if the grid s 1 < . . . < s N is chosen such that each of these components contains at least one point of the form (s k , V k j ) and at all such points the equality (8) does not hold. Taking in to account that the number of elementary components is finite and the number of points of the curve A in which the equality (8) holds is finite as well, it is easy to satisfy these conditions by choosing a sufficiently small grid step size along s and slightly shifting the node s k if the equality (8) holds at the point (s k , V k j ).

Stability analysis.
To study the stability of a given steady state U of model (2) with fixed values of the parameters, we write the solution near this steady state in the form U (t) = U +εU (t)+O(ε 2 ) where ε is a small parameter in absolute value. Substituting this solution into (2) and assuming that the obtained equalities hold for all arbitrarily small values of ε, for the function U (t) we obtain the following system of linear differential equations: where are constant square matrices of order 4. The steady state U is asymptotically stable if any solution of the system (15) of the form where λ is a complex number, and U is a constant 4-component non-zero vector, monotonically decreases as t → +∞, that is when the real part of any such λ is negative. Substituting (16) into (15), we obtain the following nonlinear eigenvalue problem: where A(λ) = λI−L 0 −e −λτ L τ −e −λτ A L τ A , and I means the identity matrix of order 4. Thus, the study of the asymptotic stability of a steady state U reduces to solving the nonlinear eigenvalue problem (17) and checking if all the found eigenvalues lie strictly in the left half-plane. Algorithms that solve the complete nonlinear eigenvalue problem of the form (17) are not known. Moreover, this problem generally has an infinite number of eigenvalues. Therefore, we will approximate this problem with a rational eigenvalue problem, which, in turn, can easily be reduced to a polynomial one. Using the leading eigenvalues of the polynomial eigenproblem we will compute approximate leading eigenvalues of the original eigenproblem (17) which can then be refined by any local Newton-type method.
We introduce a new spectral parameter µ = e δλ , where δ is a given sufficiently small positive number. Using the formula that approximates the first derivative in the implicit second order scheme BDF2 [19], we obtain the following approximation: The non-linear eigenproblem (17) can be approximated by the following rational eigenproblem: Note that along with the value of δ, the scheme used for approximating the spectral parameter λ is a parameter of the algorithm. If necessary, one can use a scheme of higher order of accuracy. Multiplying (18) by µ m A this equation can be reduced to the following polynomial eigenproblem: By supplementing the eigenproblem with the identities µ j U = µ j U , j = 0, . . . , m A − 1, we obtain a system which can be written in the following form: Here X i = µ m A −i U , and M is block matrix of block order m A with blocks of order 4. All blocks of this matrix are zero, except for subdiagonal blocks M j+1,j = I (j = 1, . . . , m A − 1) and four blocks M 11 = C 1 , M 12 = C 2 , M 1m = C m and M 1m A = C m A locating in the first block row. Thus, we have reduced the approximate solution of the nonlinear eigenproblem (17) to the computation of the eigenvalues of matrix M . To compute these eigenvalues one can use the standard QR-algorithm [18]. Approximate eigenvalues of the problem (17) are expressed in terms of the eigenvalues of matrix M as follows: λ = ln(µ)/δ.
To analyze the stability of a steady state, it is necessary to find the eigenvalue of the problem (17) with the maximum real part. The algorithm described above allows us to find approximate eigenvalues of the problem (17). To increase reliability, we will choose a certain number of leading eigenvalues, i.e., those with maximum real parts, and refine each of them using the method of successive linear problems [14], which consists in the following iterative process: Here, j is the iteration number, and ν j is the maximal in absolute value eigenvalue of the linear matrix pencil νB(λ j ) − A(λ j ) where In conclusion, we note that the order of matrix M in (19) is inversely proportional to δ and proportional to the number of equations in the model under consideration. If it turns out to be too large for solving the complete eigenvalue problem in a reasonable time, then we can solve a partial eigenvalue problem computing for reliability a dozen eigenvalues which are maximum in absolute value.
5. An example. To demonstrate the use of the proposed algorithms for computing the steady states of the model (1) and studying their stability let us consider the values of parameters as specified in Table 5 from [4] and choose the accuracy ρ for computing the steady states equal to 10 −14 , the number of grid nodes N along b p equal to 100, and the value of δ equal to 10 −2 (which leads to a matrix M of order 2240). The parameter b p = 7.73 · 10 −5 will be varied in the range from 10 −5 to 10 −3 .
As a result of using the algorithms described in Section 3, four different rows of the form (14) were found. The elementary components built with these rows are shown in the upper left part of Figure 2. Solid lines indicate stable steady states and dashed lines indicate unstable ones. The stability of each row was studied separately for each s k . It can be seen that for 10 −5 ≤ b p ≤ b * p , where b * p ≈ 6.7 · 10 −4 , there are four steady states: with V = 0, V ≈ 10 2 , V ≈ 10 7 and V = V * ≈ 4.82 · 10 7 , which will be referred to as states I, II, III, and IV, respectively. At the point (b * p , V * ), the equality (8) Figure 3 illustrates the dependence of the leading eigenvalues in each of the steady states on the parameter b p . It shows 15 -16 leading (with the largest real parts) eigenvalues for steady states I, II at b p = 10 −5 , 5.1 · 10 −4 and 10 −3 and for steady states III, IV at b p = 10 −5 , 3.5 · 10 −4 , 6.7 · 10 −4 . Markers "x", "+" and "o" denote eigenvalues at the first, second and third parameter values, respectively. One can see that there are only four different eigenvalues in state I. These multiple eigenvalues do not change when the parameter is varied. The real part of the rightmost of these eigenvalues is greater than zero and therefore state I is not asymptotically stable for all considered values of the parameter.
In state II, the four rightmost eigenvalues change little when the parameter is varied. The remaining eigenvalues shift to the left as the value of the parameter increases. The real part of the rightmost eigenvalue is greater than zero, therefore, state II is also not asymptotically stable for all considered values of the parameter.
In state III, the leading eigenvalues in general shift to the left as the parameter value increases, but the real part of the rightmost eigenvalue remains greater than zero, so this steady state is also not asymptotically stable.
In state IV, as the value of the parameter increases, the eigenvalues first shift to the left and then to the right. At the initial value of the parameter b p , the real part of the rightmost eigenvalue is less than zero. As the value of the parameter increases, this eigenvalue does not change, but the second eigenvalue shift to the right and it becomes the rightmost one at b p ≈ 6·10 −4 . Moreover, for b p = 6.7·10 −4 , this eigenvalue already lies in the right half-plane. Thus, state IV is stable at all considered values of the parameter, except the last one. Note that at the last value of the parameter, the eigenvalues in states III and IV visually coincide, just like these states themselves. . Leading eigenvalues in steady states I, II, III and IV at b p = 10 −5 ("x"), 3.5 · 10 −4 ("+"), 6.7 · 10 −4 ("o").
6. Analysis of bistability. In Section 5 we considered the case when only one of the four existing steady states was stable, the state with the highest viral load. For systems describing the dynamics of virus infections, it is of primary importance from the point of choosing the appropriate control strategy to identify the domains in the parameter space in which two or more steady states of the model can be simultaneously stable.
For the reference values of the model parameters presented in [4] (Table 5) except for parameter β = 3.35 (virus growth rate) which was varied in the range from 10 −3 to 5, the steady states and their stability were analyzed with the results shown in Figure 4. It follows that the system has four steady states with V = 0, V ≈ 1.3 · 10 2 , V ≈ 2.3 · 10 7 and V ≈ 4.82 · 10 7 , which we number as I, II, III and IV, respectively. Steady state IV is asymptotically stable for all values of the varied parameter β from the specified range, whereas steady state II is stable only in the interval 10 −3 ≤ β < β * ≈ 1.718. Therefore, the system is bistable for the values of parameter β between 10 −3 and β * and the rest model parameters being equal to the values specified in Table 5 [4]. As bistability is a relevant property of normal feedback regulation of functioning physiological systems, it is important to examine computationally whether a nonlinear dynamical system being monostable in some parameter domain can be transferred to a bistable functioning mode. One possibility is based on changing the parameters of the underlying regulatory processes. The here developed computational methods allow one to search for and implement the required parametric shifts. As a practical example, consider the set of parameter values for which the model is monostable (Figure 2). By a sequential decrease of the virus growth rate β (as shown in Figure 5) the system can be transferred to the parameters domain where it is bistable. The values of viral load V in the steady states shown in Figure  5 were obtained for the same values of parameters considered in Figure 2, except for parameter β which was varied in some range.
Overall, we have deduced that the nonlinear mathematical model can be either monostable or bistable depending on the model parameter values and the system can be transferred to a bistable mode by a monotone variation of control parameters. The possibility of tracing the stability of the second steady state provides the tool for quantifying the required parameter shift of the system such that target criteria for the second steady state are met. For the scenario shown in Figure 5, if the aim is to maximally reduce the equilibrium viral load then the parametric shift should be the one presented in Figure 5D. 7. Conclusions. In this study by the example of the model of experimental LCMV infection we proposed and justified novel algorithms for computation and stability analysis of steady states of nonlinear mathematical models with time delays which are extensively used in mathematical immunology for modeling the dynamics of infectious diseases. The algorithms permit the detection of the mono-or bistability. These are key properties of dynamical systems for gaining insight into the kinetic regulation of the system behavior and the development of virus infections. We showed that the parameter space of the model can be subdivided into regions of mono-and bistability. Using parametric shift of the processes underlying the virus-host organism interaction, the infection can be transferred from one mode of functioning to another. The shift of the system to a bistable domain of parameter space offers a possibility to control the phenotype of infectious disease through the use of optimal disturbances of the unfavorable stable steady state. However, if the nonlinear system can't be transferred to a bistable mode with a favorable stable steady state, the control strategy is to be a life-long feedback stabilization or programme control (Pontryagin's maximum principle) which represent more complex tasks.
The computational technology of bistability analysis developed in our study makes it possible • to study the mechanisms of emergence of various dynamics of human virus infections (e.g., HIV and hepatitis) and the animal infection with LCMV; • to identify the potential targets for therapies at the level of immunophysiological processes; • to construct effective control strategies for infectious disease in the framework of phenotypically-based or personalized types of combination therapies.
This technology is applicable for any time delay model used in mathematical immunology except for the algorithm described in Section 3.1. For example, for the model of hepatitis B virus infection [23,9], the computation of all steady states for fixed parameter values is not reduced to solving a single scalar equation for viral load (as in the case of LCMV infection model considered in this paper) but to solving a system of two scalar equations for two model variables and requires the joint use of the algorithm from Section 3.1 and the tracing algorithm described in Section 3.2. At the same time, the methods for tracing of the steady states along the model parameters and for their stability analysis do not require any significant modifications.