Analysis and Control of the Inflammatory Immune Response Model

Model equations [12]

The variables PV, NV, dv, and CV represent the bacterial pathogen population pro-inflammatory cytokines, tissue damage, and the anti-inflammatory mediators. Up and ua are the control parameters, which are the pro-inflammatory and the anti-inflammatory therapy of the system.

The model equations are

$\begin{array}{l}fn6=\frac{n{v}^6}{(1+{(\frac{cv}{c\inf })}^2)}\\ fnvdv=\frac{(nv+kcnd(dv))}{(1+{(\frac{cv}{c\inf })}^2)}\\ rv=\frac{(knp(pv))+(knn(nv))+(knd(dv))}{(1+{(\frac{cv}{c\inf })}^2)}\text{ (1)}\\ fn=\frac{nv}{(1+{(\frac{cv}{c\inf })}^2)}\end{array}$

$\begin{array}{l}\frac{d\left(pv\right)}{dt}=kpg(pv)\left(1-\left(\frac{pv}{p\inf }\right)\right)-\frac{\left(kpm(sm)pv\right)}{\left(\mu m+\left(kmp(pv)\right)\right)}-kpn(fn)pv\text{ (2)}\\ \frac{d\left(nv\right)}{dt}=\frac{snr(rv)}{\left(\mu nr+rv\right)}-\mu n(nv)+up\\ \frac{d(dv)}{dt}=\frac{fn6}{(xd{n}^6+fn6)}-\mu d(dv)\\ \frac{d\left(cv\right)}{dt}=sc+\frac{kcn(fnvdv)}{\left(1+fnvdv\right)}-\mu c(cv)+ua\end{array}$

The base parameter values are

kpm = 0.6; kmp = 0.01; sm = 0.005; kpg = 1;pinf = 20e+06; kpn = 1.8;knp = 0.1;knn = 0.01;snr = 0.08;

knd = 0.35; xdn = 0.06; cinf = 0.28;sc = 0.0125;kcn = 0.04;kcnd = 48; up = 0; ua = 0;

$\mu nr=0.12;\mu n=0.05;\mu c=0.1;\mu d=0.02;\mu m=\mathrm{0.002.}$

kpg is the pathogen growth rate, knp represents the activation of phagocytes by pathogen, Kpn represents the effect of pathogens on phagocytes, kpm and kmp are monod rate constants that describe the reduction of phagocytes, kcnd controls the effectiveness of activated phagocytes versus damage in the production of the anti-inflammatory mediator, while knn represents the activation of phagocytes by already activated phagocytes.

Bifurcation analysis

Bifurcation analysis is performed using the MATLAB software MATCONT, which locates branch points limit points and Hopf bifurcation points [13,14]. Consider a set of ordinary differential equations

$\frac{dx}{dt}=f(x,\alpha )\text{ (3)}$

x ∈ Rⁿ With a bifurcation parameter be α. Since the gradient is orthogonal to the tangent vector,

The tangent $z=[z_1,z_2,z_3,z_4,\mathrm{....}{z}_{n+1}]$ must satisfy

Az = 0 (4)

A is given by

$A=[\partial f/\partial x|\partial f/\partial \alpha ]\text{ (5)}$

Where is the Jacobian matrix? For both limit and branch points, the matrix [∂ƒ / ∂x] must be singular. The n + 1th component of the tangent vector Z_n+1 = 0 for a limit point (LP) and for a branch point (BP) the matrix must be singular. At a Hopf bifurcation point,

$\det (2f_x(x,\alpha )@I_n)=0\text{ (6)}$

@ indicates the bialternate product and I_n is the n-square identity matrix. Hopf bifurcations cause limit cycles and should be eliminated because limit cycles make optimization and control tasks very difficult. More details can be found in Kuznetsov [15,16] and Govaerts [17].

Multiobjective Nonlinear Model Predictive Control (MNLMPC)

The procedure developed by Flores Tlacuahuaz, et al. [18] is used for performing the MNLMPC calculations. Let the objective function variables ${\displaystyle \sum _{t_{i=0}}^{t_i=t_f}{q}_j(t_i)}$ (j = 2.n) for a problem involving a set of ODE

$\frac{dx}{dt}=F(x,u)\text{ (7)}$

Where t_f is the final time value, and n the total number of objective variables and u the control parameter? First, the single objective optimal control problem, independently and individually optimizing each of the variables ${\displaystyle \sum _{t_{i=0}}^{t_i=t_f}{q}_j(t_i)}$ is solved. Leading to the values $q_j^{*}$ . Then the multiobjective optimal control (MOOC) optimization problem that will be solved is

$\begin{array}{l}\min ({\displaystyle \sum _{j=1}^n({\displaystyle \sum _{t_{i=0}}^{t_i=t_f}{q}_j(t_i)}-q_j^{*}){)}^2}\text{ (8)}\\ subjectto\frac{dx}{dt}=F(x,u);\end{array}$

