Modeling Tumor Disease and Sepsis by Networks of Adaptively Coupled Phase Oscillators

In this study, we provide a dynamical systems perspective to the modelling of pathological states induced by tumors or infection. A unified disease model is established using the innate immune system as the reference point. We propose a two-layer network model for carcinogenesis and sepsis based upon the interaction of parenchymal cells and immune cells via cytokines, and the co-evolutionary dynamics of parenchymal, immune cells, and cytokines. Our aim is to show that the complex cellular cooperation between parenchyma and stroma (immune layer) in the physiological and pathological case can be qualitatively and functionally described by a simple paradigmatic model of phase oscillators. By this, we explain carcinogenesis, tumor progression, and sepsis by destabilization of the healthy homeostatic state (frequency synchronized), and emergence of a pathological state (desynchronized or multifrequency cluster). The coupled dynamics of parenchymal cells (metabolism) and nonspecific immune cells (reaction of innate immune system) are represented by nodes of a duplex layer. The cytokine interaction is modeled by adaptive coupling weights between the nodes representing the immune cells (with fast adaptation time scale) and the parenchymal cells (slow adaptation time scale) and between the pairs of parenchymal and immune cells in the duplex network (fixed bidirectional coupling). Thereby, carcinogenesis, organ dysfunction in sepsis, and recurrence risk can be described in a correct functional context.


