Least-Squares Fitting

One of the two fitting methods available in CRModel.fit() is least-squares estimation (LSE). For Markov Chain Monte Carlo (MCMC), see MCMC Fitting.

Background

To obtain point estimates of the model parameters, we simulate the model and interpolate the resulting species and (optionally) metabolite trajectories to the experimental measurement times using scipy.interpolate.PchipInterpolator. We then compute scaled residuals between model predictions and observations.

The objective function is the sum of squared scaled errors (SSE) across all observed species and, optionally, metabolites,

\[\mathrm{SSE}(\theta) = \sum_i \sum_k \left(\frac{\hat{X}_i(t_k;\theta) - X_i(t_k)}{s_{X,i}}\right)^2 + \sum_j \sum_k \left(\frac{\hat{S}_j(t_k;\theta) - S_j(t_k)}{s_{S,j}}\right)^2\,,\]

where \(\theta\) denotes the vector of model parameters, \(t_k\) are the experimental measurement times, \(X_i(t_k)\) and \(S_j(t_k)\) are the observed concentrations of species \(i\) and resource \(j\), \(\hat{X}_i(t_k;\theta)\) and \(\hat{S}_j(t_k;\theta)\) are the corresponding model predictions, and \(s_{X,i}\) and \(s_{S,j}\) are scaling factors for species and metabolite residuals. Missing observations are ignored in the construction of the residuals.

We minimize this objective using the Trust Region Reflective (TRF) algorithm implemented in scipy.optimize.least_squares, with parameter bounds to enforce biologically meaningful values.

Example

We load experimental time-series data from the BT_WC_export.csv file and fit a consumer-resource model (CRM) to the data.

from mgrowthctrl.utils.data import Dataloader
from mgrowthctrl.models.crm.model import CRModel

# Load dataset
data = Dataloader()
data.load_local_data(
    "examples/datasets/BT_WC_export.csv",
    x_selector=r"BT counts",
)

# Initialize model from data
model = CRModel.from_single_species_data(
    df=data.df,
    time_col=data.time_col,
    x_col=data.X_names[0],
    s_cols=data.S_names,
)

# Fit model parameters
model.fit(
    df=data.df,
    time_col=data.time_col,
    x_cols=data.X_names,
    s_cols=data.S_names,
)

The function CRModel.from_single_species_data() provides a mechanistic initialization of a single-species CRM. It is not recommended using it for multi-species settings.

For a more detailed walkthrough, see Get Started.

Backend API

mgrowthctrl.models.crm.fit.crm_fit_mechanistic( df, time_col: str, x_cols: List[str], s_cols: List[str], simulate_fn: Callable[[Dict[str, ndarray], ndarray, ndarray, ndarray], Dict[str, ndarray]], params0: Dict[str, ndarray] | None = None, lb_factor: float = 0.0, ub_factor: float = 10.0, X0: ndarray | None = None, S0: ndarray | None = None, include_metabolites: bool = True, weight: str = 'median', n_time: int = 300, fit_a: bool = False, ): High-level least-squares fit for CRM parameters.