The innate immune system is the phylogenetically oldest part of the immune system. It is composed of several humoral and cellular components and has evolved in parallel with the development of multicellular life within a period of 2.4 billion years, which corresponds to 75% of the total evolutionary time (Storch et al., 2013; Delves et al., 2016). Pathogens first come into contact with the innate immune system, which alone can render harmless over 99% of all potential threats. In addition to destroying bacteria, it is also capable of very efficiently attacking and destroying endogenous cells infected by viruses, thus stopping viral replication.

The function of the innate immune system is maintained constantly throughout the lifetime with almost the same level of response by spatially mobile cells throughout the organism. The communication in order to identify an infection and its location, the initiation of an acute phase response and the simultaneous control of inflammatory response including its extend is provided by cytokines and other mediators. These cytokines and mediators are distributed in the organism through the blood stream. If they meet cells with corresponding receptors, these can respond to the cytokine signals. Cytokine sources that do not originate from the immune system are considered as perturbations and can usually alter the balance of the inflammatory response in a proinflammatory direction. Unregulated sources of cytokines not originating from the immune system are adipose tissue, acute and chronic inflammation, and concomitant diseases. Lifestyle factors such as physical inactivity or smoking also influence the cytokine dynamics. Cytokine polymorphisms are responsible for the resilience of the innate immune system upon perturbations and thus for the high individual inflammatory host response (Rosendal et al., 2004; Egger, 2005; Tisoncik et al., 2012; Schulte et al., 2013; Hotchkiss et al., 2016; Xia et al., 2016; Elisia et al., 2017; Thomas, 2020).

The pathophysiological situation is complicated by the fact that many pathogens (bacteria, fungi, viruses, endogenous material) can trigger an inflammatory response. Moreover, the innate immune system consists of many interacting components, there are many inflammation triggering pathways, and the signaling pathways and cytokines have a high redundancy and additionally a pronounced pleiotropy. Inflammation is usually localized, encapsulated and healed by destruction and phagocytosis of destroyed cells. Via the cytokines IL-1, IL-6 and TNF-α released locally in the inflammation focus by macrophages, lymphocytes, fibroblasts and endothelial cells, a multistage defense process is started. Further, cytokines stimulate the anterior pituitary to synthesize cortisol in the adrenal cortex. Cortisol stimulates hepatocytes to synthesize cytokine receptors, which can then receive the cytokine signals and produce acute phase proteins (APP). Besides, temperature elevation occurs due to central nervous system stimulation and leukopoiesis is enhanced in the bone marrow. Acute phase proteins comprise a variety of proteins that restrict the inflammatory process. The functionally distinct proteins are produced and released step by step according to the course of the inflammatory response and are controlled by feedback mechanisms. At the center of inflammation, inhibition of the inflammatory response does not occur due to the stoichiometric ratio of acute phase proteins to proinflammatory cytokines. In the bloodstream, the ratio is reversed, acute phase proteins are clearly dominating, and they neutralize proinflammatory cytokines and can thus prevent the start of systemic inflammation.

If a local focus of inflammation cannot be adequately localized by the acute phase response and if its supplies in the blood are depleted by consumption, the proinflammatory cytokines, mediators, and immune cells have the potency to damage or destroy the function of organs far from the focus of inflammation. Reactive oxygen species (ROS) and other proinflammatory cytokines are released via cytokine-induced activation of polymorphonuclear leukocytes (PMNs) and macrophages in the bloodstream and their interaction with endothelium. This process creates the initial condition for the Systemic Inflammatory Response Syndrome (SIRS) (Egger, 2005).

The mechanism of damage in the systemic inflammatory response syndrome and sepsis is, on the one hand, the lack of oxygen availability to the parenchyma due to disruption of the microcirculation by intravascular coagulation triggered by inflammation. In addition or alternatively, cytokines may induce a shutdown of mitochondrial cellular respiration. Cellular oxygen utilization now occurs only via aerobic glycolysis. Which process dominates in which phase of the disease, in which organ, or in which patient is not yet known. Cytokines have cardiotoxic and central nervous system toxic effects.

The long-term outcome after surviving sepsis or septic shock is poor. Late effects include myopathic, neuropathic, and cognitive changes, worsening of pre-existing conditions, and increased mortality. Long-term survival is reduced regardless of pre-existing conditions, and 74% of patients are deceased two years after illness. Causes include re-infection with sepsis, cardiovascular disease and tumors. The increased vulnerability after survived sepsis is attributed to the dysregulated inflammation during the acute phase of the disease with the tissue damage that occurred in the acute process and the continuing inflammation (Prescott et al., 2016; Mostel et al., 2019; Brunkhorst et al., 2020; Schmidt et al., 2020).

