seairmobility: Mobility-Based SEAIR Epidemic Models

Description

The seairmobility package provides tools for simulating and analysing mobility-based SEAIR compartmental epidemic models. Each individual in the population carries a fixed mobility trait m \in (0,1) that scales both their susceptibility and their infectiousness through a rank-one transmission kernel. The infectious period is split into an asymptomatic stage with relative infectiousness \alpha and a symptomatic stage with mobility-reduction factor \delta.

Details

The main functions are

seair_params

Build a parameter list.

seair_init

Build initial conditions.

seair_solve

Solve the SEAIR system.

seair_aggregate

Aggregate solution over mobility.

R0_seair

Basic reproduction number.

final_size

Final epidemic size under vanishing seed.

fit_mobility

Parametric least-squares fit of the mobility distribution from an observed symptomatic time series.

Author(s)

Maintainer: Weinan Wang ww@ou.edu

References

Jiang, N., Chu, W., and Li, Y. (2025). Modeling, inference, and prediction in mobility-based compartmental models for epidemiology. SIAM Journal on Applied Mathematics, 85(5), 2355–2375. doi:10.1137/24M1691557.

Basic reproduction number of the mobility-based SEAIR model

Description

For a mobility distribution f on (0,1), the basic reproduction number at the disease-free equilibrium is

\mathcal{R}_0 \;=\; \beta\,c\,\langle m^2 f \rangle,

where c is as in c_effective and \langle m^2 f \rangle = \int_0^1 m^2 f(m)\,dm.

Usage

R0_seair(params, f, m_grid = NULL)

Arguments

Value

The basic reproduction number \mathcal{R}_0.

Examples

pars <- seair_params(beta = 1.5, sigma = 0.3, kappa = 0.2,
                     gamma_A = 0.1, gamma_I = 0.13,
                     alpha = 0.5, delta = 0.3)

## f given as a function (normalised automatically by the caller).
f <- function(m) 6 * m * (1 - m)    # Beta(2, 2) density on [0, 1]
R0_seair(pars, f)

## f given as values on a grid.
m <- seq(0, 1, length.out = 101)
R0_seair(pars, f(m), m_grid = m)

Beta-mixture density on [0, 1]

Description

Evaluates a weighted mixture of Beta densities

f(m) \;=\; \sum_{k=1}^{K} w_k\,\mathrm{Beta}(m; a_k, b_k),

renormalised on [0,1]. Weights are coerced to be nonnegative and to sum to one.

Usage

beta_mixture_density(m, weights, shape1, shape2)

Arguments

Value

Numeric vector of densities at m.

Examples

m <- seq(0, 1, length.out = 101)
beta_mixture_density(m,
                     weights = c(0.5, 0.3, 0.2),
                     shape1  = c(2, 5, 9),
                     shape2  = c(8, 3, 2))

Effective infectiousness duration of the SEAIR chain

Description

Computes the scalar

c \;=\; \frac{\alpha}{\kappa + \gamma_A} \;+\; \frac{\delta\,\kappa}{(\kappa + \gamma_A)\,\gamma_I},

the integrated contribution of a unit of newly infected mass to the force of infection, collapsing the E \to A \to I stage chain into a single number. The latency rate \sigma does not appear in c.

Usage

c_effective(params)

Arguments

Value

A positive scalar.

Examples

pars <- seair_params(beta = 1.5, sigma = 0.3, kappa = 0.2,
                     gamma_A = 0.1, gamma_I = 0.13,
                     alpha = 0.5, delta = 0.3)
c_effective(pars)

Final epidemic size under the vanishing-seed regime

Description

Solves the scalar fixed-point equation

M_R^\infty \;=\; \int_0^1 m f(m)\,(1 - e^{-\beta c M_R^\infty m})\,dm

for the first mobility moment of the recovered compartment at the end of the epidemic (with I_0 \equiv 0, R_0 \equiv 0, S_0 = f), and returns the total final epidemic size

R_\infty \;=\; 1 - \int_0^1 f(m)\,e^{-\beta c M_R^\infty m}\,dm.

