Multiplex Bi-Virus Reaction-Diffusion Models

• MBRD models are a family of frameworks that describe the spread of two interacting pathogens on multiplex networks with shared node sets and distinct layer dynamics.
• They combine within-node reaction processes such as transmission and recovery with layer-specific diffusion operators to capture phenomena like superinfection and coinfection.
• The models reveal critical thresholds, multiple equilibria, and diffusion-driven pattern formation, offering insights into epidemic control and mitigation strategies.

Multiplex Bi-Virus Reaction-Diffusion Models (MBRD) are a family of models for the spatio-temporal spread of two interacting pathogens or contagions on node-aligned multiplex networks, combining within-node reaction terms—such as transmission, recovery, virulence, competition, super-infection, or co-infection—with layer-specific transport operators. In the literature, MBRD spans several mathematically distinct but compatible formalisms: competitive continuous-time bi-virus SIS systems on directed multiplex infection graphs, exact fluid-limit ODEs on multipartite metapopulations, graph and PDE reaction-diffusion systems with supra-Laplacians and cross-diffusion, threshold-based co-diffusion models on multiplex social networks, and shared-resource epidemic systems with exclusive infection constraints [(Liu et al., 2019); (Santos et al., 2013); (Asllani et al., 2014); (Chang et al., 2018); (Janson et al., 2020); (Yu et al., 21 Aug 2025); (Yu et al., 3 Sep 2025)].

1. Network representations and state spaces

A recurrent structural assumption is that all layers share the same node set while allowing layer-specific connectivity. In the continuous-time competitive SIS backbone, virus $1$ and virus $2$ spread on their own directed, strongly connected infection graphs, with irreducible nonnegative infection matrices $B^{(1)},B^{(2)}\in\mathbb{R}^{n\times n}$, diagonal healing matrices $D^{(1)},D^{(2)}$, and heterogeneous node-wise rates. Infection fractions are $x_i(t),y_i(t)\in[0,1]$, susceptibility is $s_i(t)=1-x_i(t)-y_i(t)$, and no co-infection per individual is allowed, so $x_i(t)+y_i(t)\le 1$. The positively invariant domain is

$D=\{(x,y)\mid x\ge 0,\; y\ge 0,\; x+y\le 1\}.$

The neighbor graph is the union of the two infection graphs (Liu et al., 2019).

In the exact multipartite fluid-limit formulation, the node set is partitioned into islands $V_1,\dots,V_M$, with complete bipartite connectivity between connected islands and no intra-island edges. For $K$ strains, the limiting fractions $2$0 satisfy $2$1, again encoding exclusive infection through a shared susceptible capacity (Santos et al., 2013).

The explicit MBRD-SI and MBRD-CI frameworks use multiplex metapopulation networks with layers $2$2, $2$3, $2$4, and, in the co-infection case, $2$5. Here $2$6 denotes susceptible density, $2$7 pathogen-1 infected density, $2$8 pathogen-2 infected density, and $2$9 coinfection density. Diffusion acts on $B^{(1)},B^{(2)}\in\mathbb{R}^{n\times n}$0, $B^{(1)},B^{(2)}\in\mathbb{R}^{n\times n}$1, and $B^{(1)},B^{(2)}\in\mathbb{R}^{n\times n}$2 layers through graph Laplacians $B^{(1)},B^{(2)}\in\mathbb{R}^{n\times n}$3, while the $B^{(1)},B^{(2)}\in\mathbb{R}^{n\times n}$4 layer has no edges and no diffusion. Cross-diffusion appears as off-diagonal diffusion blocks in the $B^{(1)},B^{(2)}\in\mathbb{R}^{n\times n}$5 equation, driven by $B^{(1)},B^{(2)}\in\mathbb{R}^{n\times n}$6 and $B^{(1)},B^{(2)}\in\mathbb{R}^{n\times n}$7 layer Laplacians (Yu et al., 21 Aug 2025, Yu et al., 3 Sep 2025).

A distinct multiplex interpretation arises in shared-resource models. There, each virus has a person-person layer $B^{(1)},B^{(2)}\in\mathbb{R}^{n\times n}$8 and a person-resource layer defined by resource-to-node infection rates $B^{(1)},B^{(2)}\in\mathbb{R}^{n\times n}$9, node-to-resource contamination weights $D^{(1)},D^{(2)}$0, and resource cleaning rate $D^{(1)},D^{(2)}$1. The augmented matrices

