A qualitative mathematical model of the immune response under the effect of stress

Although it is common knowledge that stress can hurt our health, scientific results regarding this are fairly recent. The first studies attempting to show this connection between the psychological and the physical state were developed in the first half of the 20th century.^1,2 In the second half of the century, the corresponding field of psychoneuroimmunology (PNI) established the interplay of the central nervous system (CNS), the endocrine system, in other words the hormonal messenger system of the body, and the immune system.^3–6 The first goal of this field is to specify a good definition of stress to be able to work with and thus also a way to categorize it.⁷ The most common way to do this is mainly by the duration, e.g., one may contrast very brief public speaking to extremely long-time care of a spouse.⁸ Furthermore, it is important to quantify stress factors, e.g., how and to what extent we can see sports as positive or negative stressor in the same way as psychological ones.⁹ Furthermore, PNI wants to understand the biochemistry of how behavioral and psychological effects, in particular stress, are reflected in our bodies, i.e., how stress “gets inside the body.”^8,10,11 Although there is still much to learn, we have now a basic understanding of the hormonal reactions of our body to stress and how they affect the immune system, e.g., by specific receptors on cells responsible for the immunity.⁶ In addition, there are several studies regarding immunological and medical consequences of stress in animals^12,13 as well as humans,^14–17 where a key goal is to understand and counteract negative consequences. In this context, the difficulties of extensive clinical trials and the lack of a commonly accepted definition, categorization, or measurements of stressors are most evident. This is why Segerstrom and Miller in 2004 performed a meta-analytical study of numerous results going back 30 years.⁸ Combining all the results and studies we can say that, in general, stress has a negative effect on the immune system by either decreasing the number of, e.g., killer cells or simply disturbing the equilibrium between the different components. Nevertheless a brief stress, like public speaking, can also have a positive effect and even longer short-term stressors like examination periods for students can have a beneficial effect during but cause decreased immunity afterwards.^18,19

Although quite a number of studies in PNI exist, there is not yet a simplest qualitative nonlinear dynamics model available for the interplay between the immune system and stress. There are many other scientific disciplines pertaining to human health and behavior, where designing and studying low-dimensional benchmark models has already helped tremendously. This applies to smaller scales in the context of single-neuron models,²⁰ chronobiology,²¹ genomics,²² or cancer research,²³ as well as to larger scales such as collective motion,²⁴ crime hotspots,²⁵ or epidemiology;²⁶ we also refer to Ref. 27 for a broader view on the contributions mathematical models can make in various contexts to health. One key advantage of a reduced modeling approach via the essential nonlinearities is that key effects are isolated and detailed parameter studies become feasible. In particular, the toolbox of nonlinear dynamics for low-dimensional models is extremely powerful^28–30 and has led to a clear identification of new mechanisms in complex systems, e.g., we understand different excitability patterns in neurons much better via low-dimensional dynamics models.^20,31,32 Yet, it is not clear what the elementary nonlinear mechanisms in the context of PNI, and, in particular, of the interplay between stress and the immune system, are.

Having some fundamental understanding of how our immune system reacts to stress, we want to find a simple mathematical model of differential equations, which is able to represent the basic qualitative interaction between an infection and our immune system taking into account the effect of stress.

Since we are only interested in a qualitative analysis of the system, we are going to work with normalized dynamical variables $x$⁠, representing the level of activity of the immune system, and $y$⁠, as the level of infection or sickness, with values ranging between 0 and 1. The stress is included by a parameter $s$ also normalized between 0 and 1, where 0 corresponds to no stress at all and 1 is the maximal level of stress limited by the assumption that a person cannot feel infinitely strong stress.

To specify the dynamics of the infection, we are going to build upon principles of mathematical biology,^33,34 where in multiple sub-disciplines such as ecology, neuroscience, or biomechanics, one aims to generic polynomial nonlinearities to represent basic effects. Based upon population dynamics, we consider logistic growth and Allee effect for the immune system intrinsic dynamics. In this regard, our first assumption is that, as long as our immune system has some activity, the infection should have two stable states, the “sick”-state and the “healthy”-state. Furthermore, since we are constantly coming in contact with different viruses and bacteria, we are interested in having a slightly raised baseline instead of the “healthy”-state being $y=0$⁠. Our second assumption is that, if there is no immune system or if it is too weak, then the infection should spread up to its maximal capacity. Combining everything, we arrive at the equation

where we used a square root instead of the standard multiplicative term $(1−y)$ from population models. Finally, we choose the parameter values $r=2$⁠, $y0=y1=0.1$⁠, and $q=0.95$⁠.

In contrast to the previous derivation of the dynamics for $y$⁠, our approach for finding an appropriate equation for the immune system is based purely on the geometry necessary to capture the effects we aim to model. We start with the case without stress. It is clear that without stress, we should have a unique stable equilibrium with $0≪x<1$ and $0<y≪1$ corresponding to a normal healthy state. Since too much stress generally induces an infection, this fixed point should only exist for $s<s0$ with some threshold level of stress $0<s0<1$⁠.