When the basic reproduction number is less than or equal to one, the function returns R_\infty = 0.

Usage

final_size(params, f, m_grid = NULL, tol = 1e-10)

Arguments

Value

A list with R_inf (total final size), MR_inf (first moment), R0 (basic reproduction number), and converged (logical).

Examples

pars <- seair_params(beta = 1.5, sigma = 0.3, kappa = 0.2,
                     gamma_A = 0.1, gamma_I = 0.13,
                     alpha = 0.5, delta = 0.3)
f <- function(m) 6 * m * (1 - m)
final_size(pars, f)

Final epidemic size under general proportional initial conditions

Description

Solves the scalar equation of Theorem 4.3 (final-size relation),

M_R^\infty + \int_0^1 m S_0(m)\, \exp\!\bigl(-\beta m [ c\,(M_R^\infty - M_R^0) + (\delta/\gamma_I - c)\,M_I^0 ]\bigr) dm \;=\; \int_0^1 m f(m)\,dm,

allowing a nonzero initial symptomatic fraction I_0 = I_{\mathrm{seed}} f, with E_0 = A_0 = R_0 = 0.

Usage

final_size_general(params, f, m_grid = NULL, I_seed = 1e-04, tol = 1e-10)

Arguments

Value

A list with R_inf, MR_inf, R0, and converged.

Examples

pars <- seair_params(beta = 1.5, sigma = 0.3, kappa = 0.2,
                     gamma_A = 0.1, gamma_I = 0.13,
                     alpha = 0.5, delta = 0.3)
f <- function(m) 6 * m * (1 - m)
final_size_general(pars, f, I_seed = 0.001)

Fit a mobility distribution to an observed symptomatic time series

Description

Fits a Beta-mixture mobility distribution f by minimising the sum of squared errors between an observed aggregate symptomatic time series \langle I \rangle(t) and the simulated output of the mobility-based SEAIR model at the same times. All transmission and stage parameters are held fixed.

Usage

fit_mobility(
  times,
  I_obs,
  params,
  K = 3,
  m_grid = seq(0, 1, length.out = 51),
  I_seed = 1e-04,
  start = NULL,
  n_restarts = 5,
  control = list(maxit = 500),
  verbose = FALSE
)

Arguments

Details

The unconstrained optimisation parameter is, for each component k = 1, \ldots, K, a triple (\log a_k, \log b_k, \log \tilde{w}_k) with mixture weights recovered as w_k = \tilde{w}_k / \sum_j \tilde{w}_j. The first weight's log is fixed to 0 to remove the scale ambiguity, so the parameter vector has length 3 K - 1.

Value

A list with weights, shape1, shape2 (the fitted Beta-mixture parameters), f_hat (density on m_grid), m_grid, loss (final SSE), I_fit (simulated aggregate symptomatic time series at times), and optim_result.

Examples

pars <- seair_params(beta = 1.5, sigma = 0.3, kappa = 0.2,
                     gamma_A = 0.1, gamma_I = 0.13,
                     alpha = 0.5, delta = 0.3)
m <- seq(0, 1, length.out = 21)
f_true <- dbeta(m, 3, 3)
init <- seair_init(m, f_true, I_seed = 1e-4)
times <- seq(0, 20, by = 5)
sol <- seair_solve(init, pars, times)
I_obs <- seair_aggregate(sol)$I

fit <- fit_mobility(times, I_obs, pars,
                    K = 1, m_grid = m,
                    start = log(c(3, 3)),
                    n_restarts = 0,
                    control = list(maxit = 20))
fit$loss

Plot aggregate compartment trajectories

Description

Convenience plot of the five aggregate compartment trajectories produced by a SEAIR solution.

Usage

plot_seair(sol, which = c("S", "E", "A", "I", "R"), ...)

Arguments

Value

Invisibly returns the aggregate data frame.

Examples

