Complex Diseases

1. Introduction

What is a disease? How does it emerge? Can we define a theoretical framework to explain how diseases persist and transform? The answers to these questions have evolved throughout history, reflecting advances in medicine and biology as well as changes in a wide variety of human societies and environments [1].

Traditionally, human conditions have been understood to arise from two distinct causes: internal failures of regulation (related to development, aging, or homeostasis) or external factors (such as pathogens, toxins, or environmental stressors). In practice, however, many important disorders cannot be fully understood based solely on such causes. Instead, disease is increasingly viewed as the interaction between intrinsic susceptibility and external perturbations, linking pathology to organismal, ecological, and evolutionary dynamics [2], which challenges reductionist explanations centered on single genes, pathogens, or physiological defects.

A promising alternative approach involves framing human pathology not as a localized malfunction, but as a complex disease. That is, an emergent dynamical state of interconnected systems driven by nonlinear interactions across multiple scales, from molecules and cells to organisms, populations, and environments.

This reality highlights the need for theoretical models that capture the fundamental properties of such non-trivial collective systems. As illustrated in Figure 1, complex systems and network science provide quantitative tools for mapping and analyzing these properties, from multiscale interactions and feedback loops to systemic instabilities in disease [3,4].

From this perspective, we explore how concepts such as attractors, criticality, homeostasis, and resilience provide a common framework for understanding diverse complex diseases [5,6,7,8]. Specifically, we discuss how mathematical models can help characterize pathological states and the transitions between health and illness. In addition, we show how these models reflect common dynamical principles underlying disease emergence, persistence, and recovery across comorbidity, cancer, infectious diseases, and neurological disorders.

2. Diseases as Networks

Organ systems form a highly integrated architecture that maintains physiological function through multiscale, self-regulated coordination mechanisms [9,10]. This is a recurrent pattern in biology, from gene regulation and protein interactions to brain connectomes and ecosystems [11]. However, this integration also represents an important vulnerability, since it allows local failures to propagate across multiple body systems [12].

Chronic obstructive pulmonary disease, for instance, not only affects the lungs: it induces systemic inflammation and cardiovascular complications [13]. Similarly, diabetes begins as a metabolic disorder but progressively alters renal, neural, and vascular systems [14]. Even Alzheimer’s disease (AD), the most common neurodegenerative condition, involves widespread effects beyond the nervous system, such as metabolic and inflammatory alterations [15]. Hence, considering these interactions across organs and functions is key to understanding how diseases emerge, persist, and must be treated [16].

Such feedbacks are not limited to organic comorbidity, but can also describe how neurological and behavioral disorders manifest: diffuse signs, like anxiety or stress, do interact in non-trivial ways and can induce major changes in the patient’s state of mind, even when latent causes appear to have subsided [17,18,19,20]. A clear example of this is major depression (MD), a mental health disorder that affects nearly 17% of the population [21].

Network science provides a natural framework to study these phenomena from a system-agnostic perspective. An effective dynamical model can be constructed from the following principle: symptoms are not independent manifestations of isolated problems, but interacting components within an underlying dynamical system [22]. A symptom here is understood as an observable (e.g., fever, fatigue, or pain) that becomes active when the latent biological processes drive it across a threshold and thus expresses as a phenotype. Hence, symptom interactions reflect statistical co-activation patterns driven by the underlying biological processes.

Let $\mathcal{G}=(V,E)$ be a graph whose nodes represent symptoms connected through edges $e_{ij}$, which can be binary or weighted (Figure 1a). We define a set of symptoms $\mathcal{S}=\{S_1,\dots ,S_n\}$, where $S_i\in \{0,1\}$ indicates the absence or presence of a given symptom. Symptom interactions are encoded in a weighted adjacency matrix $W=\left\{w_{ij}\right\}$ defining the interaction graph (Figure 1b). The total input received by symptom i is

$$h_i\left(t\right)=c\sum _{j\in {\mathsf{\Gamma }}_i}{w}_{ij}{S}_j\left(t\right)+\eta ,$$

(1)