A schematic illustration of organic tissue consisting of parenchymal cells and immune cells is shown in Figure 1. Panel A depicts the initial configuration of a tissue element. The tissue element consists of the epithelial parenchyma, the basal membrane and the stroma. The parenchyma is the organ-specific functional layer. The basal membrane separates the parenchyma and stroma and is made of collagen of type IV which is a network-forming collagen underlying epithelial and endothelial cells. In the stroma, blood supply, lymphatic drainage and immune response occur. The stroma consists of an extracellular matrix and embedded cells that do not form a solid association. The extracellular matrix is structurally composed of collagen, glycoproteins, proteoglycans and water. Cells in the stroma are resident fibroblasts and fat cells, and mobile cells (macrophages, mast cells, granulocytes, and plasma cells). Panel B shows the functional interactions in the two-layer network model of the parenchymal layer and the immune layer.

Figure 1 shows the functional structure of the tissue element, in particular the reactants interacting during sepsis. With the blood supply via capillaries, pro- and inflammation-inhibiting molecules are delivered to the stroma of each organ. They originate from the primary focus of infection (pathogen-associated molecular patterns, damage-associated molecular patterns, cytokines), from the liver (acute phase proteins) and from the innate immune system (macrophages, polymorphonuclear leukocytes). The concentration of all reactants changes as the inflammatory response progresses. They initially interact with the endothelium of the capillaries. With the influx of pro- and inflammation-inhibiting reactants, the overall system (Figure 1A) tries to maintain a local inflammation-inhibiting equilibrium. Blood flow and oxygen supply, especially to the parenchyma, must be ensured.

An ongoing blood flow and oxygen supply is achieved by the individual and locally adapted information processing of all cells of the innate immune system (macrophages, polymorphonuclear leukocytes), of the stroma (endothelial cells, fibroblasts), the specific activation of platelets, the pleiotropy of cytokines, i.e., their concentration- and pattern-dependent reaction patterns, and the acute phase proteins produced and released in the liver via cytokines in a time-delayed manner. All cells involved are potential cytokine sources.

The pathophysiological positive response pattern is the maintenance of the inflammation-inhibitory balance. The pathological situation is the initiation of disseminated intravascular coagulation, interruption of blood flow, oxygen diffusion pathways prolonged by fluid influx into the stroma, and breakdown of parenchymal oxygen supply. In parallel and in addition, cytokines interact with the parenchyma and reduce parenchymal function via impairment of mitochondrial cellular respiration. This process may develop an autocatalytic characteristic with the involvement of reactive oxygen species, ending in organ failure.

The unified disease model is centered around the nonspecific immune system, which includes disease-specific initial conditions and infection-driven cytokine dysregulation. For the analysis of an emergening sepsis, we consider a volume element of tissue consisting of parenchyma, basal membrane and stroma, see Figure 1A. In (Sawicki et al., 2022), we have introduced a functional model to describe the dynamic interaction of parenchyma (organ tissue) and stroma (immune layer). The network layer of parenchymal cells (superscript 1) are represented by N phase oscillators ϕ i 1 , i = 1, … , N and the network layer of immune cells (superscript 2) are presented by N phase oscillators ϕ i 2 . The coupling weights in the parenchymal layer are considered to be partly fixed and partly adaptive while in the immune layer the coupling weights are completely adaptive. We model the communication through cytokines which mediate the interaction between the parenchymal cells by the coupling weights κ i j 1 , and those between the immune cells by coupling weights κ i j 2 . Note that ϕ i 2 and κ i j 2 represent the collective dynamics of all dynamical units of the stroma, see Figure 1B. Hence, this set of variables can be regarded as collective dynamical variables used in our functional modeling approach. The use of phase oscillators for the functional modeling of the interacting parenchymal cells and immune cells is motivated by the fact that phase oscillator networks are a paradigmatic model for collective coherent and incoherent dynamics. The healthy state is assumed to be characterized by regular periodic, fully synchronized dynamics of the phase oscillators. Healthy and pathological cells differ by their metabolic activity, i.e., pathological cells shut down their mitochondrial cellular respiration and switch to aerobic glycolysis. Therefore they are less energy-efficient and thus have a modified cellular metabolism and reduced function, which is reflected in our phase oscillator model by a different frequency, and the system splits into multifrequency clusters.