INTRODUCTION
Tumors and sepsis are diseases of different genesis. They have very different time scales and the therapies are completely different. But the age incidence, the risk factors, the outcome, and the temporal recurrence behavior are similar. This justifies the attempt to describe tumor disease and sepsis with a unified disease model. A shift of the paradigm is proposed by choosing the nonspecific innate immune system as the reference point for both diseases Löser (2020). The innate immune system interacts with mutant cells and pathogens. Its role is to maintain the integrity of the organism by actively eliminating foreign organisms, degrading the organism's own damaged cells, and activating and coordinating wound healing. It can ward off about 99% of all infections Bomans et al. (2018); it arose with the beginning of multicellular life, and has grown and protected life to this day Male et al. (2012); Rich and Chaplin (2019). It is ubiquitously present in the organism, is responsive without significant dead time, and is based on dynamically balanced activator-inhibitor mechanisms. Its regulation is essentially decentrally organized. The innate immune system includes humoral and cellular components, the endothelium, and the tissue stroma. With the multitude of components interacting in the innate immune system, it is organized in a complex way and has a broad response spectrum. The innate immune system is not completely resistant to disturbances; it can be brought to dysregulation. Tumor disease and sepsis are two prominent examples of this behavior. Although tumor disease and sepsis correspond to different medical conditions, after the initial stages a relatively uniform course develops Weinberg (2014), Singer et al. (2016). The variability of tumor disease results from the initial genetic state of the tumor cells, their subsequent mutations, epithelialmesenchymal transition, and interactions with the innate immune system Longo (2011). Recently, in case of cancer the crosstalk of tumor cells with immune cells has been investigated in more detail Zhang et al. (2021). In sepsis, the innate immune system is activated by infection and the clinical course is determined by the individual patient's initial condition. In both cases, the innate immune system can be brought to dysregulation and develop its own clinical picture such as cachexia, coagulation disorders, and organ failure. In the case of tumor disease, this is compounded by space-occupying lesions, tissue invasion and destruction through the nonphysiological production of proliferation factors, cytokines, and chemokines. Possibly, not only the cytokine concentration, the cytokine mix but also the gradient of their increase is responsible for the extent of dysregulation of the innate immune system and thus for the disease consequences Morán et al. (2013), Altan-Bonnet and Mukherjee (2019).
The 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 of carcogenesis, as has been reviewed elsewhere Vineis et al. (2010), but it rather describes tumor and sepsis and recurrence risk in a dynamical functional context. Complex networks are an ubiquitous paradigm in nature and technology, with a wide field of applications ranging from physics, chemistry, biology, neuroscience, to engineering, and socio-economic systems. Of particular interest are adaptive networks, where the connectivity changes in time, for instance, in chemical or biochemical systems Jain and Krishna (2001), where the reaction rates adapt dynamically depending on the variables of the system, or in neuronal synaptic plasticity Markram et al. (1997), Abbott and Nelson (2000), Meisel and Gross (2009), Lücken et al. (2016), in epidemics Gross et al. (2006), and in biological or social systems Gross and Blasius (2008). Another focus of recent research in network science are multilayer networks, which are systems interconnected through different types of links De Domenico et al. (2013), Boccaletti et al. (2014), Kivelä et al. (2014), De Domenico et al. (2015. A special case of multilayer networks are multiplex topologies, where each layer contains the same set of nodes, and only pairwise connections between corresponding nodes from neighboring layers exist Zhang et al. (2015), Maksimenko et al. (2016), Andrzejak et al. (2017), Leyva et al. (2017), Sawicki et al. (2018), Nikitin et al. (2019), Omelchenko et al. (2019), Rybalova et al. (2019), Drauschke et al. (2020), Berner et al. (2021a), Sawicki et al. (2021), Shepelev et al. (2021).
In this article, we propose a two-layer network model for carcinogenesis and sepsis based upon the interaction of parenchymal cells and immune cells via cytokines and the coevolutionary dynamics of parenchymal, immune cells, and cytokines. Parenchyma is the bulk of functional substance in an organ, in contrast to the stroma, which refers to the structural tissue of organs or structures. In many organs the parenchyma consists of epithelial cells. We stress that our model is not a detailed model of organs or of biochemical processes but a functional model of dynamic interactions. Thus, cytokines are not modeled in terms of concentrations but rather in terms of the cytokine-mediated information flow between nodes within each layer and between the parenchymal and the immune layer, describing the cytokine activity. In the following, we refer to the stroma as the immune layer. The article is organized as follows: In Section 2, we provide a brief overview of the physiology of tumor disease and sepsis. The functional network model is introduced in Section 3. Here, we discuss all variables, parameters, and their physiological meaning. Moreover, the methods and measures used for the subsequent numerical analysis are introduced and highlighted in their physiological context. In Section 4, we discuss various dynamical scenarios of the model simulations of tumor disease that are observed in the presence of pathological cells. Section 5 presents analogous computer simulations for sepsis. The results are summarized in Section 6.

Initial Situation
The function of the innate immune system is ensured in a siteindependent manner by the coordinated activity of humoral and cellular components. Coordination also occurs via cytokines. With new cytokine sources building up over individual lifetimes, the stability of regulatory behavior in the innate immune system changes. Contributing systemic factors include inflammaging Campisi (2014), Calder et al. (2017), activation of the coagulation system Tragl (1999), Franceschi and Campisi (2014), increase of body fat Fasshauer et al. (2004), Gaillard (2007), Chovatiya and Medzhitov (2014), lack of exercise Elisia et al. (2017), Fulop et al. (2018), age-related normal fibrosis Beneke (1971), Nemetschek (1971), chronic inflammation as local factors Virchow (1978), or smoking via initiation of local hypoxia Brunkhorst et al. (2018). The innate immune system exhibits pro-inflammatory activation. The number of individual variables suggest a wide range of activation levels, which still depend on cytokine gene polymorphisms, among other factors. Physiological cytokine production and additional pro-inflammatory cytokine activity (pattern, concentration, emitters involved, gradient of increase) determine the activation level of the innate immune system.

Tumor Disease
In tumor disease, mutant cells interact with macrophages of the innate immune system localized in the stroma. Mutated cells are genetically altered and genetically unstable parenchymal cells. They are a new cell entity with altered regulatory behavior. In interaction with the normal stroma (immune layer), e.g. in young healthy individuals, mutant cells do not find survival conditions. They persist silently or trigger apoptosis via tissue surveillance Löser (2018). Mutant cells in an area pre-damaged by chronic inflammation receive cytokines from the activated stroma, undergo an epithelial-mesenchymal transition Lamouille et al. (2014), Chockley and Keshamouni (2016), can proliferate, disseminate, and reprogram normal macrophages into tumorassociated macrophages Wu and Zhou (2009), Karlsson et al. (2017), Mantovani et al. (2017). They help to ensure oxygen and nutrient supply to the tumor. Tumor cells, tumor-associated macrophages, and normal stromal cells communicate via cytokines. The communication strength depends on tumor cell genetics and stroma activation. It increases with the tumor cell mass and the mass of specific tumor stroma. A mutually activating circular process starts between the tumor cells and the innate immune system. A specific cytokine microclimate is formed around the developing tumor, the size of which grows with the increase in tumor mass. Cytokines enter the lymph nodes with the lymph and enable the tumor cells that have floated there to proliferate. Lymphogenic metastasis begins, continuing downstream in additional lymph nodes according to the same pattern. Organ-specific metastatic patterns are formed Walther (1948). The circulating cytokine concentration increases and triggers the general symptoms of tumor disease such as cachexia Arends et al. (2015) and finally the lethal coagulopathy, lung failure, organ failure or lethal spaceoccupying lesions and tissue destruction. The postulate of Lewis Thomas Thomas (1972) also applies to tumor disease, according to which the interactions with the organism triggered by the tumor cells are responsible for the disease, i.e., tumor growth, enabling malignant cells to become invasive and destructive, recurrence, and finally induce the lethal general symptoms. Hematogenous metastasis may start when the disseminated tumor cells are adequately supplied with growth factors (cytokines) via further activation of the innate immune system.

Sepsis
The initial conditions for the transition of an infection into sepsis are analogous to the survival conditions of mutated cells in the tissue, a pre-activated innate immune system. Sepsis is triggered when the innate immune system no longer succeeds in locally fixing invading pathogens, but allows them to enter the organism. There they interact with the already pre-activated innate immune system. The immune system is activated systemically. Antibiosis or surgical sanitation can reduce or eliminate the pathogen load. Depending on the type of reaction, the pro-inflammatory response is stopped or it continues to escalate in varying degrees Seymour et al. (2019), for which also cytokine gene polymorphisms are responsible Majetschak and Schade (2001), Hotchkiss et al. (2016), Thomas (2020a), Thomas (2020b). Proinflammatory cytokines act on endothelial cells and the coagulation system. Microthrombi amplify inflammation through secretion of growth factors and cytokines Hotchkiss et al. (2016). The endothelial permeabilization barrier is opened, fluid retention in tissues occurs, blood pressure drops, and must be stabilized with fluid administration Brunkhorst et al. (2018). Prolonged oxygen diffusion pathways and the resulting tissue hypoxia trigger cytokine release causing a "second hit", resulting in metabolic changes in the parenchyma Bomans et al. (2018). Organ failure, particularly of the kidneys, lungs, and liver is imminent. Endogenous molecules with immunogenic effects, especially mitochondrial DNA Franceschi et al. (2018), released e.g. after trauma, burns, pancreatitis or surgery, can have the same effect as exogenous pathogens.

Relapse
Tumors and sepsis have a recurrence behavior, the frequency of which correlates with the stage of the primary disease. In both diseases, the activated innate immune system is an initial prerequisite for the primary disease, which is brought to dysregulation in the course of the disease. After tumor removal and elimination of infection, the activation status reduces but remains at least at the level before disease onset. Thus, cytokines are produced for remaining tumor cells as a result of pro-inflammatory activation Gastpar (1982), Eichinger and Lechner (2004), Arends et al. (2015), Lippman (2016), which promote their proliferation and allow the organization of the tumor stroma. The process kinetics is determined by the proliferation rate of tumor cells, the availability of cytokines, and by the ability of tumor cells to organize their own stroma with connection to the blood supply. After clearance of the systemic infection, the innate immune system has an activation status at least equivalent to that before septic shock. The infection itself appears to have triggered irreversible elements of trained immunity Bomans et al. (2018). Surviving patients after septic shock thus die in 50% within the first two post-sepsis years from persistent and secondary nosocomial infections, tumors, or cardiac failure Prescott et al. (2016).

Schematic Model for Tumor Disease and Sepsis
Organic tissue consisting of parenchymal cells and immune cells is shown schematically in Figure 1. We depict the initial and final configurations for tumor disease and sepsis on a tissue element. The tissue element consists of the epithelial parenchyma, the basal membrane, and the stroma. The parenchyma is the organspecific functional layer. The basal membrane separates the parenchyma and stroma and is made of collagen of type IV. 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).
The initial process of carcinogenesis is shown in Figures 1A-C. In Figure 1A, the parenchyma is normal, the stroma is inflammatorily activated. In Figure 1B, a mutation of a tumor cell has occurred in the parenchyma, the inflammatory activation of the stroma continues. Finally in Figure 1C, the mutant cells proliferate, break through the basal membrane and migrate into the blood and lymphatic vessels, and cytokines are emitted. Figures 1D-F shows the systemic effects of the final process of tumor disease and sepsis: Figure 1D depicts a tissue element that is not directly affected by tumor, metastasis, or primary inflammation. The stroma exhibits inflammatory activation. In Figure 1E, a systemic cytokine storm is occurring. Cytokine production by the stromal cells is additionally stimulated. The parenchyma responds to the primary and secondary cytokine storm by uncoupling mitochondrial respiration and switching to aerobic glycolysis. The energy supply is no longer sufficient for organ-specific cellular functions and the organ fails. Figure 1F shows a variant of Figure 1E. The cytokine storm activates the endothelium, blood coagulation is activated, oxygen transport breaks down, and hypoxia arises. The stromal cells respond with cytokine release. Both together, the hypoxia and the cytokine Frontiers in Network Physiology | www.frontiersin.org January 2022 | Volume 1 | Article 730385 release, lead to the collapse of organ-specific cell functions in the parenchyma and the organ fails.

