Superinfection Model MBRD-SI

• The paper introduces MBRD-SI as a modeling principle that captures mechanisms where later infection events overlap, exclude, or displace earlier ones.
• The methodology integrates stochastic-population processes, generating functions, and reaction–diffusion systems to derive thresholds, extinction criteria, and pattern formations.
• Applications span from modeling hypnozoite accrual in Plasmodium vivax to competitive exclusion in Hib and network-based pathogen displacement for public health insights.

The Super-Infection Model (MBRD-SI) denotes, in the cited literature, a set of mathematically distinct frameworks for representing infection processes in which later infection events overlap with, exclude, or displace earlier ones. The designation has been used for a general stochastic-population model of Plasmodium vivax with explicit hypnozoite accrual and blood-stage superinfection, for a multitype branching-process approximation of its early invasion regime, for a mixed stochastic–deterministic model of single-variant bottlenecks in Haemophilus influenzae bacteremia, and for multiplex bi-virus reaction–diffusion systems in which pathogen 2 can super-infect hosts already infected by pathogen 1 (Mehra et al., 2023, Mehra et al., 24 Feb 2025, Shao et al., 2016, Yu et al., 21 Aug 2025, Yu et al., 3 Sep 2025).

1. Scope of the designation and common structure

Across these formulations, “superinfection” does not denote a single universal mechanism. In the P. vivax setting, it refers to multiple concurrent blood-stage broods produced by temporally proximate mosquito inoculations and/or hypnozoite activation events. In the Hib bottleneck model, it refers to exclusion mediated by a race between phenotypic switching and innate immune recruitment. In the multiplex network formulation, it denotes one-sided hierarchical displacement: pathogen 2 can infect hosts already infected by pathogen 1, but simultaneous co-infection is not allowed (Mehra et al., 2023, Shao et al., 2016, Yu et al., 21 Aug 2025).

Framework	Core state description	Superinfection mechanism
P. vivax host–vector model	humans indexed by hypnozoite load $i$ and multiplicity of broods $j$; mosquitoes by infection status	infective bites add one primary brood and a geometric hypnozoite batch; activation adds relapse broods
Hib bottleneck model	strains A and B in “crossing” $C$ and “growing” $G$ phenotypes, plus innate immunity $I$	the first strain to switch to $G$ recruits immunity that suppresses the competitor
Multiplex bi-virus reaction–diffusion model	node-wise $S_i$, $I_i$, $J_i$ on three layers	pathogen 2 displaces pathogen-1-infected hosts through a takeover term proportional to $\sigma \beta_2$

A recurrent structural feature is the use of hybrid state spaces that separate mechanistically distinct infection stages. The malaria models resolve latent liver-stage parasites, circulating broods, and vector infection; the Hib model resolves barrier crossing, rapid bloodstream growth, and non-clonal innate immune recruitment; the multiplex models separate susceptible, pathogen-1-infected, and pathogen-2-infected populations across distinct diffusion layers. This suggests that “MBRD-SI” functions less as a single canonical equation set than as a model-class label for mechanistic superinfection dynamics.

2. Host–vector superinfection in Plasmodium vivax

The general P. vivax framework models the coupled human–mosquito system as a density-dependent Markov population process with countably many types, built on a within-host open network of infinite-server queues with inhomogeneous batch arrivals (Mehra et al., 2023). Within a single human host, hypnozoites occupy four states,

$j$0

where $j$1 denotes a hypnozoite currently present in the liver, $j$2 an activated hypnozoite causing a current relapse, $j$3 a previously activated hypnozoite whose relapse has cleared, and $j$4 a hypnozoite that died without activating. Independently, primary infections occupy states $j$5 and $j$6.

For a single hypnozoite inoculated at time $j$7, the continuous-time Markov chain has rates

$j$8

with absorbing states $j$9 and $C$0. When $C$1, the state probabilities at time $C$2 are

$C$3

Each infective mosquito bite is a Poisson event at time-dependent rate $C$4. It immediately creates one primary blood-stage infection and deposits a geometrically distributed batch $C$5 of hypnozoites with mean $C$6,

$C$7

Hypnozoites in $C$8 activate independently at rate $C$9, die at rate $G$0, and activated relapses clear at rate $G$1. Blood-stage infections also clear at rate $G$2. Seasonality enters through non-homogeneous Poisson arrivals, or equivalently through time-varying transmission parameters such as $G$3, $G$4, and $G$5.

At population scale, humans are indexed by hypnozoite load $G$6 and multiplicity of broods $G$7. If $G$8 is the number of humans with reservoir size $G$9 and brood multiplicity $I$0, then brood clearance, hypnozoite death, hypnozoite activation, mosquito-to-human infection, mosquito death/replacement, and human-to-mosquito infection all appear as explicit transitions of the countably typed Markov process. Rescaling by mosquito population size yields a functional law of large numbers limit in the form of an infinite-dimensional semilinear ODE system for the host fractions $I$1 and infected mosquito prevalence $I$2. The host component coincides exactly with the Kolmogorov forward equations of a Markovian network of infinite-server queues under inhomogeneous batch arrivals.