We consider a general multiplex network with two layers each consisting of N identical adaptively coupled phase oscillators: ϕ ̇ i 1 = ω 1 − 1 N ∑ j = 1 N a i j 1 + κ i j 1 sin ϕ i 1 − ϕ j 1 + α 11 − σ ⁡ sin ϕ i 1 − ϕ i 2 + α 12 , κ ̇ i j 1 = − ϵ 1 κ i j 1 + sin ϕ i 1 − ϕ j 1 − β , (1) ϕ ̇ i 2 = ω 2 − 1 N ∑ j = 1 N κ i j 2 ⁡ sin ϕ i 2 − ϕ j 2 + α 22 − σ ⁡ sin ϕ i 2 − ϕ i 1 + α 21 , κ ̇ i j 2 = − ϵ 2 κ i j 2 + sin ϕ i 2 − ϕ j 2 − β , (2) where ϕ i μ ∈ [ 0,2 π ) represents the phase of the i − th oscillator (i = 1, … , N) in the μ − th layer (μ = 1, 2), ω ^μ are the natural oscillator frequencies of the oscillators in the μ − th layer. The interaction between the oscillators within each layer is determined by the intralayer connectivity weights a i j 1 ∈ [ 0,1 ] (fixed interaction within an organ) and κ i j μ ∈ [ − 1,1 ] (adaptive interaction mediated by cytokines). We assume that the parenchymal layer has both fixed and adaptive couplings, while the immune layer has only adaptive coupling. Further the interactions within the layer depend on the phase lag parameters α ¹¹ and α ²².

In this work our focus is on the interaction between the two layers and their synchronization. In particular, we analyze the onset of desynchronization in the parenchymal layer induced by an activated immune layer. The interaction of the layers is controlled by two main parameters, the interlayer coupling weight σ and the interlayer phase lag parameters α ¹² and α ²¹. Between the layers the interlayer coupling weights σ ≥ 0 are fixed and symmetric for both directions of interaction. The phase lags can be considered to model interaction time delays (Sakaguchi and Kuramoto, 1986; Madadi Asl et al., 2018).

The adaptation rates 0 < ϵ ^μ ≪ 1 separate the time scales of the slow dynamics of the coupling weights and the fast dynamics of the oscillatory system. The adaptation rate of the parenchymal layer ϵ ¹ is assumed to be slow compared to the adaptation rate of the immune layer ϵ ², i.e., ϵ ¹ ≪ϵ ² to account for the faster reaction of the immune cells, see also (Sawicki et al., 2022). Thus we have two classes of adaptive coupling weights modeling two different cytokine mechanisms on two different timescales. Consequently by choosing two significantly different values for ϵ ¹ and ϵ ², a system with multiple times scale dynamics is obtained, i.e., “slow-fast-faster” dynamics (ϵ ¹ ≪ϵ ² ≪ 1) (Desroches et al., 2012; Kuehn, 2015).

From a neuroscience perspective, the phase lag parameter β of the adaptation function sin ( ϕ i μ − ϕ j μ − β ) can also be called plasticity parameter (Aoki and Aoyagi, 2009) which accounts for different adaptation rules that may occur. Depending on the value of β the adaptation rule can be symmetric, i.e., with a cosine shape (β = π/2), or causal, i.e., with a sine shape (β = π). Symmetric as well as causal relationship are well-known forms for spike timing-dependent plasticity in neuroscience (Maistrenko et al., 2007; Caporale and Dan, 2008; Popovych et al., 2013; Lücken et al., 2016; Röhr et al., 2019). The shape of the adaptation function for different choices of the parameter β is provided in Figure 2. By varying β from 0.4π to π, we can see that the maximum of the coupling term − sin(Δϕ − β), where Δ ϕ ≡ ϕ i μ − ϕ j μ , shifts from Δϕ = −0.1π to Δϕ = 0.5π. Thus, for β = 0.5π we have a Hebbian adaptation rule where the coupling term gives a maximum positive feedback for synchronization (fire together, wire together), while for β ≠ 0.5π the feedback is asymmetric, i.e., maximum positive feedback occurs for some phase lag ϕ i μ − ϕ j μ = β − 0.5 π . Thus the adaptation lag β = 0.5π seems to be most favorable for synchronization. For β = π the coupling term is zero for synchronization, negative for ϕ i μ < ϕ j μ and positive for ϕ i μ > ϕ j μ , i.e., the coupling weight κ i j μ , and hence the input from node j to node i, is increased if ϕ i μ > ϕ j μ , i.e., if the ith oscillator is advancing the jth, and vice versa. The parameter β plays an essential role in the model because it governs the adaptivity rule of the cytokines. It will be called age parameter, since it mimics a systemic sum parameter which accounts for different influences such as physiological changes due to age, inflammaging, systemic and local inflammatory baseline, adiposity, pre-existing illness, physical inactivity, nutritional influence, etc.