Dynamical Two-Layer Network Model
The unified disease model with the reference point given by the innate immune system is the basis for our model, which includes disease-specific initial conditions, mutant cells for tumor disease, and infection-driven cytokine dysregulation. A volume element of tissue consisting of parenchyma, basal membrane and stroma is used as a model for tumor disease and sepsis, describing the functional interactions between parenchyma (organ tissue) and stroma (immune layer). We represent the network layer of parenchymal cells (superscript 1) by N phase oscillators ϕ 1 i , i 1, . . . , N, with partly fixed and partly adaptive coupling weights Junqueira et al. (1995), and the network layer of immune cells (superscript 2) by N adaptively coupled phase oscillators ϕ 2 i . The communication through cytokines which mediate the interaction between the parenchymal cells is modeled by the coupling weights κ 1 ij , and those between the immune cells by coupling weights κ 2 ij . 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. As discussed in detail in Section 3.3, healthy cells and tumor cells differ by their metabolic activity, i.e., tumor cells are less energy-efficient and thus have a faster cellular metabolism, which is reflected in our phase oscillator model by a higher frequency. Thus in the healthy homeostatic equilibrium state all parenchymal cells have the same lower frequency, while the pathological state splits into two clusters with different frequencies, healthy, and unhealthy, i.e., a multifrequency cluster. The healthy state is assumed to be characterized by regular periodic, fully synchronized dynamics of the phase oscillators, i.e., all cells show the same collective frequency of cellular metabolism. The pathological state is described by multifrequency clusters with different frequencies, i.e., the pathological cells in the parenchyma attain a higher frequency, while the healthy cells are still frequencysynchronized with the "healthy" frequency. This loss of synchrony reflects the fact that in case of tumor disease, the malignant mutation basically leads to a loss of performance of the parenchymal cells. As a consequence, they are no longer fully coordinated, which leads to a loss of proliferation and apoptosis control Longo (2011), Weinberg (2014. Such a pathological condition leads also to an alteration of the metabolic activity of the immune cells Coussens and Werb (2002), Heerboth et al. (2015), Chockley and Keshamouni (2016), Porporato (2016), Razak et al. (2018), Greten and Grivennikov (2019). In case of sepsis, i.e., a systemic inflammation, cells belonging to the innate immune system produce cytokines in an unregulated way affecting the parenchyma or microcirculation leading to organ failure Singer et al. (2016), Matsumoto et al. (2018) which is also associated in our model with the loss of synchrony.
A general multiplex network with two layers each consisting of N identical adaptively coupled phase oscillators is described by where ϕ μ i ∈ [0, 2π) represents the phase of the ith oscillator (i 1, . . . , N) in the μth layer (μ 1, 2), ω 1 i ≡ ω i are the natural oscillator frequencies of the parenchymal cells which are distributed according to a probability distribution ρ( where r is the fraction of pathological parenchymal cells relative to the number of all parenchymal cells N, δ is the Dirac delta function, and ω p and ω h are the natural frequencies of pathological and healthy parenchymal cells, respectively. The value of ω 2 ≡ ω is the natural frequency of the immune cells. The interaction between the oscillators within each layer is determined by the intralayer connectivity weights a 1 ij ∈ [0, 1] (fixed interaction within an organ) and κ μ ij ∈ [−1, 1] (adaptive interaction mediated by cytokines), where the parenchymal layer has both fixed and adaptive couplings, while the immune layer has only adaptive coupling. Between the layers the interlayer coupling weights σ ≥ 0 are fixed and symmetric for both directions of interaction. The parameter α is a phase lag of the interaction modeling a time-delay Sakaguchi andKuramoto (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 ϵ 1 is assumed to be slow compared to the adaptation rate of the immune layer ϵ 2 , i.e., ϵ 1 ≪ ϵ 2 to account for the faster reaction of the immune cells Morán et al. (2013), Altan-Bonnet and Mukherjee (2019). Thus we have two classes of adaptive coupling weights modeling two different cytokine mechanisms on two different timescales. As a consequence of choosing two significantly different values for ϵ 1 and ϵ 2 , we obtain a system with multiple times scale dynamics, i.e., "slow-fast-faster" dynamics (ϵ 1 ≪ ϵ 2 ≪ 1) Kuehn (2015).
The phase lag parameter β of the adaptation function sin(ϕ μ i − ϕ μ j − β), also called plasticity rule in the neuroscience terminology Aoki and Aoyagi (2009), describes different adaptation rules that may occur. For instance, for β ( +)π/2, a symmetric rule Hoppensteadt and Izhikevich (1996), Seliger et al. (2002), Aoki (2015), Röhr et al. (2019) is obtained where the coupling κ ij decreases (increases) between any two systems with close-by phases. If β 0, the link κ ij will be strengthened if the ith oscillator is advancing the jth. Such a causal relationship is typical for spike-timing dependent plasticity in neuroscience Maistrenko et al. (2007), Caporale and Dan (2008), Popovych et al. (2013), Lücken et al. (2016). The matrix elements a 1 ij ∈ {0, 1} of the adjacency matrix A in the parenchymal layer are chosen as a 1 ij 1 if i ≠ j (global coupling). For normalization of the coupling terms the coupling sums in Eqs. 1, 2 are multiplied by the normalization factor 1 N .

Methodology and Measures
Heterogeneous dynamics, e.g., multifrequency clusters and tumor growth, may arise through the dynamic interaction of parenchymal cells, immune cells, and localized cytokine activity in a self-organized self-adaptive manner, even if the system parameters in the layers are chosen uniformly, i.e., homogeneous Berner et al. (2019a), Berner et al. (2020a). Heterogeneity can enter the parameters of oscillator networks in various ways. Very prominent are heterogeneities in the natural frequencies of the oscillators or in the connectivity structure Acebrón et al. (2005). Heterogeneous frequencies ω i are used to model pathological parenchymal cells which occur by spontaneous random mutations. We assume that all healthy cells possess the same natural frequency ω i ω h , and a fraction r of randomly chosen pathological cells possess the natural frequency ω i ω p . The characterization of tumor cells via their metabolic properties Longo (2011), Weinberg (2014) has revealed a difference in the metabolic activity between tumor and healthy cells Warburg et al. (1924). Warburg demonstrated that cells after a malignant mutation obtain their energy via aerobic glycolysis, which provides only 4 mol ATP/mol glucose, in contrast to normal cells whose metabolism is based upon mitochondrial breathing yielding 36 mol ATP/mol glucose. Thus tumor cells have a metabolic efficiency of only about 10% of healthy cells, and therefore they need more glucose. Other general differences in metabolic activity between tumor cells and healthy cells include proliferation rate, apoptosis rate, initiation of angiogenesis, tissue invasion, or metastatic capacity. Here, quantitative differences exist between different organ tumors. Ultimately, any change in metabolic performance, regardless of which metabolic branch is affected (structural, energetic, cell division, or special metabolism), can be conceived as a frequency change.
Using numerical simulations, we study whether the healthy state (frequency synchronized) is persistent also in the presence of a few pathological cells. Under certain conditions depending on various parameters (age, inflammaging, chronic inflammation, other basic diseases, obesity, smoking, lack of exercise), an unregulated cytokine expression and hence a possibly lethal tumor can occur. In these cases, the healthy (synchronized) state is not persistent anymore against the perturbation by the heterogeneity (tumor cells). For our study, the cytokine adaptivity parameter (which we call age parameter) β and the fraction r are considered as the main model parameters to account for the various system conditions. In the case of sepsis, we introduce a fixed initial perturbation of the cytokine activity in the immune layer representing a systemic immune response, while keeping the natural frequencies uniform. Similar to the tumor disease, we study the effect of this initial system perturbation on the emergence of the healthy state in dependence of the age parameter β. For our simulations we have used a Runge-Kutta method of order 4 with a fixed stepsize of Δt 0.05 and simulation time of T s 2000 time units where we have discarded the first 1,000 time units to account for transient dynamics.
In order to quantitatively characterize the dynamic collective state of the two-layer network, we introduce several measures. First, we use the mean phase velocities of the oscillators j in both layers μ 1, 2 with averaging time window T, and the spatially averaged mean phase velocity for each layer ω μ 1 N N j 1 〈 _ ϕ μ j 〉. Furthermore, for both layers μ 1, 2 we calculate the ensemble average s μ (ensemble size N E with ensemble elements E) of the standard deviation σ χ ( and the ensemble average of the corresponding normalized standard deviation . If the latter quantities are non-zero, they indicate the formation of multifrequency clusters, where the respective layer splits into clusters with different frequencies, which is indicative of a pathological state. In systems of adaptively coupled phase oscillators one often encounters frequency-synchronized dynamics. If the frequency is the same, the phases may still exhibit different behavior. 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 relative 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. (2020b), Berner et al. (2021c). In our set-up a splay state is interpreted as a more vulnerable collective state (synchronization of parenchyma and immune layer) where small perturbations can quickly lead to partial or complete desynchronization. Table 1 gives an overview of 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. In more detail, the quantities of the mathematical model which are listed in Table 1 have the following physiological correspondence: 1) Variable ϕ i mimics the metabolic activity of a parenchymal or immune cell i as a universal oscillatory phase variable.

Physiological Interpretation of Model Parameters
Here our aim is not to describe physiological processes on a detailed biochemical level, e.g., C-reactive protein (CRP) production after inflammatory activation, etc. 2) Variable κ ij describes the cytokine activity mediating information flow between cells j and i. It stands for the susceptibility of communications between cells via cytokines. Cytokines have pleiotropic properties, i.e., they can exert different, even contrary, effects depending on the initial situation and in combination with other cytokines. Therefore we model them by an adaptive coupling strength.
Reaction strength of cytokines is due to individually different gene polymorphisms. 3) Parameter α models the phase lag of the intralayer cell-to-cell coupling. Its origin is the time delay from the activation of a defined metabolic branch until the mediator release; within the inflammatory cascade, characteristic time delays exist for all reaction sub-steps starting from the triggering of the inflammatory reaction, for example, the reaction time of the individual liver cells starting from the input signal until the maximum blood concentration of interleukin IL-6 (2 h) and CRP (12 h) is reached. 4) Parameter β governs the adaptivity (or plasticity) rule of the cytokines. It mimicks a systemic sum parameter which may account for different influences such as physiological changes due to age in the extracellular matrix, inflammaging, systemic and local inflammatory baseline, adiposity, pre-existing illness, physical inactivity, nutritional influence, and others. In case of tumor disease, it can include the malignancy grade of tumor cells. For the sake of brevity, we call this parameter the age parameter. 5) Parameter ω denotes the natural frequency of the basic metabolic activity of a single cell. For instance, it can stand for the minor CRP production in a normal healthy state. 6) Parameters ϵ 1 and ϵ 2 denote the inverse relaxation times (half-life) of the cytokines in the parenchymal and immune layer, respectively. 7) Parameter r denotes the fraction of mutant (pathological) parenchymal cells in the tissue element, and is implemented in the model as the fraction of cells with a deviating (pathological) natural frequency ω p . 8) Parameter a ij denotes the genetically fixed intercellular communication pathways between parenchymal cells, they are wired by fixed cell-to-cell connections. In contrast, communication of parenchymal cells via cytokines and of cells within the immune layer runs in an open communication channel, which is controlled selfadaptively, after Shannon (1948). 9) Parameter σ denotes the interlayer coupling strength between the parenchymal and immune cells. It can be due to a mass transfer, e.g., cytokine expression of macrophages in the immune layer, and signal transfer into the parenchyma: macrophages invade the extracellular matrix and penetrate the basal membrane to the parenchyma. 10) The mean phase velocity 〈 _ ϕ i 〉 is a measure which describes the collective frequency of the cellular metabolism of cell i due to all interactions with other cells in both layers. It denotes the system performance of a single metabolic branch of the parenchymal or immune cells of a tissue element according to their activation state (normal healthy state or proinflammatory activation with increased CRP production). 11) The ensemble-averaged standard deviation of the mean phase velocities s 1 , s 2 is a measure which characterizes the inhomogeneity of metabolic activity within the parenchyma and the stroma (immune layer), respectively. It assumes nonzero values if the respective layer is not frequencysynchronized and splits into multifrequency clusters, e.g., a healthy cluster with one frequency and a pathological cluster with another frequency, and thus is a measure of pathogenicity in case of the parenchymal layer, or activation in case of the immune layer.