$D^{(1)},D^{(2)}$2

embed the resource as an additional layer-like state while preserving exclusive infection at population nodes (Janson et al., 2020).

2. Competitive SIS backbone on multiplex graphs

The canonical continuous-time bi-virus competitive SIS dynamics on a network of $D^{(1)},D^{(2)}$3 groups are

$D^{(1)},D^{(2)}$4

or, in vector form,

$D^{(1)},D^{(2)}$5

Writing $D^{(1)},D^{(2)}$6 and $D^{(1)},D^{(2)}$7, this is equivalently

$D^{(1)},D^{(2)}$8

Under $D^{(1)},D^{(2)}$9 with $x_i(t),y_i(t)\in[0,1]$0, the set $x_i(t),y_i(t)\in[0,1]$1 is positively invariant (Liu et al., 2019).

This system inherits the single-virus SIS threshold structure. For

$x_i(t),y_i(t)\in[0,1]$2

if $x_i(t),y_i(t)\in[0,1]$3, then $x_i(t),y_i(t)\in[0,1]$4 is the only equilibrium and is globally asymptotically stable on $x_i(t),y_i(t)\in[0,1]$5; if $x_i(t),y_i(t)\in[0,1]$6, there are exactly two equilibria, $x_i(t),y_i(t)\in[0,1]$7 and a unique endemic equilibrium $x_i(t),y_i(t)\in[0,1]$8, with $x_i(t),y_i(t)\in[0,1]$9 asymptotically stable on $s_i(t)=1-x_i(t)-y_i(t)$0 (Liu et al., 2019).

The multipartite fluid-limit bi-virus system has the closely related form

$s_i(t)=1-x_i(t)-y_i(t)$1

which is an exact macroscopic limit of a normalized Markov jump process for large islands. In this representation, the adjacency-driven inflow is modulated by the same occupancy factor $s_i(t)=1-x_i(t)-y_i(t)$2, so the exclusivity constraint is structurally identical even though the microscopic derivation differs (Santos et al., 2013).

A common misconception is that all MBRD models are Laplacian diffusion systems. In the 2019 bi-virus model, “diffusion” is represented implicitly by adjacency-weighted infection terms $s_i(t)=1-x_i(t)-y_i(t)$3 rather than by explicit spatial Laplacians; the paper states that there are no explicit spatial Laplacian diffusion terms. The multipartite fluid-limit model likewise yields adjacency-gated infection inflows rather than gradient-driven fluxes [(Liu et al., 2019); (Santos et al., 2013)].

3. Explicit reaction-diffusion and cross-diffusion formulations

The explicit MBRD frameworks introduced in 2025 separate reaction kinetics from graph diffusion and cross-diffusion. For the super-infection model MBRD-SI, the supra-Laplacian block structure is

$s_i(t)=1-x_i(t)-y_i(t)$4

where $s_i(t)=1-x_i(t)-y_i(t)$5 are self-diffusion rates and $s_i(t)=1-x_i(t)-y_i(t)$6 are cross-diffusion coefficients acting only on the susceptible equation (Yu et al., 3 Sep 2025).

The node-wise MBRD-SI equations are

$s_i(t)=1-x_i(t)-y_i(t)$7

Superinfection enters through the replacement terms $s_i(t)=1-x_i(t)-y_i(t)$8: pathogen $s_i(t)=1-x_i(t)-y_i(t)$9 steals hosts from pathogen $x_i(t)+y_i(t)\le 1$0 (Yu et al., 21 Aug 2025).

MBRD-CI adds a dedicated coinfection compartment $x_i(t)+y_i(t)\le 1$1, transmission modifiers $x_i(t)+y_i(t)\le 1$2, coinfection virulence $x_i(t)+y_i(t)\le 1$3, and recovery $x_i(t)+y_i(t)\le 1$4. Diffusion again acts on $x_i(t)+y_i(t)\le 1$5, $x_i(t)+y_i(t)\le 1$6, and $x_i(t)+y_i(t)\le 1$7, while $x_i(t)+y_i(t)\le 1$8 remains reaction-only. The effective incidence denominator is $x_i(t)+y_i(t)\le 1$9, because $D=\{(x,y)\mid x\ge 0,\; y\ge 0,\; x+y\le 1\}.$0 and $D=\{(x,y)\mid x\ge 0,\; y\ge 0,\; x+y\le 1\}.$1 include $D=\{(x,y)\mid x\ge 0,\; y\ge 0,\; x+y\le 1\}.$2 (Yu et al., 21 Aug 2025, Yu et al., 3 Sep 2025).