This will provide the values of u at various times. The first obtained control value of u is implemented, and the rest are discarded. This procedure is repeated until the implemented and the first obtained control values are the same or if the Utopia point where ( ${\displaystyle \sum _{t_{i=0}}^{t_i=t_f}{q}_j(t_i)}=q_j^{*}$ for all j) is obtained.

Pyomo [19] is used for these calculations. Here, the differential equations are converted to a Nonlinear Program (NLP) using the orthogonal collocation method. The NLP is solved using IPOPT [20] and confirmed as a global solution with BARON [21].

Sridhar [22] proved that the MNLMPC calculations to converge to the Utopia solution when the bifurcation analysis revealed the presence of limit and branch points. This was done by imposing the singularity condition on the co-state equation [23]. This makes the constrained problem an unconstrained optimization problem, and the only solution is the Utopia solution. More details can be found in Sridhar [23].

Results

When up was the bifurcation parameter, a limit point was found at (pv,nv,dv,cv,up) values ((75.123799, 2.493135, 50, 0.524251, 0.046348) (Figure 1a).

Figure 1a: Bifurcation diagram (with up as bifurcation parameter).

When ua is bifurcation parameter, two limit points were found at (pv,nv,dv,cv,ua) values of((0.100172, 0.226022, 40.042712, 7.458674, 0.704163); (0.10125, 0.019855, 0.007619, 0.772682, 0.063055) (Figure 1b).

Figure 1b: Bifurcation diagram (with ua as bifurcation parameter).

When kpg was the bifurcation parameter, a branch point was located at (pv,nv,dv,cv,kpg) values of (0.000000, 1.552090, 49.999999, 0.524251, 2.120064) (Figure 1c) .

Figure 1c: Bifurcation diagram (with kpg as bifurcation parameter).

For the MNLMPC ua and up are the control parameters, and ${\displaystyle \sum _{t_{i=0}}^{t_i=t_f}pv(t_i)},{\displaystyle \sum _{t_{i=0}}^{t_i=t_f}dv(t_i)}$ were minimized individually, and led to values of 0 and 10.6686069466. The overall optimal control problem will involve the minimization of ${({\displaystyle \sum _{t_{i=0}}^{t_i=t_f}pv(t_i)-0})}^2+{({\displaystyle \sum _{t_{i=0}}^{t_i=t_f}yv(t_i)}-\text{10}\text{.6686069466})}^2$ was minimized subject to the equations governing the model. This led to a value of zero (the Utopia point). The MNLMPC values of the control variables, ua and up were 1.504 and 1.472. The MNLMPC profiles are shown in Figures 2a-2c. The control profiles of ua and up exhibited noise and this was remedied using the Savitzky-Golay filter to produce the smooth profiles uasg; upsg (Figure 2d,2e).

Figure 2a: (MNLMC pv vs. t).

Figure 2b: (MNLMC dv vs. t).

Figure 2c: (MNLMC ua uasg vs. t).

Figure 2d: (MNLMC up upsg vs. t).

Figure 2e: (MNLMC cv nv vs. t).

Discussion of results

Theorem

If one of the functions in a dynamic system is separable into two distinct functions, a branch point singularity will occur in the system.

Proof

Consider a system of equations

$\frac{dx}{dt}=f(x,\alpha )\text{ (9)}$

x ∈ R_n. Defining the matrix A as

$A=\left[\begin{array}{l}\frac{\partial f_1}{\partial x_1}\frac{\partial f_1}{\partial x_2}\frac{\partial f_1}{\partial x_3}\frac{\partial f_1}{\partial x_4}\mathrm{..........}\frac{\partial f_1}{\partial x_n}\frac{\partial f_1}{\partial \alpha }\\ \frac{\partial f_2}{\partial x_1}\frac{\partial f_2}{\partial x_2}\frac{\partial f_2}{\partial x_3}\frac{\partial f_2}{\partial x_4}\mathrm{..........}\frac{\partial f_2}{\partial x_n}\frac{\partial f_2}{\partial \alpha }\\ \mathrm{...........................................................}\\ \mathrm{...........................................................}\\ \frac{\partial f_n}{\partial x_1}\frac{\partial f_n}{\partial x_2}\frac{\partial f_n}{\partial x_3}\frac{\partial f_n}{\partial x_4}\mathrm{..........}\frac{\partial f_n}{\partial x_n}\frac{\partial f_n}{\partial \alpha }\end{array}\right]\text{ (10)}$

α Is the bifurcation parameter. The matrix A can be written in a compact form as

$A=[\frac{\partial f_p}{\partial x_q}.|\frac{\partial f_p}{\partial \alpha }]\text{ (11)}$

The tangent at any point x; ( $z=[z_1,z_2,z_3,z_4,\mathrm{....}{z}_{n+1}]$ ) must satisfy

Az = 0 (12)

The matrix $\left\{\frac{\partial f_p}{\partial x_q}\right\}$ must be singular at both limit and branch points. The n + 1^Th component of the tangent vector Z_n+1 = 0 at a limit point (LP) and for a branch point (BP) the matrix $B=\left[\begin{array}{l}A\\ z^T\end{array}\right]$ must be singular. Any tangent at a point y that is defined by $z=[z_1,z_2,z_3,z_4,\mathrm{....}{z}_{n+1}]$ ) must satisfy

Az = 0 (13)

For a branch point, there must exist two tangents at the singularity. Let the two tangents be z and w. This implies that

Az = 0

AW = 0 (14)

Consider a vector v that is orthogonal to one of the tangents (say z). v can be expressed as a linear combination of z and w ( $v=\alpha z+\beta w$ ). Since Az = Aw = 0; Av = 0 and since z and v are orthogonal,

$z^{T}v=0$ . Hence $Bv=\left[\begin{array}{l}A\\ z^T\end{array}\right]v=0$ which implies that B is singular where $B=\left[\begin{array}{l}A\\ z^T\end{array}\right]$

Let any of the functions fi be separable into 2 functions φ₁, φ₂ as

f_i = φ₁φ₂ (15)

At steady-state $f_i(x,\alpha )=0$ and this will imply that either φ₁ = 0 or φ₂ = 0 or both φ₁. and φ₂ must be 0. This implies that two branches φ₁ = 0 and φ₂ = 0 will meet at a point where both φ₁ and φ₂ are 0.

At this point, the matrix B will be singular as a row in this matrix would be

$[\frac{\partial f_i}{\partial x_k}|\frac{\partial f_i}{\partial \alpha }]\text{ (16)}$

However,

$\begin{array}{l}[\frac{\partial f_i}{\partial x_k}={\varphi }_1(=0)\frac{\partial {\varphi }_2}{\partial x_k}+{\varphi }_2(=0)\frac{\partial {\varphi }_1}{\partial x_k}=0(\forall k=1.,,n)\text{ (17)}\\ \frac{\partial f_i}{\partial \alpha }={\varphi }_1(=0)\frac{\partial {\varphi }_2}{\partial \alpha }+{\varphi }_2(=0)\frac{\partial {\varphi }_1}{\partial \alpha }]=0\end{array}$

This implies that every element in the row $[\frac{\partial f_i}{\partial x_k}|\frac{\partial f_i}{\partial \alpha }]$ would be 0, and hence the matrix B would be singular. The singularity in B implies that there exists a branch point.

The branch point occurred at (pv,nv,dv,cv,kpg) values of (0, 1.55209, 49.999999, 0.524251, 2.120064). Here, the two distinct functions can be obtained from the first ODE in the model.

$\begin{array}{l}\frac{d\left(pv\right)}{dt}=(kpg\left(pv\right)\left(1-\frac{pv}{p\inf }\right)-\frac{kpm\left(sm\right)pv}{\left(\mu m+kmp\left(pv\right)\right)}-kpn\left(fn\right)pv\text{ (18)}\\ fn=\frac{nv}{(1+{(cv/cinf)}^2))}\end{array}$

The two distinct equations are

$\begin{array}{l}pv=0\\ (kpg\left(1-\frac{pv}{p\inf }\right)-\frac{kpm\left(sm\right)}{\left(\mu m+kmp\left(pv\right)\right)}-kpn\left(fn\right)=0\text{ (19)}\\ fn=\frac{nv}{(1+{(cv/cinf)}^2))}\end{array}$

Since, pv = 0, nv = 1.55209. Kpg = 2.120064, kpm = 0.6, kmp = 0.01, sm = 0.005, µm = 0.002, CV = 0.524251, cinf = 0.28, both distinct equations are satisfied, validating the theorem.

The presence of the limit point is beneficial because it allows the MNLMPC calculations to attain the Utopia solution, validating the analysis of Sridhar [22].

The limit and branch points indicate two issues: 1) for a given value of a parameter, there can exist two different profiles, and 2) there exists a singular point that enables one to obtain the best possible outcome even though there may be two conflicting objectives.

Conclusion

Bifurcation analysis and multiobjective nonlinear control (MNLMPC) studies on an inflammatory immune response model. The bifurcation analysis revealed the existence of limit and branch points. The limit and branch points (which cause multiple steady-state solutions from a singular point) are very beneficial because they enable the Multiobjective nonlinear model predictive control calculations to converge to the Utopia point (the best possible solution) in the models. A combination of bifurcation analysis and Multiobjective Nonlinear Model Predictive Control (MNLMPC) for an inflammatory immune response model is the main contribution of this paper.

Data availability statement

All data used is presented in the paper.

Conflict of interest

The author, Dr Lakshmi N Sridhar has no conflict of interest.