TUMOR DISEASE
In this section, we present exemplary computer simulations of our model to demonstrate different dynamic scenarios which this model can produce already in its simplest form. The model has not been refined or optimized with respect to the parameters, but our concern here is to display simulations which can describe principally different evolutions and outcomes of tumor disease. We assume N 200, ϵ 1 0.03, ϵ 2 0.3, a 1 ij 1 (global coupling in the parenchymal layer), σ 0.3 (interlayer coupling), natural frequencies ω h ω 2 0 (co-rotating frame of healthy parenchymal cells), ω p 1 (pathological cells describing spontaneous mutation, chosen randomly for R oscillators where r R/N), and the coupling delay parameter α is chosen the same for both layers such that the frequencies of the healthy parenchymal and immune layers are similar. The parameter β is used to model the influence of age, inflammation status, environment etc. and also the degree of malignancy of the tumor. Our system exhibits various dynamic patterns resembling healthy and pathological states. The results are depicted in Figure 2, where the mean phase velocities 〈 _ ϕ μ j 〉 are plotted in the right panel, and the absolute differences of the mean phase velocities with respect to a reference oscillator in the parenchyma 〈 _ ϕ 1 100 〉 and in the immune layer 〈 _ ϕ 2 300 〉, respectively, are plotted in the left panel. Figure 2A shows a healthy state with a low age parameter β 0.45π and only two tumor cells (r 1%) in the parenchymal layer. Although the two tumor cells in the parenchyma have different frequencies, the frequencysynchronized state with the healthy cells can be maintained due to the coupling. The coupled cooperative dynamics leads to a completely in-phase synchronized state in both layers. Note that the collective frequencies in the parenchymal layer and the immune layer are different (right panel) since the coupling terms are different. The same synchronized state of the parenchyma holds in Figure 2B for increased age parameter β 0.55π and increased number of tumor cells (r 4%), but first indications of pathological behavior are visible in the changed activity of the immune layer, showing lower frequencies of the associated immune cells coupled to the tumor cells (right panel, j 393, . . . , 400). For higher age parameter β 0.60π and the same number of tumor cells (r 4%), Figure 2C shows a scenario of an emergent pathological state, i.e., a two-frequency cluster with higher frequency of the pathological cells (red), in the parenchymal layer. A corresponding cluster of lower frequencies is observed in the immune layer ( Figure 2C) as well. With increasing age parameter β 0.66π and tumor size r 7.5%, Figure 2D shows an even more pathological state, where the frequency of the split-off pathological cluster (red) in the parenchyma is distinctly more different from the healthy cluster (blue), while the immune layer is no longer supporting a twocluster state but has adjusted in frequency to the parenchyma (right panel). This indicates an essential change in the dynamic state which is frequency synchronized close to zero, which indicates the absence of strong coupling contributions to the frequency. This frequency synchronization between parenchymal and immune layer is maintained if the large age parameter β 0.66π is kept, but the number of tumor cells is chosen as zero ( Figure 2E), and a healthy completely frequency-synchronized state is obtained. In the following, we will characterize this new dynamic state in more detail. Figure 3 shows details of these scenario for the same values of age parameter β and tumor size r. The left and right columns show snapshots of the cytokine matrices κ 1 ij (parenchymal layer) and κ 2 ij (immune layer), respectively. The second column shows snapshots of the instantaneous phases ϕ μ j , and the third column depicts space-time plots of the phases ϕ μ j (t) visualizing the oscillations. In the healthy state in Figures 3A,B the cytokine matrices are uniform and temporally constant within each cluster. The snapshots of the phases ϕ μ j (second column) and the space-time plots ϕ μ j (t) (third column) show homogeneous inphase oscillations in the parenchyma and the immune layer, respectively, though with different collective frequencies. In the pathological state in Figure 3C the small pathological cluster breaks off at j 193, . . . , 200. For even larger age parameter β 0.66π and tumor size r 7.5% in Figure 3D the phase dynamics changes qualitatively. We obtain a pathological two-frequency  cluster with strong frequency difference, where the healthy part is no longer in-phase synchronized, but becomes a splay state. With increasing age parameter β the frequency difference of the healthy cluster and the tumor cluster becomes larger, i.e., the tumor cells become autonomous. Thus larger age parameter can be associated with less favorable conditions of tumor disease (higher age, higher tumor malignity). A different scenario in dependence on the tumor size r with increased age parameter β leads from a mixed 2-frequency cluster to a healthy state, where the healthy cluster is a splay state ( Figure 3E). This state is different from the frequency clusters in Figures 3A-C, where each frequency cluster is phase synchronized. In Figures 3D,E the frequency of the healthy cluster in the parenchymal layer is fixed equal to that of the immune layer, but the phases of the parenchyma can no longer be kept in-phase, thus a strong overall perturbation of the healthy cluster is visible. The frequency of the pathological cluster is much more strongly separated from that of the healthy cluster. When the tumor size r is reduced for low age parameter β, below a certain r the splay state does no longer occur and the in-phase state is recovered Figure 3A. In contrast, for high age parameter β the splay state occurs even without tumor and multifrequency clustering in a healthy state (see Figure 3E). The two different pathological scenarios in Figures 3C,D) are connected with different responses of the immune layer and might indicate different malignity of the tumor. A map of regimes in the parameter plane of (β, r) is shown in Figures 4A,B, where the ensemble average s μ of the standard deviation of the mean phase velocities is plotted for the parenchymal and immune layer, respectively. Bright colors correspond to large deviations of the frequencies (mean phase velocities), and hence to pathological two-cluster states. Their enhancement with increasing age parameter and tumor size can be clearly seen. In Figures 4C,D, a comparison between the ensemble average of the standard deviation and the normalized standard deviation is depicted. In comparison to the standard deviation, the normalized standard deviation shows higher values of the ensemble average s 1 in case of splay states (see Figure 4C for β > 0.6π). On the contrary, the small bulge of s 1 at β 0.5π is less pronounced. Figure 5 depicts the comparison between the normalized standard deviation of the parenchymal and immune layer. It clearly shows the strong increase of the frequency deviation in the parenchyma associated with the pathogenicity of the tumor with increasing age parameter, while the much smaller bulge in the immune layer is associated with the immune layer activation before the tumor becomes clearly visible.