These graph systems admit continuous-space analogs with homogeneous Neumann boundary conditions by replacing graph Laplacians with $D=\{(x,y)\mid x\ge 0,\; y\ge 0,\; x+y\le 1\}.$3. The susceptible PDE takes the form

$D=\{(x,y)\mid x\ge 0,\; y\ge 0,\; x+y\le 1\}.$4

while $D=\{(x,y)\mid x\ge 0,\; y\ge 0,\; x+y\le 1\}.$5 and $D=\{(x,y)\mid x\ge 0,\; y\ge 0,\; x+y\le 1\}.$6 satisfy standard diffusion equations and $D=\{(x,y)\mid x\ge 0,\; y\ge 0,\; x+y\le 1\}.$7 remains reaction-only (Yu et al., 3 Sep 2025).

A broader multiplex reaction-diffusion theory predates MBRD-SI and MBRD-CI. For two species $D=\{(x,y)\mid x\ge 0,\; y\ge 0,\; x+y\le 1\}.$8 on a two-layer multiplex,

$D=\{(x,y)\mid x\ge 0,\; y\ge 0,\; x+y\le 1\}.$9

$V_1,\dots,V_M$0

This formalism provides the supra-Laplacian and perturbative spectral machinery later reused in bi-virus settings (Asllani et al., 2014).

4. Thresholds, equilibria, and instability mechanisms

In the multiplex competitive SIS model, the healthy state $V_1,\dots,V_M$1 always exists and is the unique equilibrium, globally asymptotically stable on $V_1,\dots,V_M$2, if and only if

$V_1,\dots,V_M$3

If $V_1,\dots,V_M$4 and $V_1,\dots,V_M$5, there are exactly two equilibria: the unstable healthy state and a dominant virus-1 state $V_1,\dots,V_M$6, which is asymptotically stable with domain of attraction

$V_1,\dots,V_M$7

The symmetric statement holds for virus $V_1,\dots,V_M$8. If both spectral abscissae are positive, there exist at least three equilibria: the healthy state and the two single-virus endemic states (Liu et al., 2019).

Coexistence is strongly constrained by layer alignment. When both viruses share the same directed infection graph and homogeneous rates $V_1,\dots,V_M$9, $K$0, coexisting equilibria can exist only if

$K$1

At coexistence, irreducibility forces proportionality $K$2, $K$3. If $K$4, there are infinitely many coexisting equilibria parameterized by $K$5, and the Jacobian has a zero eigenvalue along the direction $K$6, so classical linear analysis is inconclusive (Liu et al., 2019).

The 2025 MBRD papers add diffusion-driven pattern formation to this equilibrium picture. At the disease-free equilibrium, $K$7 solves

$K$8

and the uniform-mode thresholds, ignoring diffusion, are

$K$9

with

$2$00

for the coinfection compartment in MBRD-CI (Yu et al., 3 Sep 2025).

Linearization around homogeneous equilibria yields cubic and quartic characteristic polynomials. For the three-morphogen MBRD-SI case, reaction-only stability requires

$2$01

Diffusion-driven instability is then diagnosed by the shifted mode matrices and the cubic discriminant

$2$02

For the four-morphogen MBRD-CI case, reaction-only stability is characterized by quartic Routh-Hurwitz inequalities, and Turing versus Turing-Hopf regimes are tied to the quartic discriminant $2$03 and the sign structure of the corresponding derivative polynomial (Yu et al., 21 Aug 2025).

Multiplex reaction-diffusion theory supplies an additional perturbative threshold. If the decoupled layers are stable and interlayer coupling is weak, the critical interlayer diffusion that triggers instability is approximated by

$2$04

This establishes that interlayer exchange can seed instabilities absent in the decoupled limit, while stronger coupling can also suppress patterns (Asllani et al., 2014). In the explicit MBRD studies, large diffusion disparities and negative cross-diffusion in the susceptible layer destabilize otherwise homogeneous equilibria on specific network modes (Yu et al., 3 Sep 2025).

5. Numerical regimes and observed phenomena

The competitive SIS backbone exhibits clear regime separation. With $2$05 and $2$06, both viruses are eradicated and trajectories converge to $2$07. With $2$08 and $2$09, virus $2$10 converges to its unique endemic equilibrium while virus $2$11 is eradicated, for all initial conditions except those with $2$12. When $2$13 and $2$14 on distinct layers, simulations show either coexistence or single-virus dominance depending on initial conditions (Liu et al., 2019).