In the following we use a simplified model, where the natural frequencies of both layers are identical and set to zero in a co-rotating frame: ω ¹ = ω ² = 0. Further we assume phase lag parameters α ¹¹ = α ²², and α ¹² = α ²¹ = α throughout the article. The matrix elements a i j 1 ∈ { 0,1 } of the adjacency matrix A in the parenchymal layer are chosen as a i j 1 = 1 if i ≠ j (global coupling).

In (Berner et al., 2019a; Berner et al., 2019b) it has been shown that complex heterogeneous dynamical states such as multifrequency clusters may emerge in a self-organized way in networks of adaptively coupled dynamical systems, for instance, phase oscillators. It is even more surprising that these states arise in systems with homogeneous sets of parameters and simple coupling structure (Kasatkin et al., 2017; Berner et al., 2019b; Berner et al., 2020b). In addition to the plethora of dynamical states, adaptivity also induces a high degree of multistability (Maistrenko et al., 2007). In this study, we build on the findings from (Sawicki et al., 2022) and extend these in order to understand certain parameter dependencies for the emergence of sepsis.

We assume that all cells possess the same natural frequency. To model the initial state for the potential occurrence of sepsis, we introduce a fixed initial perturbation of the cytokine activity in the immune layer representing a systemic immune response, see Figure 3. We study the effect of this initial system perturbation on the emergence of the healthy state, i.e., synchronization, in dependence of the age parameter β. Under certain conditions depending on various parameters summarized by β (age, inflammaging, chronic inflammation, other basic diseases, obesity, smoking, lack of exercise, gene polymorphisms) the unregulated cytokine expression can progress into the parenchyma and desynchronize it. In these cases, the healthy (synchronized) state is not resilient anymore against the perturbation of the immune layer.

Further, we analyze how this dependency changes depending on other parameters that shape the interaction between the parenchyma (layer 1) and the immune system (layer 2), namely the interlayer coupling strength σ, the form of the initial immune layer activation expressed by the size of the perturbation 1 < C < N, and the interlayer coupling phase lag α ¹² = α ²¹ = α. The latter parameter accounts for a delay in the layers’ interaction where α = 0 can be regarded as instantaneous.

In order to quantitatively characterize the dynamical collective state of the two-layer network, in particular its degree of frequency and phase synchronization, we introduce several measures. If the frequency of all adaptively coupled phase oscillators is the same, the phases may still be different. They may either be all the same (complete in-phase synchronization) or they may be phase-locked such that each phase oscillator oscillates with the same frequency but a fixed, time-independent phase difference. A special case is a splay state, where the phase differences of all oscillators average out, for instance if the phase of the jth oscillator is 2πj/N, j = 1, … , N. In systems of the form Eqs 1, 2, it is possible to find in-phase synchronization and splay states, and they may be interpreted as different quality of synchronization (Berner et al., 2020a). In our set-up a splay state is interpreted as a more vulnerable collective state where small perturbations can quickly lead to partial or complete desynchronization.

First, we introduce the mean phase velocities of the oscillators j in both layers μ = 1, 2 ⟨ ϕ ̇ j μ ⟩ = 1 T ∫ t t + T ϕ ̇ j μ t ′ d t ′ = ϕ j μ t + T − ϕ j μ t T (3) with averaging time window T, and the spatially averaged mean phase velocity (frequency) for each layer ω ̄ μ = 1 N ∑ j = 1 N ⟨ ϕ ̇ j μ ⟩ . In case of frequency synchronized states within the layers, we further consider a classical measure for the phase coherence within each layer, namely, the Kuramoto-Daido order parameter (Kuramoto, 1984; Daido, 1994). In particular, we look at the second moment of the order parameter R 2 μ as it is the most suitable characteristic for these kinds of patterns in adaptive networks as shown in (Berner et al., 2019a; Berner et al., 2020a). This measure of phase coherence is given by R 2 μ t = 1 N ∑ j = 1 N e i 2 ϕ j μ t . (4)

