Numerical Methods

1 ODE Integration

The stiff ODE systems arising from the coupled multi-system model (42 state variables) are integrated using implicit methods. The recommended solver is the Radau IIA method (5th-order, L-stable), which handles the wide range of timescales present in the model (milliseconds for enzyme kinetics to hours for cytokine dynamics) without requiring manual timescale separation. Absolute and relative error tolerances are set to $ 10^{-8}$ and $ 10^{-6}$ respectively.

For the stochastic viral reactivation model, the Gillespie stochastic simulation algorithm (SSA) is used for the latent-to-lytic transition, coupled to deterministic ODE integration for the subsequent immune response. The hybrid scheme switches from SSA to ODE at the moment of reactivation, using the SSA-generated viral load as an initial condition for the ODE system.

2 Steady-State Computation

Steady states of the ODE system satisfy \(\mathbf{f}(\mathbf{x}^*, \mathbf{\theta}) = \mathbf{0}\). These are computed using Newton’s method with analytical Jacobian. Stability is determined by the eigenvalues of the Jacobian evaluated at \(\mathbf{x}^*\): the steady state is stable if all eigenvalues have negative real part, and unstable otherwise. Bifurcation analysis (continuation of steady states as parameters vary) is performed using pseudo-arclength continuation to trace solution branches through turning points and detect Hopf bifurcations.