3. Generating functions, threshold structure, and branching extinction theory

A central analytic step in the P. vivax formulation is the generating-function characterization of the host-state distribution. With $I$3, the host generating function

$I$4

is

$I$5

which encodes the full distribution over hypnozoite loads and multiplicity of broods (Mehra et al., 2023). Setting $I$6 and $I$7 gives human blood-stage prevalence,

$I$8

Substituting this into the mosquito equation collapses the countably infinite host system to a single integrodifferential equation for $I$9, or equivalently for the force of reinfection $G$0, with kernel

$G$1

The steady-state analysis yields the threshold quantity

$G$2

where $G$3. The disease-free equilibrium is uniformly asymptotically stable if $G$4, whereas an endemic equilibrium emerges if $G$5. The endemic mosquito prevalence $G$6 satisfies

$G$7

and a sufficient condition for uniform asymptotic stability of that endemic equilibrium is also given in closed form.

The multitype branching-process approximation reformulates the early epidemic phase by tracking only infected entities: human particles $G$8 with $G$9 and infected mosquitoes $S_i$0 (Mehra et al., 24 Feb 2025). It preserves within-inoculum superinfection, because a single infected bite creates one primary brood plus a geometric hypnozoite batch, but it ignores depletion of susceptible mosquitoes and superinfection generated by multiple inocula in the same host. The extinction criterion is expressed through the mean progeny operator $S_i$1: global extinction occurs with probability $S_i$2 if and only if $S_i$3. In this formulation,

$S_i$4

with

$S_i$5

and $S_i$6. The coupling with the full epidemic process yields a total variation bound of order $S_i$7, valid until $S_i$8 human-to-mosquito and mosquito-to-human transmission events have occurred, for arbitrary $S_i$9. In the applications developed there, this approximation is used to quantify re-introduction risk after prior elimination and elimination probability after mass drug administration with hypnozoite depletion.

4. Superinfection exclusion and the Hib single-variant bottleneck

In the Hib model, superinfection is formulated as competitive exclusion in the bloodstream rather than coexistence of multiple broods (Shao et al., 2016). Neonatal rats were intranasally co-inoculated with equal mixtures of streptomycin-sensitive and streptomycin-resistant Haemophilus influenzae type b strains. Both variants were detectable in blood shortly after inoculation, but by $I_i$0 hours most infections were “pure,” with either variant dominating and mixed infections uncommon.

For each strain, bacteria occupy two phenotypic states. The $I_i$1 (“crossing”) phenotype can traverse the nasopharynx-to-blood barrier but does not grow well in blood; the $I_i$2 (“growing”) phenotype grows rapidly in blood but cannot cross the barrier. With strain A for Sm$I_i$3 and strain B for Sm$I_i$4, the state variables are $I_i$5, $I_i$6, $I_i$7, $I_i$8, and $I_i$9, where $J_i$0 is the fraction of an innate immune-cell reservoir recruited to the infection site. The dynamics are

$J_i$1

$J_i$2

with symmetric equations for $J_i$3 and $J_i$4, and

$J_i$5

Throughout, the simplifying choices are $J_i$6, $J_i$7, $J_i$8, and $J_i$9.

The baseline crossing flux is deterministic and proportional to inoculum size: $\sigma \beta_2$0 Switching $\sigma \beta_2$1 is stochastic. In the single-step model, the number of switches during $\sigma \beta_2$2 is Poisson with hazard $\sigma \beta_2$3 for strain A and $\sigma \beta_2$4 for strain B. In the multi-step extension, the transition proceeds through $\sigma \beta_2$5, each step at rate $\sigma \beta_2$6, and only the final state grows rapidly.

The exclusion mechanism is explicit: the first strain to produce a $\sigma \beta_2$7-phenotype clone grows quickly, accelerates innate immune recruitment, and thereby increases the killing term $\sigma \beta_2$8 for both strains. The competitor, if it has not yet switched or switches later, is suppressed or cleared. For the independent-action baseline, the mean time to the first switch scales as

$\sigma \beta_2$9

which makes the model highly sensitive to inoculum size. In the multi-step case,

$j$00

weakening the dependence on $j$01 but also reducing switching-time variance.

The principal quantitative result is negative: independent-action models, with $j$02, could not reproduce the weak dependence of observed outcome frequencies on inoculum size. Mixed infections require both strains at least $j$03 cells/ml; pure infections require at least $j$04 cells/ml of one strain and fewer than $j$05 cells/ml of the other. Under those criteria, the best improvements required non-independent action. Two fitted modifications were studied: a constant crossing magnitude with sublinear duration,

$j$06

and a sublinear crossing magnitude,

$j$07

The latter gave the best fit, with $j$08, implying that effective per-bacterium crossing probability decreases with inoculum size. The model therefore directly challenges the theory of independent action in this experimental system.