It takes values 0 ≤ R 2 μ ≤ 1 , where the lowest and the highest coherence correspond to 0 and 1, respectively. We recall that for R 2 μ = 0 we call a state a splay state and for R 2 μ = 1 an antipodal state (Berner et al., 2021c). A well-known example of a splay state is a state with fixed phase difference of 2π/N between neighboring oscillators on a ring network of N phase oscillators. Further we note that in-phase and anti-phase synchronized states are included in the class of antipodal states. We emphasize that splay states are still frequency synchronized, and hence are considered as healthy states, however, due to their weaker phase coherence properties, they may be considered as more vulnerable and less resilient than in-phase synchronized states.

Furthermore, for both layers μ = 1, 2 we calculate the ensemble average s ^μ (ensemble size N _E with ensemble elements E) of the standard deviation σ χ ( ω ̄ μ ) = 1 N ∑ j = 1 N ( ⟨ ϕ ̇ j μ ⟩ − ω ̄ μ ) 2 of the mean phase velocities s μ = 1 N E ∑ E σ χ ω ̄ E μ , (5) and the ensemble average of the corresponding normalized standard deviation σ χ ( ω ̄ E μ ) ω ̄ E μ . If the latter quantities are non-zero, they indicate the formation of frequency clusters, where the respective layer splits into clusters with different frequencies, which is indicative of a pathological state. The ensemble average is necessary to account for the multistable nature of the system, i.e., for random initial conditions some of the simulations may give a pathological state, while some may still give a healthy state. This is similar to the real physiological situation where some patients will develop sepsis, while some will not.

We further introduce another complementary measure to quantify the occurrence of pathological states, which we call the frequency cluster ratio. The frequency cluster ratio f ^μ is defined as the ratio between the number of frequency clusters N f μ in layer μ found for an ensemble of initial conditions and the size of the ensemble N _E , i.e., f μ = N f μ / N E . We consider an asymptotic state to be a frequency cluster (desynchronized, pathological state) if there exist one or more nodes j ∈ {1, … , N} such that ⟨ ϕ ̇ j μ ⟩ ≠ ω ̄ μ (deviating frequencies).

Table 1 summarizes the dynamical variables, parameters and measures of the model. In the right column the physiological meaning of all quantities is given in a concise manner. For more details on the pathological interpretation, we refer the reader to (Sawicki et al., 2022).

For the parameter scans presented in the subsequent sections, we simulate system (1)–2) for each set of parameters for the same ensemble of random initial conditions.

In this subsection we investigate the influence of the interlayer coupling strength σ on the emergence of sepsis. The interlayer coupling strength appears naturally as an important parameter in order to understand the mechanism acting during the progression of sepsis. In fact, proinflammatory cytokines act on endothelial cells and hence cause an increased blood vessel leakiness (Egger, 2005). As a result, more immune cells and cytokines enter the stroma, which consequently enhances the immune-parenchymal interaction.

In the following, we present simulation results in the (β, σ)-plane showing that after an initial cytokine perturbation in the immune layer either the healthy frequency-synchronized state is likely to be restored (dark shading), or the system is more likely to transition to a pathological desynchronized multifrequency cluster state (light shadings).

The top panels of Figure 4 depict the ensemble average s ^μ of the standard deviation of the spatially averaged mean phase velocities, which measures the average frequency desynchronization, corresponding to the amount of heterogeneous activity in system. A high or low degree of desynchronization represents a pathological or healthy physiological condition, respectively. Splitting into frequency clusters corresponds to a pathological state of the parenchyma (μ = 1) or activation of the immune layer (μ = 2). Figure 4 shows three regimes of the coupling strength σ for which the system behaves qualitatively different. Within the first regime (σ < 0.5), the parenchyma, Figure 4A, evolves independently of the immune system, Figure 4B. This can be concluded from the different values of the average activity s ¹ and s ². In fact, for sufficiently large age parameter β the initial perturbation of the immune layer leads to persistent desynchronization (activation) of the immune layer. As shown in previous work (Sawicki et al., 2022), the average activity s ² increases with increasing age parameter. Up to the critical value σ _c ≈ 0.5, the parenchyma synchronizes in most of the simulations independently of the age parameter, hence no organ-threatening desynchronization s ¹ occurs. It is worth mentioning that below but near the critial value σ _c , the boundaries between low and high activity in the immune layer become more complex due to the increasing interaction of the immune system with the parenchyma.

In the second regime, i.e., in the interval of approximately 0.5 < σ < 0.8, the systems starts to show interlayer phase locking, i.e., ϕ i 1 ( t ) − ϕ i 2 ( t ) ≈ Δ i ∈ [ 0,2 π ) for all times t. We observe that beyond the threshold of σ _c , the parenchyma may also desynchronize depending on the age parameter β. The average desynchronization s ¹ in the parenchyma and hence the potential for organ failure increases with increasing age parameter. For constant σ there always exists a threshold of the age parameter above which the parenchyma is dynamically able to desynchronize. With increasing σ the threshold shifts to larger values of β.