where ${\mathsf{\Gamma }}_i$ denotes the neighborhood of node (i) (Figure 1b, shaded area), c controls the strength of symptom interaction, and $\eta$ is an external stress or vulnerability parameter that increases activation.

The probability that a symptom becomes active is then

$$\mathbb{P}\left(S_i(t+1)=1\right)={\displaystyle \frac{1}{1+exp\left[-\left(h_i\left(t\right)-{\theta }_i\right)\right]}},$$

(2)

where ${\theta }_i$ is the corresponding activation threshold. Subsequently, the system’s global state can be characterized by defining the mean activity $x\left(t\right)={\displaystyle \frac{1}{n}}{\sum }_{i=1}^n{S}_i\left(t\right)$: low values of x correspond to healthy states, while high values indicate widespread activation of symptoms.

A useful approximation can be obtained by replacing the heterogeneous network with a homogeneous parameter. In this homogeneous mean-field approximation, ${\sum }_{j\in {\mathsf{\Gamma }}_i}{w}_{ij}{S}_j\left(t\right)\simeq \alpha x\left(t\right)$, so that

$$\mathbb{P}\left(S_i(t+1)=1\right)\simeq {\displaystyle \frac{1}{1+exp\left[-\left(c\alpha x\left(t\right)+\eta -\theta \right)\right]}},$$

(3)

with an effective global activation threshold $\theta$.

Thus, $\eta$ shifts the activation threshold effectively to ${\theta }_{\mathrm{eff}}=\theta -\eta$, which increases the probability of activation even when symptom reinforcement is weak. If we assume that the average symptom activity $x\left(t\right)$ relaxes continuously toward its expected activation level with characteristic timescale $\tau$, this yields

$$\tau {\displaystyle \frac{dx}{dt}}=-x+{\displaystyle \frac{1}{1+exp[-(c\alpha x+\eta -\theta \left)\right]}}.$$

(4)

The first term describes recovery, the second captures reinforcement through the symptom network, and c controls the overall stress level or symptom connectivity.

Starting from a healthy state (A), gradual changes in stress eventually drive the system to a tipping point (B), where a sudden transition to the disease state (C) occurs. Close to the transition, the dynamics can be approximated by a cubic equation, whose fixed points define the folded surface shown in Figure 1d:

$${\displaystyle \frac{dx}{dt}}\approx x(x-x_{-})(x_{+}-x):=f(x,c,\eta ).$$

(5)

Because of hysteresis, reducing stress does not restore health immediately: recovery follows a different path $C\to D\to A^{\prime }$ and only occurs after crossing a second threshold (D). Consequently, the transition to pathology becomes a critical phenomenon where small perturbations near tipping points can trigger abrupt and long-lasting changes in system behavior. This simple theoretical approximation allows studying the stochastic dynamics of transitions between healthy and pathological states, loss of resilience, and disease chronification.

Another way of looking at this system is to consider the associated potential function $V_c\left(x\right)$, which allows us to rewrite the model as derived from a potential, namely:

$${\displaystyle \frac{dx}{dt}}=-{\displaystyle \frac{d{V}_{c\eta }\left(x\right)}{dx}},V_{c\eta }\left(x\right):=-{\int }^{x}f(x^{\prime },c,\eta )d{x}^{\prime }.$$

(6)

This potential does not represent a physical energy but a stability landscape: its minima are attractors (healthy and disease states), and their depth encodes resilience to perturbations [23]. Thus, this formalism can also characterize the system’s phase transitions [24].

In our case, $V\left(x\right)=-\int x^{\prime }(x^{\prime }-x_{-})(x_{+}-x^{\prime })d{x}^{\prime }$ gives a quartic function, analogous to Figure 1d (bottom) for two given parameter values associated with healthy states (H) and disease (D). When the healthy attractor becomes shallow due to increasing stress, the system becomes metastable, and fluctuations can cause a sudden transition ($H\to D$). Because the disease attractor may remain stable even after stress is reduced, recovery requires crossing a second threshold, thereby generating hysteresis and persistent pathological states.

3. Ecology and Disease