The threshold-based co-diffusion model on a lattice/RRG multiplex reveals a different but related phenomenology. With $2$15 nodes, an $2$16 periodic lattice, an RRG-4 layer, $2$17, synchronous updates for $2$18 steps, and $2$19 Monte Carlo runs per parameter set, lower synergy makes contagions more susceptible to percolation, especially those that diffuse on lattices. Faster diffusion of one contagion with dormancy probabilistically blocks diffusion of the other in a “ring vaccination”-like manner, and within a band $2$20, lattice contagions can undergo bimodal or trimodal branching when they are the slower diffusing contagion (Chang et al., 2018).

In MBRD-SI and MBRD-CI, stationary Turing hotspots can form and grow. On lattice multiplexes such as LA4, LA12, and LA24, spotted or maze-like stationary patterns were observed, with morphology controlled by diffusion contrasts, negative cross-diffusion, and layer degrees. In SI, the average amplitude $2$21 at $2$22 on LA4-LA4-LA4 follows a power-law decay with superinfection strength $2$23, with fit $2$24, $2$25, $2$26. In CI, the average amplitude at $2$27 peaks around $2$28, indicating a non-monotone dependence of patterning on co-transmission. The “$2$29 threshold” increases approximately linearly with $2$30, with fit $2$31, $2$32, $2$33 (Yu et al., 3 Sep 2025).

Topology strongly modulates spread. Barabási-Albert networks consistently reach saturation faster than Watts-Strogatz networks for matched average degrees, while reducing infected mobility—implemented as low average degree in infected layers—consistently delays saturation in both SI and CI. In SI, delaying pathogen-2 saturation is best achieved when the $2$34 layer has low average degree, for example in $2$35 and $2$36 among the tested settings (Yu et al., 3 Sep 2025).

6. Sensitivity, control, interpretations, and limitations

For single-virus endemic equilibria with $2$37 and $2$38, linearization of the equilibrium equation yields

$2$39

and the matrix on the left has a strictly negative inverse. Consequently,

$2$40

elementwise. In the dominant-state bi-virus regime, $2$41 inherits the same monotonicity with respect to $2$42 and $2$43 (Liu et al., 2019).

The limitations of decentralized proportional control are sharp. For a single virus, $2$44 yields a structurally equivalent SIS system with $2$45, so the origin is unstable and a unique nontrivial endemic equilibrium exists. For the bi-virus system with $2$46 and $2$47, linearization at $2$48 remains unstable; proportional local feedback cannot stabilize the healthy state (Liu et al., 2019).

The shared-resource framework provides constructive mitigation strategies. If one chooses

$2$49

then $2$50 for all $2$51 yields asymptotic eradication, while $2$52 for at least one $2$53 yields exponential eradication. Applying this to both viruses guarantees convergence to the healthy state. In the bi-virus case, under $2$54 and $2$55, healing-rate design can enforce

$2$56

entrywise, making $2$57 the only locally asymptotically stable equilibrium; the paper interprets this as using one virus to eradicate the other (Janson et al., 2020).

Beyond epidemiology, the same formal machinery has been mapped to information diffusion, malware spread, and urban transportation. In threshold-based co-diffusion, synergy is encoded through a multivariate Hill kernel and dormancy acts as a one-directional immunity against spreading only. In MBRD-SI and MBRD-CI, superinfection corresponds to competitive replacement, while coinfection parameters $2$58 encode non-interaction, mutual enhancement, one-side enhancement or inhibition, and mutual inhibition (Chang et al., 2018, Yu et al., 21 Aug 2025).

Several limitations are explicit in the literature. The 2019 bi-virus model has no explicit spatial Laplacian. The 2025 co-infection model assumes no diffusion for the $2$59 layer. The multipartite fluid-limit theory fixes the number of islands while taking a large-island limit. The threshold social-contagion model uses algorithmic adoption rules rather than a deterministic closed-form threshold equation. This suggests that “MBRD” is best understood not as a single canonical equation, but as a technically coherent class of multiplex two-contagion systems linked by exclusive or partial occupancy constraints, layer-specific transport, and spectral or mode-based criteria for eradication, dominance, coexistence, and pattern formation [(Liu et al., 2019); (Santos et al., 2013); (Chang et al., 2018); (Yu et al., 21 Aug 2025); (Yu et al., 3 Sep 2025)].