In the third regime of the interlayer coupling strength (σ > 0.8), the threshold of β above which the parenchyma may desynchronize does not shift further to larger values, but remains approximately fixed. Hence, we observe a clear separation in terms of the age parameter between parameter regions with healthy and regions with pathological dynamics.

In order to support our conclusions drawn from the frequency desynchronization measure s ^μ , we also plot the ratio f ^μ of simulations yielding frequency clusters divided by the total number of simulations N _E for an ensemble of N _E = 50 random initial conditions in the bottom panels of Figure 4. We observe that indeed a high value of s ^μ correlates with a higher probability of finding a frequency cluster. Therefore, both measures can be used interchangeably.

In Figure 5, we plot representative asymptotic states for different values of σ and β. They correspond to parameter values marked by letters A, B, C, D, E in Figure 4. The left and right columns show snapshots of the cytokine activity matrices κ i j 1 (parenchymal layer) and κ i j 2 (immune layer), respectively. The second column shows the mean phase velocities (average frequencies) ⟨ ϕ ̇ j μ ⟩ of the oscillators. The third column shows snapshots of the instantaneous phases ϕ j μ , and the fourth column depicts space-time plots of the phases ϕ j μ ( t ) visualizing the oscillations. We observe that depending on the choice of parameters different dynamical states emerge. In Figure 5A an in-phase synchronized state is presented showing that the system is capable of evolving into a healthy state after an initial perturbation of the immune layer. All mean phase velocities in the parenchyma and in the immune layer (collective frequencies) are the same (second column), and the oscillators in each layer are in phase (third column). The space-time plot shows spatially homogeneous periodic oscillations. The adaptive coupling weights both in the parenchymal and the immune layer are homogeneous and all weights are equal to unity (left and right columns). Another completely healthy state is shown in Figure 5C’ where instead of an in-phase synchronized state a splay state is formed in both layers (third column), i.e., the order parameter R 2 μ = 0 for both layers, but the frequencies are still the same (second column). The space-time plot (fourth column) shows traveling waves, rather than spatially homogeneous oscillations as in panel A. In (Sawicki et al., 2022), we have speculated that this type of synchronized states can be interpreted as a vulnerable state emerging in coexistence with pathological states. Indeed, Figure 5C shows a pathological frequency cluster state for the same parameters but different initial conditions. Here both the parenchyma and the immune layer exhibit a two-frequency cluster state, where a smaller cluster with lower frequency splits off from the large cluster (marked by a small red circle in the second column). The small clusters can also be clearly seen in the snapshots of the phases (third column), in the perturbations of the space-time pattern (fourth column), and in the lighter red color in the cytokine matrices (left and right columns). Figures 5D,E show states that can still be regarded as healthy from the perspective that the parenchymal nodes are in synchrony, while the immune layer remains activated after the initial perturbation and exhibits small clusters of deviating frequencies. This shows up in the mean phase velocity profiles (second column), in the snapshots of the phases (third column), in the space-time plot (fourth column) and in the cytokine matrix of the immune layer (right column). These states demonstrate a high degree of parenchymal resilience to the persistent activation of the immune layer. A pathological state is also presented in Figure 5B. Here, the parenchymal layer shows desynchronization and a frequency cluster (small red circle), as well, which may be considered as the starting point of an organ failure.

The splay states with R 2 μ = 0 for both layers (panel C’) represent a special class of healthy states. In particular, due to their structure, the oscillators in this state effectively decouple and are potentially more vulnerable to external perturbations. Moreover, as also shown in (Sawicki et al., 2022), these states may coexist with frequency clusters. In order to quantify this observation, we plot the probability of finding a splay state in dependence on σ and β in Figure 6. By comparing Figures 4A,B with Figure 6, we see that the regions for the existence of splay states have large overlap with the region of pathological cluster states of the paranchyma (yellow hatched area). It should, however, be noted that for intermediate values of the interlayer coupling strength and the age parameter there exists a large region in parameter space for which frequency clusters are very likely, whereas almost no splay states can be found.

Figure 7A presents a cut through the parameter plane of Figure 4C at coupling strength σ = 1. It shows that the probability of a frequency cluster, i.e., a pathological sepsis state, sharply rises with age parameter β above approximately β > 0.5π. This curve compares favorably with empirical data of patients which gives the number of cases of sepsis per 100 000 inhabitants in Germany as a function of age, presented in Figure 7B.