SEPSIS
In our model we focus on the stage of sepsis which is characterized by generalized inflammation. Depending on the individual subject and the local cellular situation of cytokine and fluid influx, either de-escalation of the inflammatory reaction and restoration of homeostasis, or organ failure occurs. Sepsis is usually preceded by a pre-septic perturbation of the parenchyma, e.g., by a wound which is infected by germs. This perturbation is terminated after a while by blocking off the wound by blood coagulation and eventually healing. Under "normal" conditions, the system returns to a healthy state; however, if sepsis occurs, an inflammatory immune response triggered by the infection spreads across the whole body (cytokine storm). As a consequence, the immune activity may invade large parts of the body through blood vessels and lead to severe organ failure and death. Whether sepsis terminates in a septic shock with severe consequences for the patient depends crucially on the ability of the immune and parenchymal system to regain homeostasis. As in the case of the tumor disease, the ability for returning to a healthy state, i.e., a frequency-synchronized state, is subsequently analyzed using the systemic sum parameter (age parameter) β.
For the simulations, we assume that a dysregulation of the cytokine activity in the immune system has already occurred due to a systemic immune response. The dysregulation of the cytokine activity is modeled by an initial perturbation of the cytokine activity matrix κ 2 ij of the immune layer. For the latter, we consider FIGURE 5 | Ensemble average s μ of the normalized standard deviation of the mean phase velocities in dependence of β for a fixed value of r 7.5% for the parenchymal (solid curve, diamonds) and immune layer (broken curve, dots), respectively for an ensemble size of N E 200. A cubic Savitzky-Golay filter has been applied to smooth the curves. Other simulation parameters as in Figure 2.   (0), ϕ μ j (0). Within each layer μ the nodes are sorted first by 〈 _ ϕ 1 j 〉, then by ϕ 1 j . Simulation parameters: σ 1, ω 1 ω 2 0. Other parameters as in Figure 2 and the initial conditions as in Figure 6. Frontiers in Network Physiology | www.frontiersin.org January 2022 | Volume 1 | Article 730385 a separation into two clusters with high activity between nodes from the same cluster (κ 2 ij 1) and no activity between nodes from different clusters (κ 2 ij 0). All other initial conditions are chosen randomly, as in the simulations for tumor disease in Section 4. An example initial condition is presented in Figure 6. Note that we arbitrarily fix the cluster sizes of the initial cytokine activities in the immune layer for the rest of this paper. Furthermore, we increase the value of the interlayer coupling strength σ compared to our simulations of the tumor disease. This choice is natural in order to understand the mechanism in action during the progression of sepsis. As described in Section 2, proinflammatory cytokines act on endothelial cells and hence cause an increased blood vessel leakiness. As a result, more immune cells and cytokines enter the stroma which in consequence increases the immune-parenchymal interaction.
In the following, we present simulation results for different choices of the age parameter β showing that after an initial cytokine perturbation in the immune layer either the healthy frequency-synchronized state is restored, or the whole system goes to a desynchronized or multifrequency cluster state. Figure 7 shows five representative dynamical scenarios that are induced by pathological initial cytokine activity in the immune layer. For an age parameter β 0.5π in Figure 7A, we observe that the system relaxes to a healthy state, i.e., frequency synchronized state, after an initial immune layer perturbation. Hence, the pro-inflammatory response is stopped. A transient of 1,000 time units is discarded before the temporal average over T 1,000 is taken in the mean phase velocities 〈 _ ϕ μ j 〉. In addition, Figure 8A shows that in the final state the phases in both layers are in-phase synchronized representing a resilient healthy state. The same scenario can be observed as well for an increased value of the age parameter β 0.6π in Figures 7B, 8B. While for β 0.5π it is unlikely to obtain a pathological state, i.e., a desynchronized state, from any random initial conditions of κ 1 ij and the phases ϕ μ j , the probability of a pathological state for an ensemble of random initial conditions κ 1 ij (0), ϕ μ j (0) increases for β 0.6π. In Figures 7C, 8C, we depict a dynamical scenario where the initial immune layer perturbation induces a desynchronization of the parenchymal layer. We clearly observe the presence of a two-frequency cluster. In this situation, the initial activated immune response can not be compensated by the coupled system and pushes the parenchyma away from a homeostatic state that may have severe consequences for the organic tissue, compare with the discussion for tumor disease in Section 4. For even higher values of the age parameter β 0.7π the probability of desynchronization increases. Depending upon the random initial conditions of κ 1 ij and the phases ϕ μ j , also for this value of β, we may obtain a frequency-synchronized (healthy) or a desynchronized (pathological) state. In Figures 7D, 8D, we display a frequency-synchronized state for β 0.7π. However, this particular state is not in-phase synchronized any more, as it is the case for the healthy states in Figures 8A,B. Figure 8D shows this state which possesses a splay distribution of the phases. Splay distributions of the phases may be interpreted as more vulnerable and less resilient against further perturbations by pathological cells in the parenchymal layer or a "second hit" phenomena known for sepsis. The more likely dynamical scenario of an emergent pathological two-frequency cluster for the same value of β 0.7π is displayed in Figures 7E, 8E.
For all presented scenarios in Figures 7, 8, we note that the frequencies of the parenchymal and the immune layer are locked. The emergence of the dynamical phenomenon of locking Pikovsky et al. (2001) can be explained by the increased value of the interlayer coupling strength σ compared to the simulation of the tumor disease, see Section 4. Locking can be interpreted as the dynamical manifestation of the physiological observation that the immune system is taking control of the whole system in case of sepsis. For smaller interlayer coupling strength this locking cannot be observed, and hence desynchronization induced by a pathologically activated immune system is not possible.
The dynamical scenarios discussed in this section are very similar to the observations presented for tumor disease. Additionally, in agreement with the results found for tumor disease, the desynchronization of the parenchymal layer due to an initially activated immune layer is observed for higher age parameter β with increasing probability. We note that also here resilient and vulnerable healthy states may coexist with pathological desynchronized states. The increased desynchronization probability with increasing β is reflected in the increasing ensemble average of the normalized standard deviation of the mean phase velocities in Figure 9.