A key feature underlying the hallmarks of cancer [27,28] is intratumoral heterogeneity, which generates the genetic and phenotypic diversity fueling adaptation, metastasis, and therapy resistance [29]. Tumors behave as evolving ecosystems composed of diverse interacting cell populations (Figure 1e,f) where competition, cooperation, and niche construction shape disease progression [30]. Thus, cancer exemplifies the interplay between ecology, evolution, and development in disease dynamics [31].

Understanding how this diversity arises and persists is therefore a central ecological problem. Mathematical models have shown how interactions among tumor, host, sensitive, and resistant cell populations, together with spatial constraints, can maintain heterogeneity and drive tumor evolution [26,32,33]. These studies have also inspired ecology-based therapeutic strategies that explicitly account for these evolutionary and ecological interactions.

Since genetic and phenotypic diversity in tumors can be extremely high [29], low-dimensional models might miss certain emergent properties of the disease. Instead, a natural step to model tumors as communities of many interacting phenotypes is to describe them as multispecies ecosystems using generalized Lotka–Volterra (gLV) equations (Figure 1g, [34,35]):

$${\displaystyle \frac{d{c}_{\mu }}{dt}}=r_{\mu }{c}_{\mu }\left(1-{\displaystyle \frac{c_{\mu }}{K_{\mu }}}+{\displaystyle \sum _{\nu \ne \mu }^{N\left(t\right)}}{A}_{\mu \nu }{c}_{\nu }\right).$$

(7)

Here, $r_{\mu }$ and $K_{\mu }$ define the growth and ecological adaptation of phenotype $\mu$ within the tumor microenvironment, while $A_{\mu \nu }$ captures interactions among phenotypes [36,37]. Importantly, the interactions involving competition [38], predation [39], and facilitation [40] can vary widely in sign and magnitude across tumors, which leads to highly heterogeneous ecological networks.

Such diversity of cancer cells and their interactions has major consequences. When interactions are strong and heterogeneous, gLV models predict rugged landscapes with multiple attractors [34], allowing abrupt transitions between ecological states, including drug-resistant configurations [35] (Figure 1h). This helps explain why maximum-dose therapies often fail: treatment may not only eliminate sensitive cells, but also perturb the ecosystem and drive transitions toward resistant states.

Furthermore, recent work in bacterial ecology shows that coexistence among many competing species can persist under strong competitive constraints [41,42], potentially explaining the remarkable heterogeneity observed in tumors. Overall, these results suggest that tumors display emergent properties characteristic of complex ecosystems, limiting the effectiveness of reductionist “magic-bullet” therapies.

4. Thresholds and Epidemic Spreading

Although infectious diseases may be thought of as continuous and gradual processes, a key characteristic of epidemic and pandemic outbreaks is that they occur only when transmission exceeds critical thresholds. Since these sharp transitions determine whether the pathogen spreads among susceptible populations, infectious diseases constitute a paradigmatic example of tipping point dynamics with profound biological, economic, and social implications [43,44]. Thus, studying these thresholds is paramount to epidemiology, public health, and disease prevention through vaccination.

The susceptible-infected-susceptible (SIS) model provides a formal framework to study infections that do not confer long-lasting immunity. This has been used not only for biological infections, but also for other contagion-like processes in networks [45].

Consider a contact network (Figure 1j) described by an adjacency matrix $A_{ij}$, where $A_{ij}=1$ if the nodes i and j interact and $A_{ij}=0$ otherwise. Let ${\rho }_i\left(t\right)$ denote the probability that node i is infected at time t. In the most general case, transmission and recovery rates may vary between nodes and links. Denoting by ${\beta }_{ij}$ the transmission rate from node j to node i, and by ${\mu }_i$ the recovery rate of node i, the SIS mean-field dynamics is (Figure 1k):

$${\displaystyle \frac{d{\rho }_i}{dt}}=-{\mu }_i{\rho }_i+\left(1-{\rho }_i\right){\displaystyle \sum _{j=1}^N}{\beta }_{ij}{A}_{ij}{\rho }_j$$

(8)