In this section, we have numerically analyzed the dependence of sepsis on the interlayer coupling strength and the age parameter after an initial perturbation of the immune system. We have identified three regimes with qualitatively different dynamics. First, below a critical coupling strength, the healthy state is preserved for all values of the age parameter. Second, above the critical coupling strength the probability of sepsis sharply rises with increasing age parameter β above a threshold of β, and the threshold itself increases with increasing coupling strength. In the third regime this threshold saturates at a fixed value of β. This means that in a certain intermediate coupling range stronger coupling to the immune layer can preserve the healthy state even at larger age parameter, but eventually the age threshold cannot be shifted further, and the pathological state cannot be avoided. It also implies that an interlayer coupling weight slightly above a critical value could be potentially threatening for patients with a wide range of age parameters, in particular also “younger” patients, i.e., with smaller values of β. This threat, however, shifts to higher values of β as the coupling strength between the layers is increased. Remarkably, our simulations show that, depending upon the initial conditions, healthy states coexist with pathological states for the same parameter values, indicating that the outcome of sepsis after an initial perturbation of the immune system cannot be straightforwardly predicted.

In this section, we analyze the dependence of sepsis on the interlayer phase lag parameter α. In particular, we investigate the robustness of our results from the previous subsection with respect to this parameter. Phase lags have been used to account for interaction delays (Madadi Asl et al., 2018; Sawicki, 2019) and are known to be critical for the emergence of complex dynamics (Omel’chenko et al., 2010; Omelchenko et al., 2013; Omel’chenko, 2018; Omel’chenko and Knobloch, 2019; Gerster et al., 2020; Schöll, 2020; Schöll, 2021). Motivated by the results presented in the previous subsection for the case α = 0, we choose an interlayer coupling strength σ for which sepsis may occur. Therefore, we set σ = 1 throughout this subsection.

In Figure 8, we show the ensemble average s ^μ as a measure the average frequency desynchronization for both layers (top left and top right, respectively). In the bottom panels we plot the corresponding ratio f ^μ of simulations yielding frequency clusters divided by the total number of simulations N _E for an ensemble of N _E = 50 random initial conditions. The behavior of the parenchyma (μ = 1, left panel) and the immune layer (μ = 2, right panel) is practically the same. From the figure, we see that for small values of α, the threshold in the age parameter for the occurrence of sepsis is only slightly changed. It should be noted that with increasing but small interlayer layer phase lag the β threshold does not change much, but the transition from the healthy state to the pathological state becomes sharper, i.e., the frequency cluster ratio increases more sharply. A dramatic change of the behavior occurs slightly below α = π/4, which is also the value of the phase lag where in single-layer networks complex partial synchronization patterns of chimera-type are found (Omel’chenko et al., 2010; Omelchenko et al., 2013). For larger values of α > π/4, we observe that the dependence upon β flips, and the desynchronized (activated) state occurs with some probability for lower β, while for higher β the healthy synchronized state is observed. At approximately α ≈ 0.42π another flip occurs, and with increasing β there is once more a pronounced transition from the synchronized state to a desynchronized frequency cluster state at a distinct threshold of β, which decreases with further increasing α. This alternating behavior is due to the periodic nature of the coupling function sin ( ϕ i 1 − ϕ i 2 + α ) . It indicates that the regime which corresponds to physiological conditions and to our interpretation of β as age parameter seems to be confined to α < π/4, but within this interval the observed behavior is robust. Supplementary Figure S1 of the Supplemental Material depicts the map of regimes for a larger range of α ∈ [0, 2π]. This clearly shows the structure of the tongues of two-cluster states (bright colors), which obeys a π-periodic pattern in α.

Figure 9A shows details of the dynamics for an exemplary parameter set α = 0.2 π, β = 0.58 π, in a plot similar to Figure 5. Comparing it with Figure 5B where α = 0, but the other parameters are the same, we see that our model is robust with respect to the parameter α = 0.

This section is devoted to study the impact of the initial perturbation in the immune system corresponding to cytokine activation. For this, we vary the cluster size C of the initial condition of the adaptive coupling weight matrix presented in Figure 3. Here, we choose the two other parameters as σ = 1 and α = 0.

We see in Figure 10 that independently of the initial cluster size a transition from the healthy synchronized state to the pathological desynchronized state in the parenchyma may be observed, and surprizingly the threshold β _c is insensitive to the size C of the initial perturbation in a wide range from only a few cells to half the immune system (C/N = 0.5). The behavior of the parenchyma (μ = 1, left panel) and the immune layer (μ = 2, right panel) is very similar.

