Parametric Sensitivity Analysis for Models of Reaction Networks within Interacting Compartments

TL;DR

This work tackles parametric sensitivity analysis for reaction-network-in-interacting-compartments (RNIC), a generalization of stochastic reaction networks with dynamic, interacting cells. It systematically adapts unbiased Girsanov/Likelihood Ratio and multiple finite-difference coupling methods (CRV, CRP, local-CRP, Split) to RNIC, addressing when parameters lie in the chemical network versus the compartment process. The study demonstrates, via numerical examples, that Girsanov remains unbiased and intuitive to use, while coupling-based finite differences offer substantial variance reduction, with Split Coupling excelling for few parameter perturbations and CRP variants becoming more efficient as the parameter count grows. The results provide actionable guidance for practitioners and include freely available Matlab implementations, enabling efficient sensitivity analysis of RNIC models in biological and chemical contexts.

Abstract

Models of reaction networks within interacting compartments (RNIC) are a generalization of stochastic reaction networks. It is most natural to think of the interacting compartments as "cells" that can appear, degrade, split, and even merge, with each cell containing an evolving copy of the underlying stochastic reaction network. Such models have a number of parameters, including those associated with the internal chemical model and those associated with the compartment interactions, and it is natural to want efficient computational methods for the numerical estimation of sensitivities of model statistics with respect to these parameters. Motivated by the extensive work on computational methods for parametric sensitivity analysis in the context of stochastic reaction networks over the past few decades, we provide a number of methods in the basic RNIC setting. Provided methods include the (unbiased) Girsanov transformation method (also called the Likelihood Ratio method) and a number of coupling methods for the implementation of finite differences, each motivated by methods from previous work related to stochastic reaction networks. We provide several numerical examples comparing the various methods in the new setting. We find that the relative performance of each method is in line with its analog in the "standard" stochastic reaction network setting. We have made all of the Matlab code used to implement the various methods freely available for download.

Figures (6)

• Figure 1: Estimator variances \ref{['eq:estimatorvarianceformula']} for $\frac{d}{d\kappa_b} \mathbb E[ S_{tot}(t)]$ using various methods. For the (forward) finite difference methods, we used a perturbation of $h = 1/10$ and $n=$ 1,000 coupled paths. For the Girsanov method we used $n =$ 2000 paths. Note the vastly different scales on the $y$-axis.
• Figure 2: Estimator variances \ref{['eq:estimatorvarianceformula']} for $\frac{d}{d\kappa_E} \mathbb E[S_{tot}(t)]$ using various methods. For the (centered) finite difference methods, we used a perturbation of $h = 1/1000$ and $n=$ 1,000 coupled paths. For the Girsanov method we used $n =$ 2000 paths. Note that the Split Coupling, with both the compartment model and chemistry model being coupled, as detailed at the end of Section \ref{['sec:FD']}, has by far the lowest variance and is also the fastest (and utilizes the fewest random variables).
• Figure 3: A realization of the compartment model for the birth and death process of Example \ref{['example:gene']}.
• Figure 4: Estimator variances \ref{['eq:estimatorvarianceformula']} for $\frac{d}{d\kappa_2} \mathbb E[ \text{Total dimers}]$ using various methods. For the (centered) finite difference methods, we used a perturbation of $h = 1/10$ and $n =$ 10,000 coupled paths. For the Girsanov method we used $n=$ 20,000 paths. We leave out the local-CRP method since it is indistinguishable from the Split Coupling and the $x$-axis.
• Figure 5: Estimator variance \ref{['eq:estimatorvarianceformula']} for $\frac{d}{d\kappa_2} \mathbb E[ \text{Total dimers}]$ using various methods. Each method shares the compartment model. However, the chemical models are coupled in various ways including: Common Random Variables, local-CRP, and the Split Coupling. Note that this figure simply zooms in on the plots for Common Random Variables and the Split Coupling from Figure \ref{['fig:2345670']}, while adding the plot for the local-CRP method.

Theorems & Definitions (6)

• Definition 2.1
• Example 2.1
• Remark 1
• Remark 2
• Example 4.1
• Example 4.2: Genetic model with coagulation