($i=1,\dots ,N$). The first term represents recovery processes, while the second term accounts for infection coming from neighboring infected nodes. This formulation allows for heterogeneous susceptibilities, heterogeneous infectivities, weighted interactions, and nonuniform recovery dynamics [46]. The disease-free equilibrium is given by ${\rho }_i^{*}=0,\forall i=1,\dots ,N$.

A simpler and analytically transparent situation is obtained by assuming homogeneous spreading and recovery dynamics, namely ${\beta }_{ij}=\beta$ and ${\mu }_i=\mu$, for all nodes and links. Equation (8) then reduces to

$${\displaystyle \frac{d{\rho }_i}{dt}}=-\mu {\rho }_i+\beta \left(1-{\rho }_i\right){\displaystyle \sum _{j=1}^N}{A}_{ij}{\rho }_j.$$

(9)

The disease-free state again corresponds to ${\rho }_i^{*}=0$. Linearization gives $d\rho /dt=(\beta A-\mu I)\rho$, implying that the epidemic threshold is controlled by the largest eigenvalue ${\Lambda }_{max}\left(A\right)$ of the adjacency matrix:

$${\lambda }_c={\displaystyle \frac{\beta }{\mu }}={\displaystyle \frac{1}{{\Lambda }_{max}\left(A\right)}}.$$

(10)

For $\lambda <{\lambda }_c$, outbreaks remain subcritical and decay to the disease-free state. For $\lambda >{\lambda }_c$, the infection survives, and an endemic state emerges.

To gain insight into the epidemic threshold, consider a homogeneous network where all nodes have approximately the same (average) degree $\langle k\rangle$. In this limit, the node-dependent probabilities can be replaced by a single average infected fraction, ${\rho }_i\left(t\right)\approx \rho \left(t\right)$, leading to the mean-field equation

$${\displaystyle \frac{d\rho }{dt}}=-\mu \rho +\beta \langle k\rangle \rho (1-\rho ).$$

(11)

The fixed points are the disease-free state, ${\rho }_0^{*}=0$, and the endemic state

$${\rho }_1^{*}=1-{\displaystyle \frac{\mu }{\beta \langle k\rangle }}=1-{\displaystyle \frac{1}{\lambda \langle k\rangle }},$$

(12)

which exists only if $\lambda \langle k\rangle >1$. The epidemic threshold is therefore

$${\lambda }_c={\displaystyle \frac{1}{\langle k\rangle }}⟺{\beta }_c={\displaystyle \frac{\mu }{\langle k\rangle }}.$$

(13)

Below this threshold, infections disappear, whereas above it, a stable endemic state emerges continuously:

$${\rho }^{*}=1-{\displaystyle \frac{{\lambda }_c}{\lambda }}.$$

(14)

The infected fraction $\rho$, then, is the order parameter, while the effective spreading rate $\lambda =\beta /\mu$ acts as the control parameter that drives the transition [43].

This mechanism provides the theoretical basis for vaccination strategies that reduce transmission below the epidemic threshold. Nevertheless, real contact networks are strongly heterogeneous, since they often contain hubs and highly connected individuals that drastically lower or even eliminate the effective threshold [46]. Accordingly, epidemic spreading and vaccination outcomes may differ substantially from homogeneous mean-field predictions, which highlights the importance of network structures in epidemiology.

Box 1. Networks in health and disease

Health and disease can be represented as dynamical states of interacting molecular, cellular, and other variables of clinical interest. In network medicine, genes, proteins, metabolites, cells or phenotypes are modeled as nodes, and their biochemical, regulatory or functional dependencies are links. Accordingly, disease is not simply assigned to a single altered component, but to perturbations of a subnetwork: a disease module embedded in the interactome [47] or in a multilayer representation of multi-omic data [4,48,49]. This framework allows us to explain why mutations, environmental exposures or drugs may have system-wide effects that depend on network context, including the systemic importance of a node (e.g., centrality) and its synergetic relationships with other nodes defining a specific (dys)function (e.g., module membership). Taking into account the underlying interdependencies among relevant biological processes across multiple scales, from molecular to organ, allows one to map the propagation of state perturbations that can lead to altered states characterizing disease.