5. Multiplex reaction–diffusion MBRD-SI on networks

The network MBRD-SI formulation is a three-layer multiplex bi-virus reaction–diffusion system in which each node $j$09 carries susceptible density $j$10, pathogen-1 infection density $j$11, and pathogen-2 infection density $j$12, with no co-infected compartment (Yu et al., 21 Aug 2025). Superinfection is one-sided: pathogen 2 can infect hosts already infected by pathogen 1 and displace it. With $j$13, the node-level reaction system is

$j$14

On the multiplex, each layer has its own Laplacian $j$15, defined by

$j$16

and the full coupled system is

$j$17

Here $j$18 and $j$19 are cross-diffusion coefficients: when they are positive, susceptibles move toward lower $j$20 and $j$21 concentrations; negative values model attraction to higher infection densities.

For homogeneous equilibria, the disease-free state satisfies

$j$22

At this equilibrium, the basic reproduction numbers are

$j$23

If pathogen 1 is endemic and pathogen 2 is rare, pathogen 2 invades when

$j$24

showing explicitly that $j$25 lowers the invasion barrier for pathogen 2. The corresponding invasion condition for pathogen 1 into a pathogen-2 endemic state contains the negative displacement term

$j$26

so sufficiently large $j$27 favors replacement equilibria.

The linear instability analysis is formulated through the reaction Jacobian at a homogeneous equilibrium and the cubic characteristic polynomial

$j$28

With

$j$29

$j$30

local reaction stability requires

$j$31

Spatial instability then depends on the cubic discriminant $j$32 and on the signs of $j$33 and $j$34 for nonzero modes. The paper distinguishes Turing instability and Turing–Hopf instability through these criteria, and in the single-layer special case reduces the dispersion relation to eigenmodes of a common Laplacian. The stated interpretation is that cross-diffusion and multiplexity broaden the parameter window for instability by decoupling movement scales of $j$35, $j$36, and $j$37.

6. Numerical behavior, applications, and interpretive issues

The later dynamics paper emphasizes numerical pattern formation and hotspot growth for the network MBRD-SI, using lattices such as LA4, LA12, and LA24, as well as Watts–Strogatz and Barabási–Albert networks (Yu et al., 3 Sep 2025). A pattern-amplitude metric is introduced as

$j$38

For the parameter set

$j$39

$j$40

stationary Turing patterns appear on LA12 lattices, with stationary spotted clusters reported in the $j$41-layer at $j$42. On WS networks with the same average degree, similar pattern shapes are reported. Quantitatively, the early-time amplitude at $j$43 is fitted by

$j$44

as a function of $j$45, indicating rapid suppression of pattern amplitude as superinfection strength increases. In point-source experiments, higher $j$46 reduces pathogen-1 peak prominence and accelerates pathogen-2 saturation; BA networks generally exhibit faster diffusion and earlier pathogen-2 takeover than WS networks.

Implementation guidance is explicit in both the malaria and multiplex settings. For the infinite-dimensional malaria ODEs, truncation at finite $j$47 and $j$48, stiff solvers, positivity preservation, conservation of $j$49, and adaptive expansion of the truncation domain are recommended (Mehra et al., 2023). For the reduced IDE, convolution quadrature and efficient history arrays are suggested, with equivalence to the full collapse relying on the initial condition $j$50. In the multiplex system, Laplacians are precomputed from layer adjacencies, explicit schemes such as Euler or Runge–Kutta are used, $j$51 is chosen below a CFL-like bound tied to the largest diffusion coefficient and Laplacian spectral radius, and non-negativity may be enforced by clamping (Yu et al., 21 Aug 2025).

The application range is correspondingly broad. In malaria, the branching approximation is used for re-introduction after prior elimination and for idealized mass drug administration with radical cure, where elimination probability depends sharply on entomological parameters, blood-stage clearance, and hypnozoite survival probability $j$52 (Mehra et al., 24 Feb 2025). In the multiplex reaction–diffusion setting, the same formalism is proposed for information propagation, malware diffusion, and urban transportation networks, with transmission, recovery, virulence, takeover, diffusion, and cross-diffusion reinterpreted in those domains (Yu et al., 21 Aug 2025).

Several interpretive cautions follow directly from the models. First, MBRD-SI is not synonymous with co-infection: the multiplex superinfection model explicitly has no co-infected compartment, whereas the companion MBRD-CI model does. Second, in the Hib study, the best-fitting explanations violate independent action through sublinear crossing or saturation. Third, in the P. vivax setting, the framework is designed precisely to avoid a binary treatment of hypnozoite carriage and the pseudoequilibrium approximation for recovery rates; relapse risk depends on hypnozoite burden and multiplicity of broods is unbounded (Mehra et al., 2023, Shao et al., 2016, Yu et al., 21 Aug 2025). These contrasts indicate that the central contribution of MBRD-SI lies not in a single canonical set of equations, but in a transferable modeling principle: superinfection must be represented at the mechanistic scale where overlap, takeover, or exclusion is actually generated.