Non-minimum-phase dynamics in an inflammation model

Inflammation is the innate immune system’s response to tissue damage caused by trauma or infection. Thinking about it as the quickly receding rash after an insect bite detracts from the major antagonistic role it plays in medicine. As Baldur Tumi Baldursson of the National University Hospital of Iceland put it, `I tell my students, your work is inflammation. Practically all of internal medicine is just fighting inflammation.’

Inflammation is a complex process involving different cell types, mediator molecules and changes in the permeability of capillaries and affected tissue. This reaction has to be ramped up quickly to defend our body, but it must also be kept under strict control to prevent immune cells causing tissue damage themselves. On occasion, things go wrong and a patient is left with chronic, abnormal inflammation: inflammatory bowel disease, coeliac disease or rheumatoid arthritis are a few debilitating examples.

A mathematical model of inflammation control

These notes extend the findings of three researchers from Britain—Joanne Dunster (Reading), Helen Byrne (Oxford) and John King (Nottingham)—using control theoretical insights, and perhaps contribute a new aspect to our understanding of the problem. Their 2014 paper¹ drew attention to a shift from understanding inflammation resolution as a passive process to an active, anti-inflammatory mechanism. The authors argued that the interaction of different cell types (neutrophils and macrophages) is crucial to this process. The paper expounded on three mathematical models of increasing level of detail.

The basic script of the modelled process goes like this. The scene is some generic sterile soft tissue, which suffered mechanical stress. The players are two different white blood cell types, neutrophils and macrophages, and a generic pro-inflammatory mediator, which encompasses molecules like THF- $\alpha$ and IL-8. Neutrophils are subdivided into two subtypes: active and apoptotic.

The initiating physical damage causes a local increase in the amount of the pro-inflammatory mediator. Its higher concentration attracts neutrophils to the area. However, the active neutrophils die at a certain rate, and transform into apoptotic neutrophils. As they break down, they spill toxic load (e.g. reactive oxygen species), damaging tissue and actually increasing the level of the pro-inflammatory mediator. Macrophages accompany neutrophils in responding to the original signal, and once there, they remove apoptotic neutrophils. While neutrophils on their own prolong inflammation by the positive feedback loop, macrophages present a negative feedback loop.

Let me show the first of the paper’s three ordinary differential equation models. The concentrations of active and apoptotic neutrophil populations are denoted by $n$ and $a$ , respectively, that of macrophages by $m$ , and the concentration of pro-inflammatory mediator molecules by $c$ . The mechanical stress is represented by the input function $f(t)$ . With the necessary rate constants, the model reads

$\dot{n}(t)=\chi_n c(t)-\nu n(t)$

$\dot{a}(t)=\nu n(t) -\gamma_a a(t) - \phi m(t) a(t)$

$\dot{m}(t)=\chi_m c(t) - \gamma_m m(t)$

$\dot{c}(t)=\alpha f(t) +k_a \gamma_a \left(\frac{a(t)^2}{\beta_a^2+a(t)^2}\right)-\gamma_c c(t)$

The initial condition is $n(0) = a(0) = m(0) = c(0) = 0$ . The outcome of the inflammation is either resolution, when the four concentrations return to their initial values of zero; or chronic, self-perpetuating inflammation, when the system stabilises in a positive state, meaning that white blood cells remain active in the inflamed tissue. The authors expect any inflammation model to be able to reproduce both outcomes, this kind of bistability. This model exhibits both behaviours in numerical simulations (the outcome can be influenced by how prolonged the stimulation is), and a linear stability analysis explains the observed behaviour.

Moreover, a bifurcation diagram indicates that an increase in the rate of macrophage phagocytosis (macrophages clearing apoptotic neutrophils) increases the likelihood of healthy resolution. Additionally, independent research had suggested that the rate of neutrophil apoptosis $\nu$ also influences the outcome. However, this first model does not exhibit this trait, the existence of steady states is not influenced by $\nu$ . This motivated Dunster et al. to develop another model with an additional positive feedback loop from activated neutrophils via pro-inflammatory mediators.

The second model, where active neutrophils release pro-inflammatory mediator

The additional assumption here is that neutrophils fight pathogens (even mistakenly, when pathogens are not present) by the release of toxins, thereby increasing levels of the pro-inflammatory mediator. After a rescaling for clarity (with the new variables and constants denoted by the original symbols), Model 2 takes the following non-dimensionalised form:

$\dot{n}(t)=c(t)-\nu n(t)$

$\dot{a}(t)=\nu n(t) -\gamma_a a(t) - \phi m(t) a(t)$

$\dot{m}(t)=c(t) - \gamma_m m(t)$

$\dot{c}(t)=\alpha f(t) +k_n \left(\frac{n(t)^2}{\beta_n^2+n(t)^2}\right) + \gamma_a \left(\frac{a(t)^2}{\beta_a^2+a(t)^2}\right)- c(t)$

This recovers Model 1 when $k_n=0$ . (Here, due to the rescaling, e.g. $k_n=k_n^{\mathrm{original}}k_a/\gamma_c.)$

For physiologically relevant parameter values $k_n \ll k_a \gamma_a$ , $\nu$ can now influence model behaviour: if neutrophil apoptosis rate is reduced, then the monostable system becomes bistable. (To best see what is going on, compare the first panels of Figures 8 and 11 in the paper.)

The conclusion from Model 2 is that for the resolution of inflammation, not only a high enough rate of phagocytosis is needed, but in my opinion counterintuitively, this should be accompanied by a rate of neutrophil apoptosis that is not too low to rule out bistability. I find this counterintuitive because neutrophil apoptosis in the short term causes a spike in the generic pro-inflammatory mediator. On the other hand, in the long run the removal of neutrophils contributes to inflammation resolution.

Some control theoretical insights in systems biology

From the viewpoint of control engineering, inflammation is an interesting phenomenon because it must exhibit rapid amplification, but it should not overshoot, requiring effective feedback control. It is natural to anticipate that this discipline can tell us something new about inflammation, and I present evidence in the following.

To me, the counterintuitive behaviour was reminiscent of non-minimum-phase dynamics. A standard example is driving a car in reverse: when we turn the steering wheel, the car turns in the direction we want to go, but not before swinging in the opposite direction². I can recommend the whole review² by my former postdoc supervisor Mustafa Khammash about insights in biology from robust control.

The paper 3, which preceded the review, wrote:

Control engineers make considerable efforts to avoid non-minimum phase systems, as they are generally subject to undesirable tradeoffs […] [T]he attenuation of some disturbances must come at the expense of the amplification of others. […] [I]ncreasing feedback gains, which might otherwise help with the rejection of disturbances, can lead to increasing oscillation and, ultimately, instability.

The authors Buzi, Lander and Khammash went on to argue that a simple feedback mechanism with non-minimum-phase behaviour is not robust to perturbations and is unlikely to be prevalent in stem cell differentiation. Instead, they proposed and predicted a modified differentiation feedback architecture with lineage branching, which exhibits minimum phase behaviour. They also presented possible biological examples consistent with the latter model.

It is worth highlighting here Gunter Stein’s Bode Lecture from 1989, published in 2003⁴, as a classic on robust control and on its fundamental limitations. It is highly recommended reading.

The claim I make is that Model 2 possesses non-minimum-phase dynamics in a biologically relevant parameter regime, as I show with calculations performed with the computer algebra system Mathematica by Wolfram Research, Inc. Model 1 is a boundary case which just misses that. I do not discuss Model 3 of Dunster et al. because the calculations were too difficult even with Mathematica to reach a conclusion.

Non-minimum-phase dynamics in continuous-time linear systems

Consider a rather generic linear time-invariant (LTI) system describing the evolution of state vector $x,$ the influence of a potentially multidimensional input function $u$ , and the output function $y$ (essentially a readout), using fixed matrices of appropriate sizes $A, B, C$ :

$\dot{x}(t)=Ax(t)+Bu(t),\qquad x(0)=x_0$

$y(t)=Cx(t)$

The solution of this system is explicitly

$y(t)=Cx(t)=C\left(\mathrm{e}^{At}x(0) + \int_0^t \mathrm{e}^{A(t-\tau)}Bu(\tau)\,\mathrm{d}\tau\right)$ .

When I studied complex analysis, one of the take-home messages ingrained in us was that whenever you see a convolution, like here from input to output, you must take its Fourier transform because that transforms it into a more tractable product. Control engineering conventionally takes the similar Laplace transform, see e.g. Section 2.6 of Reference 5.

For the special case of $x(0)=0,$ the Laplace transform of the equation is the following relation between the Laplace transforms of input and output, $\hat{u}$ and $\hat{y}:$

$\hat{y}(s)=\hat{G}(s)\hat{u}(s)$ .

$\hat{G}$ is called the transfer function of the system. In this case it has the form

$\hat{G}(s)=C(sI-A)^{-1}B,$

where $I$ is an identity matrix of the same size as the square matrix $A$ , so $sI$ is a diagonal matrix with $s$ in every diagonal entry. $\hat{G}$ is a rational function (it is the ratio of two polynomials).

I won’t go into detail why, but the system is called minimum phase if both the roots (also called zeros) and poles of $\hat{G}$ are inside the left half of the complex plane (their real parts are negative). Conversely, assuming that the system is stable and its poles are in the left half-plane, it is non-minimum phase if it has a zero in the right half-plane. These systems are causal and stable but their inverses, while being causal, are unstable.

As explained in Refs. 3 (Box 2 in Methods) and 4, Bode’s integral formulae express conservation laws about how good controllers can be at compensating disturbances. Even in an ideal case, some limit on performance cannot be surpassed. In non-minimum-phase systems, there is always some frequency where disturbances are amplified instead of rejected. This is why you can reverse with your car at low speed, but at higher speeds you’d be confused and would lose control.

Calculations in the inflammation models

In the following, I show how to compute the roots of the transfer function in Mathematica, first for the more transparent Model 1. We start by linearising the system around a steady state, so we copy the Jacobi matrix from (5)¹. To follow the notation we have used so far, I denote the system matrix by $A$ instead of $J:$

A = {{-\[Nu], 0, 0, 1}, {\[Nu], -\[Gamma]a*P, -\[Phi]*a, 0}, {0, 
   0, -\[Gamma]m, 1}, {0, Q, 0, -1}}
MatrixForm[A]

$P$ and $Q$ are shorthands for longer terms. The transfer matrix requires $(sI-A)^{-1},$

MatrixForm[s*IdentityMatrix[4] - A]
sImA = Inverse[s*IdentityMatrix[4] - A]

whose entries are filled with rational functions. Now on to the choice of $B$ and $C$ .

Because we are interested in the effect of changing the rate of neutrophil apoptosis $\nu$ , I turn this into an input by selecting $B=[-\nu_2 \ \nu_2 \ 0 \ 0]^{\mathrm{T}}$ for some $\nu_2 >0$ . (Actually, $\nu_2 =1$ would also do.) By increasing input $u$ , we can increase the rate of active neutrophils $n$ transitioning into the apoptotic state $a.$

What should be the measured quantity, how to choose $C?$ Ideally, I would choose the active neutrophil concentration $n$ and the corresponding $C=[1\ 0\ 0\ 0].$ However, this choice leads to a calculation that is so convoluted even with Mathematica that I could not extract any useful knowledge. Instead, I picked the concentration of pro-inflammatory mediator $c$ and with that, $C=[0\ 0\ 0\ 1]$ . This is a good proxy of inflammation activity.

C0 = {{0, 0, 0, 1}} 
B1 = Transpose[{{-\[Nu]2, \[Nu]2, 0, 0}}]
L1 = Together[C0.sImA.B1]
zeros = Solve[L1 == 0, s]

The result is $\{\{s\to 0\}, \{s\to -\gamma_m \} \}.$ One root is negative, the other is on the boundary between the left and right half-planes.

The solution for $C=[0\ 1\ 0\ 0]$ is $\{\{s\to -1\}, \{s\to 0\}, \{s\to -\gamma_m \} \},$ and for $C=[0\ 0\ 1\ 0]$ , it is $\{\{s\to 0\}\}.$ There are no roots with positive real part.

In Model 2, the calculation is more challenging. The linearisation (9b) is entered by

A = {{-\[Nu], 0, 0, 1}, {\[Nu], -\[Gamma]a*P, -\[Phi]*a, 0}, {0, 
   0, -\[Gamma]m, 1}, {R, Q, 0, -1}}
MatrixForm[A]
sImA = Inverse[s*IdentityMatrix[4] - A]
B1 = Transpose[{{-\[Nu]2, \[Nu]2, 0, 0}}]

Once again, it is very hard to make use of $C=[1\ 0\ 0\ 0],$ and I turn to pro-inflammatory mediator concentration $c$ and $C=[0\ 0\ 0\ 1]:$

C0 = {{0, 0, 0, 1}}
L1 = Together[C0.sImA.B1]

I substitute the values of $P, Q$ and $R:$

numL1 = Collect[Numerator[L1], 
   s] /. {P -> 1 + a (\[Gamma]m/\[Phi] - a)^(-1), 
   Q -> 2 \[Gamma]a \[Beta]a^2 a (\[Beta]a^2 + a^2)^(-2), 
   R -> (2 kn \[Beta]n^2 n) (\[Beta]n^2 + n^2)^(-2)}

to solve the equation for the zeros of the transfer function

zeros = Solve[numL1 == 0, s]

One root is $\{s\to -\gamma_m \},$ but the other is not pretty: $\{s\to (k_n n (a^2 + \beta_a^2)^2 \beta_n^2 \gamma_a \gamma_m)/((a^4 k_n n \beta_n^2 + 2 a^2 k_n n \beta_a^2 \beta_n^2 + k_n n \beta_a^4 \beta_n^2 - a n^4 \beta_a^2 \gamma_a - 2 a n^2 \beta_a^2 \beta_n^2 \gamma_a - a \beta_a^2 \beta_n^4 \gamma_a) (-\gamma_m + a \phi))\}.$

This is a fraction with enumerator $S=k_n n (a^2 + \beta_a^2)^2 \beta_n^2 \gamma_a \gamma_m >0$ and denominator $T=(T_1-a T_2)T_3,$ with

$T_1=a^4 k_n n \beta_n^2 + 2 a^2 k_n n \beta_a^2 \beta_n^2 + k_n n \beta_a^4 \beta_n^2 = k_n n \beta_n^2 (a^4 + 2 a^2 \beta_a^2 + \beta_a^4)$

$a T_2=a (n^4 \beta_a^2 \gamma_a + 2 n^2 \beta_a^2 \beta_n^2 \gamma_a + \beta_a^2 \beta_n^4 \gamma_a) = a \beta_a^2 \gamma_a (n^4 + 2 n^2 \beta_n^2 + \beta_n^4)$

$T_3=-\gamma_m + a \phi$

It helped me to use the bit of code

Collect[a^4 kn n \[Beta]n^2 + 2 a^2 kn n \[Beta]a^2 \[Beta]n^2 + 
  kn n \[Beta]a^4 \[Beta]n^2 - a n^4 \[Beta]a^2 \[Gamma]a - 
  2 a n^2 \[Beta]a^2 \[Beta]n^2 \[Gamma]a - 
  a \[Beta]a^2 \[Beta]n^4 \[Gamma]a, a]
Collect[a^4 kn n \[Beta]n^2 + 2 a^2 kn n \[Beta]a^2 \[Beta]n^2 + 
  kn n \[Beta]a^4 \[Beta]n^2, kn n \[Beta]n^2]
Collect[-n^4 \[Beta]a^2 \[Gamma]a - 
  2 n^2 \[Beta]a^2 \[Beta]n^2 \[Gamma]a - \[Beta]a^2 \[Beta]n^4 \
\[Gamma]a, \[Gamma]a \[Beta]a^2 ]

to discover this. $n=\bar{n}$ and $a=\bar{a}$ are the positive steady state values of the respective variables, everything else is a fixed parameter.

According to Table 4¹, $\gamma_m /\phi \approx 10$ . From e.g. Figure 6 or 13, $\bar{a}$ < $10$ is perfectly plausible, and consequently, so is $T_3$ < $0$ . (I refer to figures instead of an explicit form because $\bar{a}$ is only known as the solution of a quartic equation.)

Furthermore, from Figure 6, $\bar{a}$ < $\bar{n}$ < $100$ is to be expected, and from Table 4, $\beta_a\approx \beta_n$ and $\gamma_a / k_n \approx 100.$ So if $\bar{n}$ < $100 \bar{a}$ , then $k_n \bar{n} \beta_n^2$ < $\bar{a}\beta_a^2 \gamma_a$ , and $T_1$ < $\bar{a}T_2$ . Then $T>0$ , and with $S>0$ , we have got non-minimum-phase dynamics. $S>0$ is dependent on $k_n>0$ , which distinguishes Model 2 from Model 1.

This was the case of $C=[0\ 0\ 0\ 1]$ in Model 2. Analogously to $C=[1\ 0\ 0\ 0],$ it is hard to say anything about $C=[0\ 1\ 0\ 0].$ Finally, with $C=[0\ 0\ 1\ 0],$ there is only one root, which is positive under the previous assumptions (and exactly those assumptions are needed).

Conclusion

The calculations show that something is strange either with the above Model 2 of inflammation by Dunster et al., or with inflammation itself. It exhibits non-minimum-phase dynamics around the positive steady state, which is in general an unwanted property. Whether this quirk has any significant role in the development of chronic inflammation leaves me curious. Is there a frequency range where inflammation is more sensitive to disturbances, and do such disturbances actually exist? Are they relevant in any pathologies? If some of you (perhaps with more experience in control engineering) has something to add, please let me know.

(Commenting requires an email address, but it will not be publicly displayed.)

References

¹Joanne L. Dunster, Helen M. Byrne, and John R. King. The resolution of inflammation: a mathematical model of neutrophil and macrophage interactions. Bulletin of Mathematical Biology, 76(8):1953–1980, 2014. doi: 10.1007/s11538-014-9987-x ($).

²Mustafa Khammash. An engineering viewpoint on biological robustness. BMC Biology, 14:22, 2016. doi: 10.1186/s12915-016-0241-x.

³Gentian Buzi, Arthur D. Lander, and Mustafa Khammash. Cell lineage branching as a strategy for proliferative control. BMC Biology, 13:13, 2015. doi: 10.1186/s12915-015-0122-8.

⁴Gunter Stein. Respect the unstable. IEEE Control Systems Magazine 23:12–25, 2003. doi: 10.1109/MCS.2003.1213600 ($). Free download.

⁵Geir E. Dullerud, Fernando G. Paganini. A course in robust control theory: a convex approach. Vol. 36. Springer-Verlag New York, 2000. doi: 10.1007/978-1-4757-3290-0.

A mathematical model of inflammation control

The second model, where active neutrophils release pro-inflammatory mediator

Some control theoretical insights in systems biology

Non-minimum-phase dynamics in continuous-time linear systems

Calculations in the inflammation models

Conclusion

References

Share this:

Like this:

Related

Leave a Reply Cancel reply