A more general network-constrained dynamics [4] can be written as

$$\begin{array}{ccc}\hfill {\dot{x}}_i\left(t\right)& =& F_i\left(x_i\left(t\right),{\mathbf{x}}_{\partial i}\left(t\right),{\theta }_i,{\mathbf{u}}_{\partial i}\left(t\right),t\right)+{\eta }_i\left(t\right)\hfill \\ & \approx & f_i\left(x_i\left(t\right),{\theta }_i,u_i\left(t\right)\right)+{\displaystyle \sum _{\alpha =1}^L}{\displaystyle \sum _{j=1}^N}{A}_{ij}^{\left[\alpha \right]}{g}_{ij}^{\left[\alpha \right]}\left(x_i\left(t\right),x_j\left(t\right),{\theta }_{ij}^{\left[\alpha \right]}\right)+{\eta }_i\left(t\right),\hfill \end{array}$$

(15)

where $\partial i$ denotes the adjacency set of unit $i=1,2,...,N$, $x_i\left(t\right)$ is its state, and ${\mathbf{x}}_{\partial i}\left(t\right)=\{x_j\left(t\right):j\in \partial i\}$ denotes the states of its neighbors. The first line is fully general: the dynamics of node i may depend on its own state, on the joint state of all its neighbors, on local and contextual parameters, and on perturbations. The second line is the pairwise additive approximation used in many current formulations, assuming that the effect of the network can be decomposed into edge-level contributions $g_{ij}^{\left[\alpha \right]}$ across layers $\alpha =1,2,...,L$, weighted by the strength of the interaction encoded by $A_{ij}^{\left[\alpha \right]}$. Here, $u_i\left(t\right)$ represents genetic, environmental, or therapeutic perturbations, ${\eta }_i\left(t\right)$ captures uncertainty, ${\theta }_i$ and ${\theta }_{ij}^{\left[\alpha \right]}$ denote, respectively, node-specific and interaction-specific parameters, which may encode molecular, cellular, disease- or patient-dependent variability.

Within this framework, network analysis can identify disease modules, biomarkers, drug targets, and repurposing opportunities, and can support patient stratification by linking multi-omic variation to pathophysiological mechanisms and clinical outcomes [50,51,52,53,54].

5. Disease as Multiscale Systems

Human disease involves processes where many components interact in intricate ways, and major changes can occur abruptly. Decomposing these processes into simpler parts and identifying a molecular cause has previously led to extraordinary results in biomedicine [55]. Nevertheless, many important conditions resist this decompositional approach because they are rooted in highly interconnected physiological systems where the underlying processes are deeply interdependent [48]. This leads to emergent properties: behaviors that arise from interactions among components and cannot be predicted by examining those components in isolation [56].

In these systems, moving from one scale to the next introduces new variables and laws. As a result, higher-level properties cannot be straightforwardly deduced from a complete description of lower-level components. Causality flows both upward and downward, since a molecular disturbance can reorganize cellular dynamics, but tissues or system-level states can equally constrain and reshape molecular processes [57].

Multiscale models provide a framework to approximate such emergent properties [58]. Instead of a single unified description, they rely on a structured hierarchy of scale-dependent formalisms, which can be connected to study how disturbances initiated at one scale propagate, amplify, or dampen as they manifest in other levels [59].

As shown in Figure 2, AD is a clear example of this kind of multi-scale process. To further illustrate this, we discuss formal representations that can be used to study its relevant mechanisms across various scales.

From the molecular perspective, the accumulation of amyloid-$\beta$ plaques and tau neurofibrillary tangles is considered central to the pathology [60]. This can be represented as reversible reactions among monomers $\left(M\right)$, dimers $\left(D\right)$, and trimers $\left(T\right)$, whose downstream dynamics govern the accumulation of plaques $\left[P\right]$ and fibrils $\left[F\right]$:

$$\begin{array}{c}\hfill {\displaystyle \frac{d\left[P\right]}{dt}}={\lambda }_1\left[D\right](K_1-\left[D\right])-{\mu }_4\left[P\right],\end{array}$$