Panels B, C of Figure 9 show details of the dynamics for C/N = 0.1 (C = 20) and C/N = 0.4 (C = 80), respectively, in a plot similar to Figure 5. Comparing it with Figure 5B where C = 40, but the other parameters are the same, we see that the asymptotic state does not depend upon the size of the initial perturbation. This finding seems to be in line with the medical observation that there is no direct relation between the cause and form of an inflammatory response and the frequency of occurrence of sepsis.

The adaptive network model Eqs 1, 2 can also be written in the form of two integral equations with an exponential kernel for the phases in the two layers ϕ i 1 and ϕ i 2 by using a Green’s function technique to eliminate the differential equations for the adaptive coupling weights κ i j μ . The solution of the general inhomogeneous differential equation κ ̇ i j = − ϵ κ i j + sin ϕ i − ϕ j − β (6) is given by the integral κ i j t = − ϵ ∫ 0 ∞ d s e − ϵ s sin ϕ i t − s − ϕ j t − s − β (7)

Hence the adaptive two-layer phase oscillator model in the co-rotating frame (ω ¹ = ω ² = 0) with α ¹¹ = α ²² = α ⁰ and α ¹² = α ²¹ = 0 is reduced to: ϕ ̇ i 1 = − 1 N ∑ j = 1 N a i j 1 − ϵ 1 ∫ 0 ∞ d s e − ϵ 1 s ⁡ sin ϕ i 1 t − s − ϕ j 1 t − s − β sin ϕ i 1 − ϕ j 1 + α 0 − σ ⁡ sin ϕ i 1 − ϕ i 2 , (8) ϕ ̇ i 2 = 1 N ∑ j = 1 N ϵ 2 ∫ 0 ∞ d s e − ϵ 2 s ⁡ sin ϕ i 1 t − s − ϕ j 1 t − s − β sin ϕ i 2 − ϕ j 2 + α 0 − σ ⁡ sin ϕ i 2 − ϕ i 1 , (9)

The adaptation function sin ( ϕ i μ − ϕ j μ − β ) shown in Figure 2 now enters as a distributed time delayed feedback which contains the whole history. For the completely synchronized (healthy) state ϕ i μ = ϕ j μ = ϕ μ this term can be integrated out, using ϵ ∫ 0 ∞ d s e − ϵ s = 1 and 1 N ∑ j = 1 N = 1 , and setting a i j 1 = 1 (for N − 1 ≈ N): ϕ ̇ 1 = − 1 + sin ⁡ β sin α 0 − σ ⁡ sin ϕ 1 − ϕ 2 , (10) ϕ ̇ 2 = − sin ⁡ β ⁡ sin α 0 − σ ⁡ sin ϕ 2 − ϕ 1 , (11)

The condition for frequency synchronization ⟨ ϕ ̇ 1 ⟩ = ⟨ ϕ ̇ 2 ⟩ yields a condition for the phase lag between the two layers 1 and 2 sin ϕ 1 − ϕ 2 = − sin α 0 2 σ (12) which agrees with the numerical simulations in Figure 5A (ϕ ¹ − ϕ ² = 0.126π). It follows from Eq. (12) that σ > | sin α 0 | 2 is a condition for the existence of the fully in-phase synchronized state in both layers, e.g., σ > 0.385 for α ⁰ = −0.28π.

For cluster states in either the immune layer, or in both layers, the situation is more complicated. If a large synchronized cluster with ⟨ ϕ ̇ i 1 ⟩ = ω L coexists with a smaller cluster of different frequency ⟨ ϕ ̇ j 1 ⟩ = ω L − Δ ω and desynchronized phases θ _j , Eqs 10, 11 must be supplemented for the large cluster (i ∈ L) and the small cluster (j ∈ S) by correction terms. These are complicated temporally oscillating functions, and condition (12) is modified for frequency synchronization of the large cluster at frequency ω _L and for the small cluster at ω _L − Δω. By temporal averaging over trigonometric functions, one may obtain rough approximations. Assuming slow adaptation ϵ ¹ and ϵ ², and inserting the phases ϕ i 1 = ω L t for i ∈ L and ϕ j 1 = ( ω L − Δ ω ) t + θ j for j ∈ S, for non-synchronous solutions i, j the fast oscillating terms in the integrals average out to zero. Thus more detailed expressions for the regime of existence of frequency cluster states as a function of σ and β may be derived.