Critical Parameters in Dynamic Network Modeling of Sepsis

1 Introduction

The Systemic Inflammatory Response Syndrome (SIRS) is a life-threatening organ dysfunction, which is induced by infectious pathogens or endogenous antigens. It is an induced disease of the innate immune system. Because of its complexity no detailed model is available. Sepsis, which is the largest subclass of SIRS, is defined as an infect-induced organ failure, where, however, only in 30–40% of all cases the pathogen can be identified. Organs far from the location of primary infection are disturbed in their proper function by the host reaction (Singer et al., 2016). The lethality of sepsis or septic shock in spite of high-performance medicine is as high as 45% in intense care units (hospital mortality) and 74% after 48 months (Schmidt et al., 2020). A gold standard for diagnosis of sepsis is still missing (Brunkhorst et al., 2020), clinical diagnosis rests upon infection-related organ dysfunction of lungs, kidney, liver, circulatory system, blood count, or central nervous system (Schmidt et al., 2020). In a particular case with severe infection like pneumonia or peritonitis and the same risk factors like age, sex, and underlying medical condition, it cannot be predicted whether the patient will survive unscathed, or whether the infection progresses and ends up lethally within a short time by multi-organ failure. There exists a highly individual inflammatory reaction of the host, and a specific therapy for pro-inflammatory dysregulation is not available (Weis et al., 2017).

The organ-damaging host reaction is caused by dysregulatory, pro-inflammatory cytokines. This condition is known as cytokine storm. The organ damage resulting from this can occur sequentially or simultaneously in several organs, and it may be mild, moderate, or severe. Clinically, the organ functioning is monitored and rated on a daily basis in terms of the four-stage Sepsis-related Organ Failure Assessment (SOFA) score. In case of SIRS not all organs are always and to the same extent disturbed. One may hypothesize that in each individual patient certain organs possess a more or less pronounced resilience against the cytokine storm. The aim of this work is to model the conditions for organ failure, the induced organ dysfunction, and resilience of organs, as well as the overall state of the organism after recovery.

A unified disease model with the innate immune system as reference point is the basis for our modeling approach in terms of nonlinear dynamics of complex networks. Note that this is not a biochemical or genetic or cellular or tissue model, but it rather describes the functional interplay of the immune system with parenchymal cells of the organs in terms of a simple generic model of coupled nonlinear oscillators which may exhibit coherent or incoherent collective dynamics. The role of synchronization is an important aspect in the field of network physiology, where multi-component physiological systems continuously interact in an integrated network to coordinate their functions (Bashan et al., 2012; Ivanov et al., 2014; Bartsch et al., 2015; Lin et al., 2016; Moorman et al., 2016). The structural organization and functional complexity of human organisms has been associated with phase synchronization as well as phase transitions (Chen et al., 2006; Xu et al., 2006; Ivanov et al., 2009; Bartsch et al., 2012) between different modes of synchronization in real physiological systems. In case of complex diseases, the progression from a healthy to sick state can be abrupt and may cause a critical transition (Chen et al., 2012; Liu et al., 2012; Liu et al., 2013a; Liu et al., 2013b; Shi et al., 2022).

In this paper, we employ a two-layer network model for sepsis based upon the interaction of parenchymal cells and immune cells via cytokines and the coevolutionary dynamics of parenchymal and immune cells and cytokines (Sawicki et al., 2022). Parenchyma is the bulk of functional substance in an organ or structure, in contrast to the stroma, which refers to the structural tissue of organs or structures, namely, the unspecific connective tissues. In many organs the parenchyma consists of epithelial cells. A simple paradigmatic model of adaptively coupled phase oscillators (Aoki and Aoyagi, 2009; Aoki and Aoyagi, 2011; Nekorkin and Kasatkin, 2016; Kasatkin et al., 2017; Berner et al., 2019a; Berner et al., 2020a; Berner et al., 2021a) is chosen as a first step to model the coupled dynamics of parenchymal cells and unspecific immune cells, which are represented by nodes of a duplex network. Multiplex networks are special multilayer networks (Kivelä et al., 2014), where the network layers are interconnected such that the layers contain the same number of nodes, and only pairwise connections between corresponding nodes from the layers exist. Duplex networks are simple representatives of multiplex networks consisting of two layers, and are known to exhibit complex synchronization scenarios (Leyva et al., 2018; Sawicki et al., 2018; Sawicki, 2019; Berner et al., 2021b). The cytokine interaction within both layers is modeled by adaptive coupling strengths between the nodes representing the parenchymal cells (slow timescale) and between the nodes of the immune cells (faster timescale, but still slower than the timescale of the cell metabolism governed by phase oscillator dynamics). We stress that our model is not a detailed model of organs, such as for instance specific biochemical models for carcinogenesis (Vineis et al., 2010), but a functional model of dynamic interactions. Thus the cytokines are not modeled as concentrations but rather by information flow between the parenchymal layer and the immune layer describing the cytokine activity. In both layers the base topology is global (all-to-all) coupling; while the coupling in the parenchymal layer has a fixed time-independent contribution and an adaptive time-dependent contribution modeling the cytokine activity, the coupling in the immune layer is only adaptive. In further work more sophisticated local dynamics, such as activator-inhibitor kinetics of FitzHugh-Nagumo type with two variables (fast activator and slow inhibitor, respectively) and more elaborate network topologies might be chosen.

In this article, we analyze the parameter dependency of an organ-damaging interaction between the dysregulated immune system and the parenchyma. In terms of our proposed functional model, we investigate the emergence of synchronization and frequency synchronization in the parenchyma, which is related to a healthy and unhealthy condition of the organ system, respectively. In (Sawicki et al., 2022), it has been shown that an initially activated immune system may induce an activation of the parenchyma, i.e., emergence of frequency clusters, or leave the parenchyma unaffected depending on a patient’s individual characteristics. The present study provides insight into the robustness of the emerging pathological state with respect to changes of parameters. We utilize a numerical analysis to find critical parameters that are crucial for the interaction of the immune system with the parenchyma. We shed further light on the question how a dysregulated immune system triggers the onset of organ failure.

The article is organized as follows: In Section 2 we provide a pathophysiological description of sepsis. In Section 3 we introduce the functional model that we employ for the analysis of sepsis. Section 4 gives a systematic survey of critical parameters of sepsis in our model simulations. Finally, in Section 5 we draw conclusions.