m <- seq(0, 1, length.out = 31)
f <- dbeta(m, 2, 2)
init <- seair_init(m, f, I_seed = 1e-4)
pars <- seair_params(beta = 1.5, sigma = 0.3, kappa = 0.2,
                     gamma_A = 0.1, gamma_I = 0.13,
                     alpha = 0.5, delta = 0.3)
sol <- seair_solve(init, pars, times = seq(0, 30, by = 2))
plot_seair(sol, which = c("S","I","R"))

Aggregate a SEAIR solution over mobility

Description

Computes the aggregate time series \int_0^1 S(m,t)\,dm, \int_0^1 E(m,t)\,dm, etc., using the trapezoidal rule on the solver's mobility grid.

Usage

seair_aggregate(sol)

Arguments

Value

A data frame with columns time, S, E, A, I, R (aggregate ratios).

Examples

m <- seq(0, 1, length.out = 31)
f <- dbeta(m, 2, 2)
init <- seair_init(m, f, I_seed = 1e-4)
pars <- seair_params(beta = 1.5, sigma = 0.3, kappa = 0.2,
                     gamma_A = 0.1, gamma_I = 0.13,
                     alpha = 0.5, delta = 0.3)
sol <- seair_solve(init, pars, times = seq(0, 30, by = 2))
agg <- seair_aggregate(sol)
head(agg)

Build initial conditions for the SEAIR system

Description

Constructs initial profiles (S_0, E_0, A_0, I_0, R_0) on the mobility grid from a user-supplied mobility density f(m) and a seed fraction I_{\mathrm{seed}} placed in the symptomatic compartment. The initial condition is "proportional": S_0 = (1 - I_{\mathrm{seed}}) f, I_0 = I_{\mathrm{seed}} f, E_0 = A_0 = R_0 = 0.

Usage

seair_init(m_grid, f_vals, I_seed = 1e-04)

Arguments

Value

A list of class "seair_init" with entries m_grid, f, S, E, A, I, R, and I_seed.

Examples

m <- seq(0, 1, length.out = 101)
f <- dbeta(m, 2, 2)
seair_init(m, f, I_seed = 1e-4)

Build a parameter list for the mobility-based SEAIR model

Description

Constructs and validates the list of transmission and stage-progression parameters used by the mobility-based SEAIR model \partial_t S = -\beta m S \mathcal{K}, \partial_t E = \beta m S \mathcal{K} - \sigma E, \partial_t A = \sigma E - (\kappa + \gamma_A) A, \partial_t I = \kappa A - \gamma_I I, \partial_t R = \gamma_A A + \gamma_I I, where \mathcal{K}(A,I) = \int_0^1 \bar{m}\,(\alpha A(\bar{m}) + \delta I(\bar{m}))\,d\bar{m}.

Usage

seair_params(beta, sigma, kappa, gamma_A, gamma_I, alpha = 1, delta = 1)

Arguments

Value

A list of class "seair_params" containing the checked parameters.

Examples

seair_params(beta = 1.5, sigma = 0.3, kappa = 0.2,
             gamma_A = 0.1, gamma_I = 0.13,
             alpha = 0.5, delta = 0.3)

Solve the mobility-based SEAIR system

Description

Discretises the mobility-based SEAIR system on the user-supplied mobility grid (using the trapezoidal rule for the non-local integral \mathcal{K}) and integrates the resulting coupled ODE system over times with deSolve.

Usage

seair_solve(init, params, times, method = "lsoda", ...)

Arguments

Value

A list of class "seair_solution" with entries times, m_grid, S, E, A, I, R (each a matrix of dimension length(times) by length(m_grid)), together with the input params and init.

Examples

m <- seq(0, 1, length.out = 31)
f <- dbeta(m, 2, 2)
init <- seair_init(m, f, I_seed = 1e-4)
pars <- seair_params(beta = 1.5, sigma = 0.3, kappa = 0.2,
                     gamma_A = 0.1, gamma_I = 0.13,
                     alpha = 0.5, delta = 0.3)
sol <- seair_solve(init, pars, times = seq(0, 30, by = 2))