Furthermore, we also want to reflect the fact that, while a moderate stress can make slightly more resistant to infections, at the same time, its sudden drop to normal levels can be detrimental increasing the likelihood of infections, e.g., this naturally occurs for short-term stressors. To capture this effect, we want to have a second unstable equilibrium with $0<y≪1$⁠. As $s$ increases, the $x$-coordinates of both fixed points should slowly start decreasing and approaching each other such that for some $0<s1<s0$ the stable equilibrium crosses the original position of the unstable one. If we would now drop the stress back to $s=0$⁠, it would in fact trigger the expected infection before converging to the healthy state.

Finally, it is reasonable to assume that an extreme (prolonged) stress can prevent the body from recovering from an infection, giving a “burn-out” state, such that for $s$ close to 1 we want to have a fixed point where $0≪y<1$⁠.

Again combining all the assumptions above, one simple polynomial nonlinear differential equation capturing these effects is given by

presenting an S-shaped nullcline with $x0=0.3$ and the parameter $k=1.3−0.3(1−s)4$ to control the shift to the right of the equilibria and $x1=0.8−0.2s2$ responsible for the disappearance of the stable healthy state. Of course, the actual numerical values in (2.2) are not the key aspect but the geometry of the dynamics and the qualitative description of the different aspects of immune reaction, infection level, and stress.

In summary, the complete system we arrive at

with the parameter values

integrating the effect of stress into the equations. The additional parameter $α$ controls the scale of change of the infection $y$ with respect to the immune system $x$⁠. Figure 1 shows the corresponding phase portraits and time series for multiple values of $s$⁠. And indeed, our model presents different dynamic regimes upon different stress levels, which we can observe in the phase portraits and time series shown in Fig. 1.

To analyze in more detail the change between the different dynamical regimes, we take a look at the corresponding bifurcation diagram²⁹ in Fig. 2.

For $s=0$⁠, we have three equilibria: an unstable spiral, a saddle, and a stable node. As $s$ increases, the stable node and the saddle collide and disappear in a saddle-node in cycle (SNIC) bifurcation for $s≈0.48$⁠. At $s≈0.95$ we find a second SNIC, where another pair of a stable node and a saddle appear. In addition to the equilibria mentioned above, in the parameter regime between the two bifurcation points we find that all orbits except the unstable equilibrium have to converge to a limit cycle as $t→∞$⁠. This can be easily shown using the Poincaré–Bendixson theorem on the unit square. Note that SNIC bifurcations are of co-dimension one, i.e., they are generic/typical in planar systems with one parameter presenting a mechanism to obtain oscillations. They have been observed already in many models ranging from neuroscience^35,36 to mechanics.³⁷ They are accompanied by a heteroclinic loop for parameter values $s<0.48$ (or $s>0.95$⁠) since the unstable manifolds of the saddle converge to the stable node. As $s$ approaches the bifurcation, the two equilibria collide, one of the heteroclinic connections disappears and the second one becomes a homoclinic orbit. By perturbing $s$ further, this homoclinic gives rise to a family of stable limit cycles for $0.48<s<0.95$⁠.

In this work, we have initialized the qualitative mathematical model development for the interaction between immune system levels, infection dynamics, and stress level. Based upon established principles of theoretical biology, we developed a simple, yet powerful, planar dynamical system to analyze the influence of stress as an external parameter. Using only elementary nonlinearities, out model can represent (I) a stable healthy equilibrium for low stress, (II) the effect of sudden stress level drop from moderate to low levels inducing a single infection, (III) periodic outbreaks of infections under high stress, (IV) the transition to a burn-out state at very high stress. The key transition mechanism we have identified is a saddle-node in cycle (SNIC) bifurcation. This bifurcation transition actually makes the model predictive, e.g., it predicts that if stress levels are close to a burn-out stage but still below, then the oscillatory sick periods become longer. Although this looks like a natural prediction, it is one that we did not anticipate to be able to extract from our system at all during the modeling process.

Having established the ability of conceptual mathematical models to contribute to the understanding of stress, we believe that further model development and connecting qualitative mathematical models to more detailed biophysical principles as well as to clinical trials, could forge an effective path mitigating, and probably even more importantly, predicting in advance, the positive and negative effects of stressors.

C.K. conceived the original idea. Both authors contributed into developing the idea further to the final form of the publication. M.E.G.H. developed the theory and performed the computations under the supervision and guidance of C.K. M.E.G.H. drafted the article. Both authors revised the manuscript critically and contributed to its final form.

M.E.G.H. work was partly supported by a Hans-Fischer Senior Fellowship of the Technical University of Munich—Institute for Advanced Study (TUM-IAS) and the TopMath Program of the Elite Network of Bavaria. C.K. acknowledges partial support of the VolkswagenStiftung via a Lichtenberg Professorship. We also thank an anonymous referee, whose comments have helped to improve the presentation of our results.

All figures were generated using MATLAB R2020b. Figure 2 was generated with the MATLAB continuation package MATCONT 7.1. The code generating Fig. 1 as well as the model used in MATCONT to generate Fig. 2 has been made available in https://doi.org/10.5281/zenodo.4584715, Ref. 39.