CONCLUSION
In this article, we have proposed a two-layer network model for carcinogenesis and sepsis based upon the interaction of parenchymal cells and immune cells via cytokines and the co-evolutionary FIGURE 9 | Ensemble average s μ of the normalized standard deviation of the mean phase velocities in dependence of β for a fixed value of σ 1 for the parenchymal and immune layer (both curves coincide) for an ensemble size of N E 200. A cubic Savitzky-Golay filter has been applied to smooth the curve. Other simulation parameters as in Figure 7.
Frontiers in Network Physiology | www.frontiersin.org January 2022 | Volume 1 | Article 730385 dynamics of parenchymal, immune cells and cytokines. Certain parallels between cancer and infectious disease have been unveiled in the medical sciences, however, little is understood about the underlying mechanism behind these very similar pathological states induced by tumor and sepsis, respectively Hotchkiss and Moldawer (2014). With this study, we propose a novel paradigm of unified functional modeling of tumor disease and sepsis from the complex dynamical network perspective by choosing the nonspecific innate immune system as the reference point for both diseases. The proposed model is not a detailed model of organs but a functional model of dynamic interactions. Here, cytokine activity is described as information flow within and between the parenchymal and the immune layer. In particular, the communication between cells of the parenchyma, and between cells of the immune layer is modeled by adaptive coupling weights. This approach is complementary to works modeling the cytokine concentrations as additional dynamical species Yiu et al. (2012). Thus our perspective accounts for the cytokine activation rather than their physical mass.
In this paper we have presented the simplest form of a twolayer network model based upon adaptively coupled phase oscillators. Although many simplifying assumptions have been made, this model can already capture essential unifying features of tumor disease and sepsis. Two important tunable model parameters have been identified, i.e., the tumor size r denoting the fraction of mutant (pathological) parenchymal cells which have a deviating natural frequency of single-cell basic metabolic activity, and an adaptivity parameter β (called age parameter) summarizing various physiological conditions like age, inflammaging, adiposity, pre-existing illness, physical inactivity, nutritional influence, and malignancy in case of tumor disease. The healthy system is modeled as a completely frequency-synchronized homeostatic state of metabolic activity in both the parenchyma and the immune layer, while in the pathological case the parenchyma splits into two frequency clusters with different collective frequency, one corresponding to the healthy part, and one corresponding to organ failure. The desynchronization of the parenchymal layer is observed for higher age parameter β with increasing probability both in case of tumor disease and sepsis.
To characterize the dynamical patterns resulting from this model, we have introduced the temporally averaged mean phase velocities of the cells as collective frequencies of metabolism, and their ensemble-averaged standard deviation s μ as a measure of pathogenicity in case of the parenchymal layer (μ 1), and as a measure of metabolic activation in case of the immune layer (μ 2). Thus s 1 ≠ 0 is an immediate indicator of a pathological state. Further, we have discussed the initial perturbation for both types of pathological conditions, i.e., tumor and sepsis. For the analysis of the tumor disease, we assume heterogeneity in the natural frequency distribution of the parenchymal layer to account for the fact that mutated parenchymal cells may possess a faster metabolism than healthy parenchymal cells. For sepsis, we consider, however, a homogeneous network of oscillator with an initially activated cytokine matrix representing a systemic immune response. Moreover, for sepsis, we assume a high coupling strength between the parenchymal and the immune layer caused by the pro-inflammatory immune response.
In case of tumor disease, our simulations show that for small values of the age parameter β, as representative for young and healthy subjects, almost all simulations lead to a healthy state even for a relatively high fraction r of mutated cells. Thus the system is robust against random mutations of parenchymal cells. For increasing age parameter, the immune layer turns out to be more activated in order to keep the parenchymal layer synchronized. Here, only for a low fraction of mutated cells r the healthy state can be maintained. For higher values of the fraction r the immune layer is not able anymore to keep the parenchyma synchronized, and a pathological desynchronized state emerges. For high values of the age parameter the immune layer is completely dysregulated, and a severe pathological state of tumor disease is almost inevitable. One might speculate that in the in-phase scenario the pathological cluster (tumor) may be removed (by surgery) so that the healthy state is restored, whereas in the splay scenario the non-tumor cluster is already so strongly perturbed that removal of the tumor may not bring the system back to the healthy state.
In analogy with the study of tumor disease, we have demonstrated various dynamical scenarios for different age parameters in case of sepsis. Also here, the age parameter β has been shown to be critical for the system to return to the healthy state. For specifically chosen parameter values, close to those chosen for the examples of tumor disease, we have visualized the emergence of different synchronized and desynchronized states. Furthermore, we have found parameter regimes for which healthy states may coexist with pathological desynchronized states, i.e., it depends upon the random initial conditions which of the two states is asymptotically reached. Further systematic studies should investigate the detailed dependence upon the initial perturbation of the cytokin activity matrix in the immune layer, i.e., the initial cluster and its size, and upon the choice of the interlayer coupling strength. We hypothesize that the interlayer coupling strength σ may be used to model functionally the progression of sepsis where an increased σ may be interpreted as an ongoing in-stream of cytokines from the blood vessels into the stroma. For both tumor disease and sepsis our model represents only a first step towards more detailed and systematic investigations.
In summary, a unified disease model has been established using the innate immune system as the reference point. Diseases such as tumorigenesis and sepsis have been modeled by studying the dynamic functional interactions between parenchyma and stroma. For this purpose, a unified system of dynamic equations has been developed, in which only three different parameters, the age parameter, the number of mutated cells, and the initial cytokine activation are varied. Thereby, carcinogenesis, organ dysfunction in sepsis, and recurrence risk can be described in a correct functional context. Our studies open a perspective for further dynamic disease modeling. For instance, more detailed models, including disease progression may be developed in future studies.