(16)

$$\begin{array}{c}\hfill {\displaystyle \frac{d\left[F\right]}{dt}}={\lambda }_3{\displaystyle \frac{\left[P\right]}{1+\left[N\right]/K_3}}-{\mu }_6\left[F\right].\end{array}$$

(17)

This model provides a basis for why amyloid-$\beta$ burden does not reliably predict clinical onset, since fixed-point analysis reveals regimes where self-limiting pathology and irreversible growth coexist as stable solutions.

But what about the spatially heterogeneous patterns that result from aggregation dynamics? Local plaque formation at the tissue scale can be modeled on a 2D discrete lattice as a density $N_t$ that grows logistically under a carrying capacity $N_{\infty }$, with initial density $N_0$ and growth rate A. In turn, disaggregation probability obeys the following law [61]:

$$P_{\mathrm{dis}}(t+1)=P_{\mathrm{dis}}\left(t\right)+w·{\displaystyle \frac{N_t-N_{t-1}}{V}}.$$

(18)

Parameters A and w cannot be directly derived from the molecular equations: they absorb the collective effects of neuroinflammation, synaptic loss, and glia reactivity, which ultimately influence whole-brain organization.

At that scale, the local dynamics of each brain region (j) can be modeled as a Stuart-Landau oscillator [62]

$${\displaystyle \frac{d{z}_j}{dt}}=z_j\left[a_j+i{\omega }_j-{\left|z_j\right|}^2\right]+{\eta }_j,$$

(19)

where the bifurcation parameter $a_j$ controls whether a region sustains oscillations $(a_j>0)$ or settles to a fixed point $(a_j<0)$. Consequently, neurodegeneration acts as a distributed reduction in $a_j$, progressively shifting the network toward a fragmented, incoherent activity regime. Such pathological changes eventually manifest in the behavioral plane, including cognitive decline, memory loss, confusion, and personality changes [63].

Cognitive deficits can be formally linked to the topological disintegration of the patient’s functional network. In particular, defining a thresholded adjacency matrix $A_{ij}\left(w_c\right)=\mathbf{1}\{w_{ij}\ge w_c\}$, the size of the largest connected component $S\left(w_c\right)=\left|C_{max}\right|/N$ serves as an order parameter for network integrity. If we assume an uncorrelated random network [64], percolation theory predicts the following critical threshold:

$$p_c={\displaystyle \frac{\langle k\rangle }{\langle k^2\rangle -\langle k\rangle }}.$$

(20)

Below this threshold, the functional network fragments abruptly. The heterogeneity term $\langle k^2\rangle -\langle k\rangle$ implies that hub-rich networks are more resilient to random failures yet disproportionately sensitive to targeted attacks on high-degree nodes [65,66]. This is precisely the failure mode that AD appears to exploit in its early stages.

6. Discussion

The perspective of complex systems suggests that human disease should not be viewed as a static entity or the result of isolated failures, but as emergent dynamics arising from interactions across multiple scales.

As discussed in the case studies above, elemental models show that abrupt transitions are key: bistability and hysteresis in symptom networks, thresholds governing epidemic outbreaks, or amyloid-$\beta$ aggregation dynamics, all reflect the importance of tipping points separating qualitatively different regimes, even in seemingly gradual processes. Health and disease, then, correspond to regions in an underlying dynamic landscape, and the onset of disease is often associated with a loss of resilience or transitions between alternative attractors [67].

Remarkably, heterogeneity can reshape thresholds in ways that homogeneous descriptions cannot capture. The ecological interactions among tumor subpopulations generate rugged landscapes with attractors that drive drug resistance, while intricate network structures can drastically reduce or even eliminate the thresholds predicted by mean-field models. Furthermore, many of these processes display emergent properties, so multiscale descriptions are required to characterize their behavior.

Addressing these challenges will require modeling extensions that connect mechanisms across scales and capture heterogeneity by integrating data from real patients, which has already begun to yield promising results in diseases such as AD and cancer [68,69].