2 Pathophysiological Description of Sepsis

2.1 Innate Immune System

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.

2.2 Relaps

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).

3 Model

In this section, we introduce the functional model that we employ for the analysis of sepsis. We introduce all parameters and variables and provide details on the measures used to analyze the system.

3.1 Schematic Sepsis Model

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

FIGURE 1. Schematic illustration of the sepsis model. (A) A tissue element is depicted, in which the basic processes of sepsis take place: shown are the cells (colored) involved such as parenchymal, fibroblast, endothelial cells, and macrophages, polymorphonuclear leukocytes and thrombocytes in the parenchyma (grey), the stroma (yellow), and the capillary blood vessel. (B) depicts the functional interactions within and between the two corresponding network layers in our model, the parenchyma and the stroma (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.

3.2 Functional Two-Layer Network Model

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 ${\varphi }_i^1$, i = 1, … , N and the network layer of immune cells (superscript 2) are presented by N phase oscillators ${\varphi }_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 ${\kappa }_{ij}^1$, and those between the immune cells by coupling weights ${\kappa }_{ij}^2$. Note that ${\varphi }_i^2$ and ${\kappa }_{ij}^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:

$\begin{array}{ll}\hfill {\dot{\varphi }}_i^1& ={\omega }^1-\frac{1}{N}\sum _{j=1}^N\left(a_{ij}^1+{\kappa }_{ij}^1\right)\sin \left({\varphi }_i^1-{\varphi }_j^1+{\alpha }^{11}\right)-\sigma \sin \left({\varphi }_i^1-{\varphi }_i^2+{\alpha }^{12}\right),\hfill \\ \hfill {\dot{\kappa }}_{ij}^1& =-{ϵ}^1\left({\kappa }_{ij}^1+\sin \left({\varphi }_i^1-{\varphi }_j^1-\beta \right)\right),\hfill \end{array}(1)$

$\begin{array}{ll}\hfill {\dot{\varphi }}_i^2& ={\omega }^2-\frac{1}{N}\sum _{j=1}^N{\kappa }_{ij}^2\sin \left({\varphi }_i^2-{\varphi }_j^2+{\alpha }^{22}\right)-\sigma \sin \left({\varphi }_i^2-{\varphi }_i^1+{\alpha }^{21}\right),\hfill \\ \hfill {\dot{\kappa }}_{ij}^2& =-{ϵ}^2\left({\kappa }_{ij}^2+\sin \left({\varphi }_i^2-{\varphi }_j^2-\beta \right)\right),\hfill \end{array}(2)$

where ${\varphi }_i^{\mu }\in [\mathrm{0,2}\pi )$ 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_{ij}^1\in [\mathrm{0,1}]$ (fixed interaction within an organ) and ${\kappa }_{ij}^{\mu }\in [-\mathrm{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 ({\varphi }_i^{\mu }-{\varphi }_j^{\mu }-\beta )$ 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 $\mathrm{\Delta }\varphi \equiv {\varphi }_i^{\mu }-{\varphi }_j^{\mu }$, 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 ${\varphi }_i^{\mu }-{\varphi }_j^{\mu }=\beta -0.5\pi$. Thus the adaptation lag β = 0.5π seems to be most favorable for synchronization. For β = π the coupling term is zero for synchronization, negative for ${\varphi }_i^{\mu }<{\varphi }_j^{\mu }$ and positive for ${\varphi }_i^{\mu }>{\varphi }_j^{\mu }$, i.e., the coupling weight ${\kappa }_{ij}^{\mu }$, and hence the input from node j to node i, is increased if ${\varphi }_i^{\mu }>{\varphi }_j^{\mu }$, 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.

FIGURE 2

FIGURE 2. Illustration of the adaptation function in dependence of the age parameter β.

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_{ij}^1\in \{\mathrm{0,1}\}$ of the adjacency matrix A in the parenchymal layer are chosen as $a_{ij}^1=1$ if i ≠ j (global coupling).

3.3 Methods of Analysis

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.

FIGURE 3

FIGURE 3. Initial conditions of sepsis: Cytokine dysregulation expressed by a cluster structure of the cytokine activity matrix κ² imposes a systemic activation of the immune layer representing the beginning of sepsis. The figure shows an initial condition used for simulations of (1)–(2) with N = 200. The initial cytokine activities in the parenchymal layer ${\kappa }_{ij}^1$ and the initial phases in both layers are randomly drawn from a uniform distribution on the interval [ − 1, 1] and the interval [0, 2π], respectively. The cytokine activities in the immune layer ${\kappa }_{ij}^2$ are initially given by a two-cluster structure where the smallest cluster has size C. The cytokine activities ${\kappa }_{ij}^2$ are 1 within and 0 between the clusters.

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

$⟨{\dot{\varphi }}_j^{\mu }⟩=\frac{1}{T}{\int }_t^{t+T}{\dot{\varphi }}_j^{\mu }\left(t^{\prime }\right)\mathrm{d}{t}^{\prime }=\frac{{\varphi }_j^{\mu }\left(t+T\right)-{\varphi }_j^{\mu }\left(t\right)}{T}(3)$

with averaging time window T, and the spatially averaged mean phase velocity (frequency) for each layer ${\bar{\omega }}^{\mu }=\frac{1}{N}{\sum }_{j=1}^N⟨{\dot{\varphi }}_j^{\mu }⟩$. 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^{\mu }$ 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^{\mu }\left(t\right)=\frac{1}{N}\left|\sum _{j=1}^N{e}^{\mathrm{i}2{\varphi }_j^{\mu }\left(t\right)}\right|.(4)$

It takes values $0\le R_2^{\mu }\le 1$, where the lowest and the highest coherence correspond to 0 and 1, respectively. We recall that for $R_2^{\mu }=0$ we call a state a splay state and for $R_2^{\mu }=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 ${\sigma }_{\chi }({\bar{\omega }}^{\mu })=\sqrt{\frac{1}{N}{\sum }_{j=1}^N{(⟨{\dot{\varphi }}_j^{\mu }⟩-{\bar{\omega }}^{\mu })}^2}$ of the mean phase velocities

$s^{\mu }=\frac{1}{N_E}\sum _E{\sigma }_{\chi }\left({\bar{\omega }}_E^{\mu }\right),(5)$

and the ensemble average of the corresponding normalized standard deviation $\frac{{\sigma }_{\chi }({\bar{\omega }}_E^{\mu })}{{\bar{\omega }}_E^{\mu }}$. 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^{\mu }$ in layer μ found for an ensemble of initial conditions and the size of the ensemble N_E, i.e., $f^{\mu }=N_f^{\mu }/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 $⟨{\dot{\varphi }}_j^{\mu }⟩\ne {\bar{\omega }}^{\mu }$ (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).

TABLE 1

TABLE 1. Physiological meaning of the dynamical variables, parameters and measures of the model (superscripts referring to layers μ = 1 and μ = 2 omitted).

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.

4 Critical Parameters for Sepsis

This section is devoted to the numerical analysis of critical parameters controlling the interaction of the parenchyma with the immune system, i.e., σ and α, and the initial activation of the immune system, i.e., activation cluster size C, see Figure 3. In the following, we analyze the impact of these parameters in addition to the age parameter β that has been found to be crucial for the description of the patient’s physiological condition (Sawicki et al., 2022).

4.1 The Interlayer Interaction Strength as a Critical Parameter for Modeling Sepsis

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.

FIGURE 4

FIGURE 4. Map of regimes: ensemble average s^μ of the standard deviation of the spatially averaged mean phase velocities (A,B) and the frequency cluster ratio (C,D) in the parameter plane of age parameter β and interlayer coupling strength σ for the parenchymal (A,C) and immune layer (B,D), respectively. Bright colors correspond to the formation of frequency clusters. Simulation parameters: N = 200, α¹¹ = α²² = −0.28π, α¹² = α²¹ = 0, $a_{ij}^1=1$, ϵ¹ = 0.03, ϵ² = 0.3, ω¹ = ω² = 0, C = 40. Ensemble size is N_E = 50. The simulation time is 2000 time units, the averaging time window 1000.

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., ${\varphi }_i^1(t)-{\varphi }_i^2(t)\approx {\mathrm{\Delta }}_i\in [\mathrm{0,2}\pi )$ 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 ${\kappa }_{ij}^1$ (parenchymal layer) and ${\kappa }_{ij}^2$ (immune layer), respectively. The second column shows the mean phase velocities (average frequencies) $⟨{\dot{\varphi }}_j^{\mu }⟩$ of the oscillators. The third column shows snapshots of the instantaneous phases ${\varphi }_j^{\mu }$, and the fourth column depicts space-time plots of the phases ${\varphi }_j^{\mu }(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^{\mu }=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.

FIGURE 5

FIGURE 5. Details of dynamics for healthy parenchymal states without clusters (A,C’,D,E) and a pathological parenchymal state with frequency clusters (B,C) for different values of β chosen as in Figure 4. The states shown represent a healthy state in (A) (β = 0.5 π, σ = 1), a pathological state in (B) (β = 0.58 π, σ = 1) and (C) (β = 0.7 π, σ = 1), where a red circle marks the small pathological cluster, a healthy but vulnerable state in (C’) (β = 0.7 π, σ = 1) and two resilient states in (D) (β = 0.5 π, σ = 0.2) and (E) (β = 0.65 π, σ = 0.45). The left and right columns show snapshots of cytokine activity matrices ${\kappa }_{ij}^1$ (parenchymal layer) and ${\kappa }_{ij}^2$ (immune layer), respectively (color coded). Second column: mean phase velocities (average frequencies) $⟨{\dot{\varphi }}_j^{\mu }⟩$ of the oscillators. Third column: snapshots of phases ${\varphi }_j^{\mu }$. The parenchymal nodes are labeled j = 1, … , 200, and the immune nodes are labeled j = 201, … , 400. Within each layer μ the nodes are sorted first by $⟨{\dot{\varphi }}_j^1⟩$, then by ${\varphi }_j^1$, respectively. Fourth column: space-time plot of phases ${\varphi }_j^{\mu }(t)$ (color coded). All parameters are chosen as in Figure 4.

The splay states with $R_2^{\mu }=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 6

FIGURE 6. Probability of finding a splay state (splay ratio), see Figure 5, from N_E = 50 random initial conditions, plotted in the parameter plane of interlayer coupling strength σ and the age parameter β for the parenchymal (A) and the immune layer (B). The yellow hatched area shows schematically the regime of pathological cluster states in the parenchyma. Data taken from simulation shown in Figure 4.

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.

FIGURE 7

FIGURE 7. Qualitative comparison of model prediction with empirical data. (A) Frequency cluster ratio for σ = 1 vs. age parameter β, for the parameters in Figure 4, where all data points were averaged over a sliding window of 4 neighboring data points. (B) Empirical data taken from (Fleischmann et al., 2016) showing the hospitalization incidence of sepsis per 100 000 inhabitants in Germany by age group for the years from 2007 to 2013.

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.

4.2 The Interlayer Phase Lag as a Critical Parameter for Modeling Sepsis

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 ({\varphi }_i^1-{\varphi }_i^2+\alpha )$. 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 8

FIGURE 8. Map of regimes: ensemble average s^μ of the standard deviation of the spatially averaged mean phase velocities (A,B) and the frequency cluster ratio (C,D) in the parameter plane of age parameter β and interlayer interaction phase lag α for the parenchymal (A,C) and immune layer (B,D), respectively. Bright colors correspond to the formation of frequency clusters. Ensemble size is N_E = 50. Simulation parameters: σ = 1, α ≡ α¹² = α²¹; all other parameters as in Figure 4.

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.

FIGURE 9

FIGURE 9. Details of dynamics for different values of α and C. (A) α = 0.2 π, (B) C/N = 0.1, (C) C/N = 0.4. The left and right columns show snapshots of cytokine activity matrices ${\kappa }_{ij}^1$ (parenchymal layer) and ${\kappa }_{ij}^2$ (immune layer), respectively (color coded). Second column: mean phase velocities (average frequencies) $⟨{\dot{\varphi }}_j^{\mu }⟩$ of the oscillators. Third column: snapshots of phases ${\varphi }_j^{\mu }$. The parenchymal nodes are labeled j = 1, … , 200, and the immune nodes are labeled j = 201, … , 400. Within each layer μ the nodes are sorted first by $⟨{\dot{\varphi }}_j^1⟩$, then by ${\varphi }_j^1$, respectively. Fourth column: space-time plot of phases ${\varphi }_j^{\mu }(t)$ (color coded). A red circle marks the small pathological cluster in the parenchyma. Simulation parameters: σ = 1, β = 0.58 π; all other parameters are as in Figure 4.

4.3 The Immune Activation as a Critical Parameter for Modeling Sepsis

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.

FIGURE 10

FIGURE 10. Map of regimes: ensemble average s^μ of the standard deviation of the spatially averaged mean phase velocities (A,B) and the frequency cluster ratio (C,D) in the parameter plane of age parameter β and the initial immune layer perturbation expressed by the cluster size C/N for the parenchymal (A,C) and immune layer (B,D), respectively. Bright colors correspond to the formation of frequency clusters. Ensemble size is N_E = 50. Simulation parameters: σ = 1, α ≡ α¹² = α²¹ = 0; all other parameters as in Figure 4.

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.

4.4 Analytic Approximations

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 ${\varphi }_i^1$ and ${\varphi }_i^2$ by using a Green’s function technique to eliminate the differential equations for the adaptive coupling weights ${\kappa }_{ij}^{\mu }$. The solution of the general inhomogeneous differential equation

${\dot{\kappa }}_{ij}=-ϵ\left({\kappa }_{ij}+\sin \left({\varphi }_i-{\varphi }_j-\beta \right)\right)(6)$

is given by the integral

${\kappa }_{ij}\left(t\right)=-ϵ{\int }_0^{\infty }ds{e}^{-ϵs}\left(\sin \left({\varphi }_i\left(t-s\right)-{\varphi }_j\left(t-s\right)-\beta \right)\right)(7)$

Hence the adaptive two-layer phase oscillator model in the co-rotating frame (ω¹ = ω² = 0) with α¹¹ = α²² = α⁰ and α¹² = α²¹ = 0 is reduced to:

${\dot{\varphi }}_i^1=-\frac{1}{N}\sum _{j=1}^N\left(a_{ij}^1-{ϵ}^1{\int }_0^{\infty }ds{e}^{-{ϵ}^{1}s}\sin \left({\varphi }_i^1\left(t-s\right)-{\varphi }_j^1\left(t-s\right)-\beta \right)\right)\sin \left({\varphi }_i^1-{\varphi }_j^1+{\alpha }^0\right)-\sigma \sin \left({\varphi }_i^1-{\varphi }_i^2\right),(8)$

${\dot{\varphi }}_i^2=\frac{1}{N}\sum _{j=1}^N\left({ϵ}^2{\int }_0^{\infty }ds{e}^{-{ϵ}^{2}s}\sin \left({\varphi }_i^1\left(t-s\right)-{\varphi }_j^1\left(t-s\right)-\beta \right)\right)\sin \left({\varphi }_i^2-{\varphi }_j^2+{\alpha }^0\right)-\sigma \sin \left({\varphi }_i^2-{\varphi }_i^1\right),(9)$

The adaptation function $\sin ({\varphi }_i^{\mu }-{\varphi }_j^{\mu }-\beta )$ shown in Figure 2 now enters as a distributed time delayed feedback which contains the whole history. For the completely synchronized (healthy) state ${\varphi }_i^{\mu }={\varphi }_j^{\mu }={\varphi }^{\mu }$ this term can be integrated out, using $ϵ{\int }_0^{\infty }ds{e}^{-ϵs}=1$ and $\frac{1}{N}{\sum }_{j=1}^N=1$, and setting $a_{ij}^1=1$ (for N − 1 ≈ N):

${\dot{\varphi }}^1=-\left(1+\sin \beta \right)\sin {\alpha }^0-\sigma \sin \left({\varphi }^1-{\varphi }^2\right),(10)$

${\dot{\varphi }}^2=-\sin \beta \sin {\alpha }^0-\sigma \sin \left({\varphi }^2-{\varphi }^1\right),(11)$

The condition for frequency synchronization $⟨{\dot{\varphi }}^1⟩=⟨{\dot{\varphi }}^2⟩$ yields a condition for the phase lag between the two layers 1 and 2

$\sin \left({\varphi }^1-{\varphi }^2\right)=-\frac{\sin {\alpha }^0}{2\sigma }(12)$

which agrees with the numerical simulations in Figure 5A (ϕ¹ − ϕ² = 0.126π). It follows from Eq. (12) that $\sigma >\frac{|\sin {\alpha }^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 $⟨{\dot{\varphi }}_i^1⟩={\omega }_L$ coexists with a smaller cluster of different frequency $⟨{\dot{\varphi }}_j^1⟩={\omega }_L-\mathrm{\Delta }\omega$ 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 ${\varphi }_i^1={\omega }_{L}t$ for i ∈ L and ${\varphi }_j^1=({\omega }_L-\mathrm{\Delta }\omega )t+{\theta }_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.

5 Conclusion

Within the framework of network physiology, we have proposed a functional model of coupled dynamical systems which is able to describe healthy states as well as pathological states related to sepsis. Sepsis is a life threatening pathological state that can potentially lead to organ dysfunction and death. By using a multilayer dynamical network approach, our model describes the collective dynamics of the parenchyma and the stroma (innate immune system) as well as their interaction.

Extending previous work on a unified description of tumor disease and sepsis (Sawicki et al., 2022), we have modeled the coevolutionary adaptive dynamics of parenchymal cells, immune cells, and cytokines. By means of the simple paradigmatic model of phase oscillators in a two-layer system, we have analyzed the emergence of organ threatening interactions between the dysregulated immune system and the parenchyma. We have demonstrated that the complex cellular cooperation between the parenchymal layer and the immune layer results either in a healthy physiological (frequency synchronized) or in a pathological (desynchronized or multifrequency cluster) state in the parenchyma. Thus we have explained sepsis by the dysregulation of the healthy homeostatic state and have provided insight into the complex stabilizing and destabilizing interaction of parenchyma and immune system. The coupled dynamics of parenchymal cells (metabolism) and nonspecific immune cells (response of the innate immune system) is represented by phase oscillators in a duplex layer. The cytokine-mediated indirect communication pathways of the different cell types involved in both layers are modeled by adaptive coupling weights between nodes representing immune cells (with fast adaptation timescale) and parenchymal cells (slow adaptation timescale), and between pairs of parenchymal and immune cells in the duplex network (fixed bidirectional coupling).

In a pathophysiological context, the different scenarios obtained in our model from an initial activation of the immune system, e.g. by inflammation, can be interpreted as inflammation without organ failure (the parenchyma stays in-phase synchronized, Figure 5A), organ failure (the parenchyma forms a two-frequency cluster state, Figures 5B,C), systemic spreading into other organ systems (large-scale desynchronization, large frequency clusters), healing or parenchymal resilience to the persistent activation of the immune layer (synchronization of the parenchyma, although the immune layer forms a two-frequency cluster state, Figures 5D,E), or relapse from a vulnerable healthy state (splay-synchronized state, Figure 5C’). As critical interaction parameters we have identified the adaptation phase lag β which determines the adaptation law and is a physiological sum parameter (called age parameter), the interlayer coupling strength σ, the interlayer coupling phase lag α, and the size C of the initial perturbation of the activated immune layer cytokine coupling matrix which describes the immune system’s initial activation caused by inflammation. An adaptation phase lag β of the order of π/2 corresponds to a cosine-like adaptation function which assumes its maximum for the healthy (synchronized) state, while a larger phase lag β is related to delays in adaptability (Figure 2). Thus β ≈ π/2 can be interpreted in a physiological context as fast adaptability which is typical of young age and good physical conditions, and favors the healthy state, while larger β is not optimal for maintaining the healthy state. Regarding the size of the initial perturbation C, it should be noted that we use special initial conditions (random initial conditions of the phases in the parenchyma, the immune layer, the weighted coupling matrix of the parenchyma, and a cluster state in the coupling matrix of the immune layer, see Figure 3) which do not correspond to the healthy state (in-phase synchronization of parenchyma and immune layer). Rather, our motivation is to map out the whole dynamic state space which is characterized by multistability between the healthy state and pathological states, and the probability of observing pathological states in an ensemble of simulations depends upon these initial conditions. Of course, by choosing the healthy fully phase-synchronized state as initial condition, one could increase the number of observed healthy states.

In extensive simulations, we have analyzed the dynamics of the sepsis model in dependence on these critical parameters, and have found that particularly the age parameter β and the interlayer interaction strength σ are important model parameters for describing the emergence of pathological states. The crucial role of the age parameter has been already described in (Sawicki et al., 2022) for the emergence of tumor disease. In this study, we have shown that depending on the age parameter and the interlayer coupling strength different dynamical regimes with clear pathophysiological meaning emerge. We have mapped out parameter regimes where an initial inflammation 1) can be regulated and the systems enters a completely healthy state (healing), 2) is persistent, i.e., can not be regulated by the immune system, but the parenchyma stays healthy (chronic inflammation), 3) leads to a dysregulation of the immune and the parenchyma and hence a pathological state (eventually organ failure). Moreover, we have compared the probability for the emergence of pathological states depending on the age parameter obtained from the simulation of our model with empirical data for the hospitalization incidence of sepsis in Germany. This comparison shows a striking similarity that needs to be investigated in further studies, however, providing first evidence for the strength of our functional modeling approach.

This study lines up with other works in the emerging field of network physiology (Ivanov, 2021). Network physiology is a rather young interdisciplinary research area bridging between physiological modeling approaches from the micro to the macro scale. In the theory and application of dynamical systems, the network perspective has revolutionized (Newman, 2003) the field over the last 20 years, as it also allows for describing interaction structures on various spatial scales. Bringing together network science, dynamical system theory and physiological modeling, network physiology is a promising framework for getting insight into systemic diseases such as sepsis. Our approach provides a first step towards a functional dynamic modeling of sepsis. The basic limitation of our model is the approximation of the metabolic dynamics of the cells by simple one-dimensional phase oscillators, and of the cell interaction in the parenchyma by a fixed structural component and an adaptive component (cytokines) describing the information exchange between the cells. Further, the cellular interaction within the immune system is described only by adaptive information exchange, and the interaction between the immune system and the organs is described by a fixed coupling strength, and no differentiation between specific organs or different immune cells is made. An extension of our results guided by a systemic viewpoint, however, will pave the way for a deeper understanding of how the systemic spreading into other organ systems in case of sepsis occurs or how a relapse could be predicted. For this, one needs to further investigate which factors are crucial for a systemic spreading of disease, learn how different organ systems are interrelated, and how the complementary perspectives from physiology, network science, and dynamical systems can be further developed in an interdisciplinary context.

Data Availability Statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Author Contributions

All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.

Funding

This work was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation, project Nos. 429685422 and 440145547) and the Open Access Publication Fund of TU Berlin.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnetp.2022.904480/full#supplementary-material

References

Aoki, T., and Aoyagi, T. (2009). Co-Evolution of Phases and Connection Strengths in a Network of Phase Oscillators. Phys. Rev. Lett. 102, 034101. doi:10.1103/PhysRevLett.102.034101

PubMed Abstract | CrossRef Full Text | Google Scholar

Aoki, T., and Aoyagi, T. (2011). Self-Organized Network of Phase Oscillators Coupled by Activity-Dependent Interactions. Phys. Rev. E 84, 066109. doi:10.1103/PhysRevE.84.066109

PubMed Abstract | CrossRef Full Text | Google Scholar

Bartsch, R. P., Liu, K. K. L., Bashan, A., and Ivanov, P. C. (2015). Network Physiology: How Organ Systems Dynamically Interact. PLoS One 10, e0142143. doi:10.1371/journal.pone.0142143

PubMed Abstract | CrossRef Full Text | Google Scholar

Bartsch, R. P., Schumann, A. Y., Kantelhardt, J. W., Penzel, T., and Ivanov, P. C. (2012). Phase Transitions in Physiologic Coupling. Proc. Natl. Acad. Sci. U.S.A. 109, 10181–10186. doi:10.1073/pnas.1204568109

PubMed Abstract | CrossRef Full Text | Google Scholar

Bashan, A., Bartsch, R. P., Kantelhardt, J. W., Havlin, S., and Ivanov, P. C. (2012). Network Physiology Reveals Relations between Network Topology and Physiological Function. Nat. Commun. 3, 702. doi:10.1038/ncomms1705

PubMed Abstract | CrossRef Full Text | Google Scholar

Berner, R., Fialkowski, J., Kasatkin, D., Nekorkin, V., Yanchuk, S., and Schöll, E. (2019). Hierarchical Frequency Clusters in Adaptive Networks of Phase Oscillators. Chaos 29, 103134. doi:10.1063/1.5097835

PubMed Abstract | CrossRef Full Text | Google Scholar

Berner, R., Mehrmann, V., Schöll, E., and Yanchuk, S. (2021). The Multiplex Decomposition: An Analytic Framework for Multilayer Dynamical Networks. SIAM J. Appl. Dyn. Syst. 20, 1752–1772. doi:10.1137/21m1406180

CrossRef Full Text | Google Scholar

Berner, R., Polanska, A., Schöll, E., and Yanchuk, S. (2020). Solitary States in Adaptive Nonlocal Oscillator Networks. Eur. Phys. J. Spec. Top. 229, 2183–2203. doi:10.1140/epjst/e2020-900253-0

CrossRef Full Text | Google Scholar

Berner, R., Sawicki, J., and Schöll, E. (2020). Birth and Stabilization of Phase Clusters by Multiplexing of Adaptive Networks. Phys. Rev. Lett. 124, 088301. doi:10.1103/PhysRevLett.124.088301

PubMed Abstract | CrossRef Full Text | Google Scholar

Berner, R., Schöll, E., and Yanchuk, S. (2019). Multiclusters in Networks of Adaptively Coupled Phase Oscillators. SIAM J. Appl. Dyn. Syst. 18, 2227–2266. doi:10.1137/18m1210150

CrossRef Full Text | Google Scholar

Berner, R., Vock, S., Schöll, E., and Yanchuk, S. (2021). Desynchronization Transitions in Adaptive Networks. Phys. Rev. Lett. 126, 028301. doi:10.1103/physrevlett.126.028301

PubMed Abstract | CrossRef Full Text | Google Scholar

Berner, R., Yanchuk, S., Maistrenko, Y., and Schöll, E. (2021). Generalized Splay States in Phase Oscillator Networks. Chaos 31, 073128. doi:10.1063/5.0056664

PubMed Abstract | CrossRef Full Text | Google Scholar

Brunkhorst, F. M., Weigand, M. A., Weigand, M. A., Pletz, M., Gastmeier, P., Lemmen, S. W., et al. (2020). S3-Leitlinie Sepsis - Prävention, Diagnose, Therapie und Nachsorge. Med. Klin. Intensivmed. Notfmed 115, 37–109. doi:10.1007/s00063-020-00685-0

CrossRef Full Text | Google Scholar

Caporale, N., and Dan, Y. (2008). Spike Timing-Dependent Plasticity: A Hebbian Learning Rule. Annu. Rev. Neurosci. 31, 25–46. doi:10.1146/annurev.neuro.31.060407.125639

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, L., Liu, R., Liu, Z.-P., Li, M., and Aihara, K. (2012). Detecting Early-Warning Signals for Sudden Deterioration of Complex Diseases by Dynamical Network Biomarkers. Sci. Rep. 2, 342. doi:10.1038/srep00342

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Z., Hu, K., Stanley, H. E., Novak, V., and Ivanov, P. C. (2006). Cross-Correlation of Instantaneous Phase Increments in Pressure-Flow Fluctuations: Applications to Cerebral Autoregulation. Phys. Rev. E 73, 031915. doi:10.1103/PhysRevE.73.031915

PubMed Abstract | CrossRef Full Text | Google Scholar

Daido, H. (1994). Generic Scaling at the Onset of Macroscopic Mutual Entrainment in Limit-Cycle Oscillators with Uniform All-To-All Coupling. Phys. Rev. Lett. 73, 760–763. doi:10.1103/physrevlett.73.760

PubMed Abstract | CrossRef Full Text | Google Scholar

Delves, P. J., Martin, S. J., Burton, D. R., and Roitt, I. M. (2016). Roitt’s Essential Immunology. 13th ed. Chichester: Wiley-Blackwell.

Google Scholar

Desroches, M., Guckenheimer, J., Krauskopf, B., Kuehn, C., Osinga, H. M., and Wechselberger, M. (2012). Mixed-Mode Oscillations with Multiple Time Scales. SIAM Rev. 54, 211–288. doi:10.1137/100791233

CrossRef Full Text | Google Scholar

Egger, G. (2005). Die Akute Entzündung. Wien New York: Springer.

Google Scholar

Elisia, I., Lam, V., Hofs, E., Li, M. Y., Hay, M., Cho, B., et al. (2017). Effect of Age on Chronic Inflammation and Responsiveness to Bacterial and Viral Challenges. PLoS One 12, e0188881. doi:10.1371/journal.pone.0188881

PubMed Abstract | CrossRef Full Text | Google Scholar

Fleischmann, C., Thomas–Rueddel, D. O., Hartmann, M., Hartog, C. S., Welte, T., Heublein, S., et al. (2016). Hospital Incidence and Mortality Rates of Sepsis: an Analysis of Hospital Episode (DRG) Statistics in Germany from 2007 to 2013. Dtsch. Arztebl. Int. 113, 159. doi:10.3238/arztebl.2016.0159

PubMed Abstract | CrossRef Full Text | Google Scholar

Gerster, M., Berner, R., Sawicki, J., Zakharova, A., Škoch, A., Hlinka, J., et al. (2020). FitzHugh-Nagumo Oscillators on Complex Networks Mimic Epileptic-Seizure-Related Synchronization Phenomena. Chaos 30, 123130. doi:10.1063/5.0021420

PubMed Abstract | CrossRef Full Text | Google Scholar

Hotchkiss, R. S., Moldawer, L. L., Opal, S. M., Reinhart, K., Turnbull, I. R., and Vincent, J.-L. (2016). Sepsis and Septic Shock. Nat. Rev. Dis. Prim. 2, 16045. doi:10.1038/nrdp.2016.45

PubMed Abstract | CrossRef Full Text | Google Scholar

Ivanov, P. C., and Bartsch, R. P. (2014). “Network Physiology: Mapping Interactions between Networks of Physiologic Networks,” in Networks of Networks: The Last Frontier of Complexity. Editors G. D'Agostino,, and A. Scala (Cham: Springer), 203–222. chap. 10. doi:10.1007/978-3-319-03518-5_10

CrossRef Full Text | Google Scholar

Ivanov, P. C., Ma, Q. D. Y., and Bartsch, R. P. (2009). Maternal-Fetal Heartbeat Phase Synchronization. Proc. Natl. Acad. Sci. U.S.A. 106, 13641–13642. doi:10.1073/pnas.0906987106

PubMed Abstract | CrossRef Full Text | Google Scholar

Ivanov, P. C. (2021). The New Field of Network Physiology: Building the Human Physiolome. Front. Netw. Physiol. 1, 1. doi:10.3389/fnetp.2021.711778

CrossRef Full Text | Google Scholar

Kasatkin, D. V., Yanchuk, S., Schöll, E., and Nekorkin, V. I. (2017). Self-Organized Emergence of Multilayer Structure and Chimera States in Dynamical Networks with Adaptive Couplings. Phys. Rev. E 96, 062211. doi:10.1103/PhysRevE.96.062211

PubMed Abstract | CrossRef Full Text | Google Scholar

Kivelä, M., Arenas, A., Barthélemy, M., Gleeson, J. P., Moreno, Y., and Porter, M. A. (2014). Multilayer Networks. J. Complex Netw. 2, 203–271. doi:10.1093/comnet/cnu016

CrossRef Full Text | Google Scholar

Kuehn, C. (2015). Multiple Time Scale Dynamics. Cham: Springer.

Google Scholar

Kuramoto, Y. (1984). Chemical Oscillations, Waves and Turbulence. Berlin: Springer-Verlag.

Google Scholar

Leyva, I., Sendiña-Nadal, I., Sevilla-Escoboza, R., Vera-Avila, V. P., Chholak, P., and Boccaletti, S. (2018). Relay Synchronization in Multiplex Networks. Sci. Rep. 8, 8629. doi:10.1038/s41598-018-26945-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, A., Liu, K. K. L., Bartsch, R. P., and Ivanov, P. C. (2016). Delay-Correlation Landscape Reveals Characteristic Time Delays of Brain Rhythms and Heart Interactions. Phil. Trans. R. Soc. A 374, 20150182. doi:10.1098/rsta.2015.0182

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, R., Aihara, K., and Chen, L. (2013). Dynamical Network Biomarkers for Identifying Critical Transitions and Their Driving Networks of Biologic Processes. Quant. Biol. 1, 105–114. doi:10.1007/s40484-013-0008-0

CrossRef Full Text | Google Scholar

Liu, R., Li, M., Liu, Z.-P., Wu, J., Chen, L., and Aihara, K. (2012). Identifying Critical Transitions and Their Leading Biomolecular Networks in Complex Diseases. Sci. Rep. 2, 813. doi:10.1038/srep00813

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, R., Wang, X., Aihara, K., and Chen, L. (2013). Early Diagnosis of Complex Diseases by Molecular Biomarkers, Network Biomarkers, and Dynamical Betwork Biomarkers. Med. Res. Rev. 34 (3), 455–478. doi:10.1002/med.21293

PubMed Abstract | CrossRef Full Text | Google Scholar

Lücken, L., Popovych, O. V., Tass, P. A., and Yanchuk, S. (2016). Noise-Enhanced Coupling between Two Oscillators with Long-Term Plasticity. Phys. Rev. E 93, 032210. doi:10.1103/PhysRevE.93.032210

PubMed Abstract | CrossRef Full Text | Google Scholar

Madadi Asl, M., Valizadeh, A., and Tass, P. A. (2018). Dendritic and Axonal Propagation Delays May Shape Neuronal Networks with Plastic Synapses. Front. Physiol. 9, 1849. doi:10.3389/fphys.2018.01849

PubMed Abstract | CrossRef Full Text | Google Scholar

Maistrenko, Y. L., Lysyansky, B., Hauptmann, C., Burylko, O., and Tass, P. A. (2007). Multistability in the Kuramoto Model with Synaptic Plasticity. Phys. Rev. E 75, 066207. doi:10.1103/PhysRevE.75.066207

PubMed Abstract | CrossRef Full Text | Google Scholar

Moorman, J. R., Lake, D. E., and Ivanov, P. C. (2016). Early Detection of Sepsis-A Role for Network Physiology? Crit. Care Med. 44, e312–e313. doi:10.1097/CCM.0000000000001548

PubMed Abstract | CrossRef Full Text | Google Scholar

Mostel, Z., Perl, A., Marck, M., Mehdi, S. F., Lowell, B., Bathija, S., et al. (2019). Post-Sepsis Syndrome - an Evolving Entity that Afflicts Survivors of Sepsis. Mol. Med. 26, 6. doi:10.1186/s10020-019-0132-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Nekorkin, V. I., and Kasatkin, D. V. (2016). Dynamics of a Network of Phase Oscillators with Plastic Couplings. AIP Conf. Proc. 1738, 210010. doi:10.1063/1.4951993

CrossRef Full Text | Google Scholar

Newman, M. E. J. (2003). The Structure and Function of Complex Networks. SIAM Rev. 45, 167–256. doi:10.1137/s003614450342480

CrossRef Full Text | Google Scholar

Omelchenko, I., Omel’chenko, O. E., Hövel, P., and Schöll, E. (2013). When Nonlocal Coupling between Oscillators Becomes Stronger: Patched Synchrony or Multichimera States. Phys. Rev. Lett. 110, 224101. doi:10.1103/physrevlett.110.224101

PubMed Abstract | CrossRef Full Text | Google Scholar

Omel’chenko, O. E., and Knobloch, E. (2019). Chimerapedia: Coherence–Incoherence Patterns in One, Two and Three Dimensions. New J. Phys. 21, 093034. doi:10.1088/1367-2630/ab3f6b

CrossRef Full Text | Google Scholar

Omel’chenko, O. E. (2018). The Mathematics behind Chimera States. Nonlinearity 31, R121. doi:10.1088/1361-6544/aaaa07

CrossRef Full Text | Google Scholar

Omel’chenko, O. E., Wolfrum, M., and Maistrenko, Y. (2010). Chimera States as Chaotic Spatiotemporal Patterns. Phys. Rev. E 81, 065201(R). doi:10.1103/PhysRevE.81.065201

CrossRef Full Text | Google Scholar

Popovych, O. V., Yanchuk, S., and Tass, P. A. (2013). Self-Organized Noise Resistance of Oscillatory Neural Networks with Spike Timing-Dependent Plasticity. Sci. Rep. 3, 2926. doi:10.1038/srep02926

PubMed Abstract | CrossRef Full Text | Google Scholar

Prescott, H. C., Osterholzer, J. J., Langa, K. M., Angus, D. C., and Iwashyna, T. J. (2016). Late Mortality after Sepsis: Propensity Matched Cohort Study. BMJ 353, i2375. doi:10.1136/bmj.i2375

PubMed Abstract | CrossRef Full Text | Google Scholar

Röhr, V., Berner, R., Lameu, E. L., Popovych, O. V., and Yanchuk, S. (2019). Frequency Cluster Formation and Slow Oscillations in Neural Populations with Plasticity. PLoS One 14, e0225094. doi:10.1371/journal.pone.0225094

PubMed Abstract | CrossRef Full Text | Google Scholar

Rosendal, C., Bouchon, A. I., Bierhaus, A., Nawroth, P. P., Martin, E., Bardenheuer, H. J., et al. (2004). Bedeutung der angeborenen immunantwort in der sepsis. Der Anaesthesist 53, 10–28. doi:10.1007/s00101-003-0626-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Sakaguchi, H., and Kuramoto, Y. (1986). A Soluble Active Rotater Model Showing Phase Transitions via Mutual Entertainment. Prog. Theor. Phys. 76, 576–581. doi:10.1143/ptp.76.576

CrossRef Full Text | Google Scholar

Sawicki, J., Berner, R., Löser, T., and Schöll, E. (2022). Modeling Tumor Disease and Sepsis by Networks of Adaptively Coupled Phase Oscillators. Front. Netw. Physiol. 1, 730385. doi:10.3389/fnetp.2021.730385

CrossRef Full Text | Google Scholar

Sawicki, J. (2019). Delay Controlled Partial Synchronization in Complex Networks. Heidelberg: Springer. Springer Theses.

Google Scholar

Sawicki, J., Omelchenko, I., Zakharova, A., and Schöll, E. (2018). Delay Controls Chimera Relay Synchronization in Multiplex Networks. Phys. Rev. E 98, 062224. doi:10.1103/physreve.98.062224

CrossRef Full Text | Google Scholar

Schmidt, K., Gensichen, J., Fleischmann-Struzek, C., Bahr, V., Pausch, C., Sakr, Y., et al. (2020). Long-Term Survival Following Sepsis: Results of a Single-Center Registry Study with 4-Year Follow-Up. Dtsch. Arztebl. Int. 117, 775. doi:10.3238/arztebl.2020.0775

PubMed Abstract | CrossRef Full Text | Google Scholar

Schöll, E. (2020). Chimeras in Physics and Biology: Synchronization and Desynchronization of Rhythms. Nova Acta Leopoldina 425, 67. doi:10.26164/leopoldina_10_00275

CrossRef Full Text | Google Scholar

Schöll, E. (2021). Partial Synchronization Patterns in Brain Networks. Europhys. Lett. 136, 18001. doi:10.1209/0295-5075/ac3b97

CrossRef Full Text | Google Scholar

Schulte, W., Bernhagen, J., and Bucala, R. (2013). Cytokines in Sepsis: Potent Immunoregulators and Potential Therapeutic Targets-An Updated View. Mediat. Inflamm. 2013, 165974. doi:10.1155/2013/165974

PubMed Abstract | CrossRef Full Text | Google Scholar

Shi, J., Kirihara, K., Tada, M., Fujioka, M., Usui, K., Koshiyama, D., et al. (2022). Criticality in the Healthy Brain. Front. Netw. Physiol. 1, 755685. doi:10.3389/fnetp.2021.755685

CrossRef Full Text | Google Scholar

Singer, M., Deutschman, C. S., Seymour, C. W., Shankar-Hari, M., Annane, D., Bauer, M., et al. (2016). The Third International Consensus Definitions for Sepsis and Septic Shock (Sepsis-3). JAMA 315, 801. doi:10.1001/jama.2016.0287

PubMed Abstract | CrossRef Full Text | Google Scholar

Storch, V., Welsch, U., and Wink, M. (2013). Evolutionsbiologie. 3rd ed. Berlin, Heidelberg: Springer Spektrum.

Google Scholar

Thomas, L. (2020). Clinical Laboratory Diagnostics. Available at: https://www.clinical-laboratory-diagnostics-2020.com/.

Google Scholar

Tisoncik, J. R., Korth, M. J., Simmons, C. P., Farrar, J., Martin, T. R., and Katze, M. G. (2012). Into the Eye of the Cytokine Storm. Microbiol. Mol. Biol. Rev. 76, 16–32. doi:10.1128/mmbr.05015-11

PubMed Abstract | CrossRef Full Text | Google Scholar

Vineis, P., Schatzkin, A., and Potter, J. D. (2010). Models of Carcinogenesis: an Overview. Carcinogenesis 31, 1703–1709. doi:10.1093/carcin/bgq087

PubMed Abstract | CrossRef Full Text | Google Scholar

Weis, S., Pletz, M. W., and Bauer, M. (2017). “Sepsis – Neue Definition, Neue Kontroversen,” in Epidemiologisches Bulletin (Robert Koch-Institut, Epidemiologie und Gesundheitsberichterstattung), 37. Chichester: Wiley-Blackwell.

Google Scholar

Xia, S., Zhang, X., Zheng, S., Khanabdali, R., Kalionis, B., Wu, J., et al. (2016). An Update on Inflamm-Aging: Mechanisms, Prevention, and Treatment. J. Immunol. Res. 2016, 8426874. doi:10.1155/2016/8426874

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, L., Chen, Z., Hu, K., Stanley, H. E., and Ivanov, P. C. (2006). Spurious Detection of Phase Synchronization in Coupled Nonlinear Oscillators. Phys. Rev. E 73, 065201. doi:10.1103/PhysRevE.73.065201

PubMed Abstract | CrossRef Full Text | Google Scholar