A Quantitative Systems Pharmacology Model Describing the Cellular Kinetic-Pharmacodynamic Relationship for a Live Biotherapeutic Product to Support Microbiome Drug Development.

Live biotherapeutic products (LBPs) are human microbiome therapies showing promise in the clinic for a range of diseases and conditions. Describing the kinetics and behavior of LBPs poses a unique modeling challenge because, unlike traditional therapies, LBPs can expand, contract, and colonize the host digestive tract. Here, we present a novel cellular kinetic-pharmacodynamic quantitative systems pharmacology model of an LBP. The model describes bacterial growth and competition, vancomycin effects, binding and unbinding to the epithelial surface, and production and clearance of butyrate as a therapeutic metabolite. The model is calibrated and validated to published data from healthy volunteers. Using the model, we simulate the impact of treatment dose, frequency, and duration as well as vancomycin pretreatment on butyrate production. This model enables model-informed drug development and can be used for future microbiome therapies to inform decision making around antibiotic pretreatment, dose selection, loading dose, and dosing duration.

Disease phenotypic and geospatial features vary across genetic lineages for Tuberculosis within Arkansas, 2010-2020.

Tuberculosis (TB) elimination in the United States remains elusive, and community-specific, localized intervention strategies may be necessary to meet elimination goals. A better understanding of the genotypic diversity of Mtb, the population subgroups affected by different TB strains, and differences in disease presentation associated with these strains can aid in identifying risk groups and designing tailored interventions. We analyze TB incidence and genotype data from all Arkansas counties over an 11-year time span from 2010 through 2020. We use statistical methods and geographic information systems (GIS) to identify demographic and disease phenotypic characteristics that are associated with different Mtb genetic lineages in the study area. We found the following variables to be significantly associated with genetic lineage (p<0.05): patient county, patient birth country, patient ethnicity, race, IGRA result, disease site, chest X-ray result, whether or not a case was identified as part of a cluster, patient age, occupation risk, and date arrived in the US. Different Mtb lineages affect different subpopulations in Arkansas. Lineage 4 (EuroAmerican) and Lineage 2 (East Asian) are most prevalent, although the spatial distributions differ substantially, and lineage 2 (East Asian) is more frequently associated with case clusters. The Marshallese remain a particularly high-risk group for TB in Arkansas.

Structural identifiability analysis of age-structured PDE epidemic models.

Computational and mathematical models rely heavily on estimated parameter values for model development. Identifiability analysis determines how well the parameters of a model can be estimated from experimental data. Identifiability analysis is crucial for interpreting and determining confidence in model parameter values and to provide biologically relevant predictions. Structural identifiability analysis, in which one assumes data to be noiseless and arbitrarily fine-grained, has been extensively studied in the context of ordinary differential equation (ODE) models, but has not yet been widely explored for age-structured partial differential equation (PDE) models. These models present additional difficulties due to increased number of variables and partial derivatives as well as the presence of boundary conditions. In this work, we establish a pipeline for structural identifiability analysis of age-structured PDE models using a differential algebra framework and derive identifiability results for specific age-structured models. We use epidemic models to demonstrate this framework because of their wide-spread use in many different diseases and for the corresponding parallel work previously done for ODEs. In our application of the identifiability analysis pipeline, we focus on a Susceptible-Exposed-Infected model for which we compare identifiability results for a PDE and corresponding ODE system and explore effects of age-dependent parameters on identifiability. We also show how practical identifiability analysis can be applied in this example.

To Sobol or not to Sobol? The effects of sampling schemes in systems biology applications.

Computational and mathematical models in biology rely heavily on the parameters that characterize them. However, robust estimates for their values are typically elusive and thus a large parameter space becomes necessary for model study, particularly to make translationally impactful predictions. Sampling schemes exploring parameter spaces for models are used for a variety of purposes in systems biology, including model calibration and sensitivity analysis. Typically, random sampling is used; however, when models have a high number of unknown parameters or the models are highly complex, computational cost becomes an important factor. This issue can be reduced through the use of efficient sampling schemes such as Latin hypercube sampling (LHS) and Sobol sequences. In this work, we compare and contrast three sampling schemes - random sampling, LHS, and Sobol sequences - for the purposes of performing both parameter sensitivity analysis and model calibration. In addition, we apply these analyses to different types of computational and mathematical models of varying complexity: a simple ODE model, a complex ODE model, and an agent-based model. In general, the sampling scheme had little effect when used for calibration efforts, but when applied to sensitivity analyses, Sobol sequences exhibited faster convergence. While the observed benefit to convergence is relatively small, Sobol sequences are computationally less expensive to compute than LHS samples and also have the benefit of being deterministic, which allows for better reproducibility of results.

Predicting the second wave of COVID-19 in Washtenaw County, MI.

The COVID-19 pandemic has highlighted the patchwork nature of disease epidemics, with infection spread dynamics varying wildly across countries and across states within the US. To explore this issue, we study and predict the spread of COVID-19 in Washtenaw County, MI, which is home to University of Michigan and Eastern Michigan University, and in close proximity to Detroit, MI, a major epicenter of the epidemic in Michigan. We apply a discrete and stochastic network-based modeling framework allowing us to track every individual in the county. In this framework, we construct contact networks based on synthetic population datasets specific for Washtenaw County that are derived from US Census datasets. We assign individuals to households, workplaces, schools, and group quarters (such as prisons or long term care facilities). In addition, we assign casual contacts to each individual at random. Using this framework, we explicitly simulate Michigan-specific government-mandated workplace and school closures as well as social distancing measures. We perform sensitivity analyses to identify key model parameters and mechanisms contributing to the observed disease burden in the three months following the first observed cases of COVID-19 in Michigan. We then consider several scenarios for relaxing restrictions and reopening workplaces to predict what actions would be most prudent. In particular, we consider the effects of 1) different timings for reopening, and 2) different levels of workplace vs. casual contact re-engagement. We find that delaying reopening does not reduce the magnitude of the second peak of cases, but only delays it. Reducing levels of casual contact, on the other hand, both delays and lowers the second peak. Through simulations and sensitivity analyses, we explore mechanisms driving the magnitude and timing of a second wave of infections upon re-opening. We find that the most significant factors are workplace and casual contacts and protective measures taken by infected individuals who have sought care. This model can be adapted to other US counties using synthetic population databases and data specific to those regions.

Predicting the second wave of COVID-19 in Washtenaw County, MI.

Marissa Renardy and Denise Kirschner University of Michigan Medical School The COVID-19 pandemic has highlighted the patchwork nature of disease epidemics, with infection spread dynamics varying wildly across countries and across states within the US. These heteroge- neous patterns are also observed within individual states, with patches of concentrated outbreaks. Data is being generated daily at all of these spatial scales, and answers to questions regarded re- opening strategies are desperately needed. Mathematical modeling is useful in exactly these cases, and using modeling at a county scale may be valuable to further predict disease dynamics for the purposes of public health interventions. To explore this issue, we study and predict the spread of COVID-19 in Washtenaw County, MI, the home to University of Michigan, Eastern Michigan University, and Google, as well as serving as a sister city to Detroit, MI where there has been a serious outbreak. Here, we apply a discrete and stochastic network-based modeling framework allowing us to track every individual in the county. In this framework, we construct contact net- works based on synthetic population datasets specific for Washtenaw County that are derived from US Census datasets. We assign individuals to households, workplaces, schools, and group quarters (such as prisons). In addition, we assign casual contacts to each individual at random. Using this framework, we explicitly simulate Michigan-specific government-mandated workplace and school closures as well as social distancing measures. We also perform sensitivity analyses to identify key model parameters and mechanisms contributing to the observed disease burden in the three months following the first observed cases on COVID-19 in Michigan. We then consider several scenarios for relaxing restrictions and reopening workplaces to predict what actions would be most prudent. In particular, we consider the effects of 1) different timings for reopening, and 2) different levels of workplace vs. casual contact re-engagement. Through simulations and sensitivity analyses, we explore mechanisms driving magnitude and timing of a second wave of infections upon re-opening. This model can be adapted to other US counties using synthetic population databases and data specific to those regions.

A Framework for Network-Based Epidemiological Modeling of Tuberculosis Dynamics Using Synthetic Datasets.

We present a framework for discrete network-based modeling of TB epidemiology in US counties using publicly available synthetic datasets. We explore the dynamics of this modeling framework by simulating the hypothetical spread of disease over 2 years resulting from a single active infection in Washtenaw County, MI. We find that for sufficiently large transmission rates that active transmission outweighs reactivation, disease prevalence is sensitive to the contact weight assigned to transmissions between casual contacts (that is, contacts that do not share a household, workplace, school, or group quarter). Workplace and casual contacts contribute most to active disease transmission, while household, school, and group quarter contacts contribute relatively little. Stochastic features of the model result in significant uncertainty in the predicted number of infections over time, leading to challenges in model calibration and interpretation of model-based predictions. Finally, predicted infections were more localized by household location than would be expected if they were randomly distributed. This modeling framework can be refined in later work to study specific county and multi-county TB epidemics in the USA.

Data-Driven Model Validation Across Dimensions.

Data-driven model validation across dimensions in mathematical and computational biology assumptions are often made (e.g., symmetry) to reduce the problem from three spatial dimensions (3D) to two (2D). However, some experimental datasets, such as cell counts obtained via flow cytometry, represent the entire 3D biological object. For purpose of model calibration and validation, it is sometimes necessary to compare these biological datasets with model outputs. We propose a methodology for scaling 2D model outputs to compare with 3D experimental datasets, and we discuss the application of this methodology to two examples: agent-based models of granuloma formation and skeletal muscle tissue. The accuracy of the method is evaluated in artificially generated scenarios.

Evaluating vaccination strategies for tuberculosis in endemic and non-endemic settings.

According to the World Health Organization, tuberculosis (TB) is the leading cause of death from infectious disease worldwide (WHO, 2017). While there is no effective vaccine against adult pulmonary TB, more than a dozen vaccine candidates are in the clinical trial pipeline. These include both pre-exposure vaccines to prevent initial infections and post-exposure vaccines to prevent reactivation of latent disease. Many epidemiological models have been used to study TB, but most have not included a continuous age structure and the possibility of both pre- and post-exposure vaccination. Incorporating age-dependent death rates, disease properties, and social contact data allows for more realistic modeling of disease spread. We propose a continuous age-structured model for the epidemiology of tuberculosis with pre- and post-exposure vaccination. We use uncertainty and sensitivity analysis to make predictions about the efficacy of different vaccination strategies in a non-endemic setting (United States) and an endemic setting (Cambodia). In particular, we determine optimal age groups to target for pre-exposure and post-exposure vaccination in both settings. We find that the optimal age groups tend to be younger for Cambodia than for the US, and that post-exposure vaccination has a significantly larger effect than pre-exposure vaccination in the US.

Global sensitivity analysis of biological multi-scale models.

Mathematical models of biological systems need to both reflect and manage the inherent complexities of biological phenomena. Through their versatility and ability to capture behavior at multiple scales, multi-scale models offer a valuable approach. Due to the typically nonlinear and stochastic nature of multi-scale models as well as unknown parameter values, various types of uncertainty are present; thus, effective assessment and quantification of such uncertainty through sensitivity analysis is important. In this review, we discuss global sensitivity analysis in the context of multi-scale and multi-compartment models and highlight its value in model development and analysis. We present an overview of sensitivity analysis methods, approaches for extending such methods to a multi-scale setting, and examples of how sensitivity analysis can inform model reduction. Through schematics and references to past work, we aim to emphasize the advantages and usefulness of such techniques.

Methods for determining key components in a mathematical model for tumor-immune dynamics in multiple myeloma.

In this work, we analyze a mathematical model we introduced previously for the dynamics of multiple myeloma and the immune system. We focus on four main aspects: (1) obtaining and justifying ranges and values for all parameters in the model; (2) determining a subset of parameters to which the model is most sensitive; (3) determining which parameters in this subset can be uniquely estimated given certain types of data; and (4) exploring the model numerically. Using global sensitivity analysis techniques, we found that the model is most sensitive to certain growth, loss, and efficacy parameters. This analysis provides the foundation for a future application of the model: prediction of optimal combination regimens in patients with multiple myeloma.

Parameter uncertainty quantification using surrogate models applied to a spatial model of yeast mating polarization.

A common challenge in systems biology is quantifying the effects of unknown parameters and estimating parameter values from data. For many systems, this task is computationally intractable due to expensive model evaluations and large numbers of parameters. In this work, we investigate a new method for performing sensitivity analysis and parameter estimation of complex biological models using techniques from uncertainty quantification. The primary advance is a significant improvement in computational efficiency from the replacement of model simulation by evaluation of a polynomial surrogate model. We demonstrate the method on two models of mating in budding yeast: a smaller ODE model of the heterotrimeric G-protein cycle, and a larger spatial model of pheromone-induced cell polarization. A small number of model simulations are used to fit the polynomial surrogates, which are then used to calculate global parameter sensitivities. The surrogate models also allow rapid Bayesian inference of the parameters via Markov chain Monte Carlo (MCMC) by eliminating model simulations at each step. Application to the ODE model shows results consistent with published single-point estimates for the model and data, with the added benefit of calculating the correlations between pairs of parameters. On the larger PDE model, the surrogate models allowed convergence for the distribution of 15 parameters, which otherwise would have been computationally prohibitive using simulations at each MCMC step. We inferred parameter distributions that in certain cases peaked at values different from published values, and showed that a wide range of parameters would permit polarization in the model. Strikingly our results suggested different diffusion constants for active versus inactive Cdc42 to achieve good polarization, which is consistent with experimental observations in another yeast species S. pombe.

Control of cell fraction and population recovery during tissue regeneration in stem cell lineages.

Multicellular tissues are continually turning over, and homeostasis is maintained through regulated proliferation and differentiation of stem cells and progenitors. Following tissue injury, a dramatic increase in cell proliferation is commonly observed, resulting in rapid restoration of tissue size. This regulation is thought to occur via multiple feedback loops acting on cell self-renewal or differentiation. Models of ordinary differential equations have been widely used to study the cell lineage system. Prior modeling studies have suggested that loss of homeostasis and initiation of tumorigenesis can be contributed to the loss of control of these processes, and the rate of symmetric versus asymmetric division of the stem cells may also be altered. While most of the previous works focused on analysis of stability, existence and uniqueness of steady states of multistage cell lineage models, in this work we attempt to understand the cell lineage model from a different perspective. We compare three variants of hierarchical stem cell lineage tissue models with different combinations of negative feedbacks and use sensitivity analysis to examine the possible strategies for the cells to achieve certain performance objectives. Our results suggest that multiple negative feedback loops must be present in the stem cell lineage to keep the fractions of stem cells to differentiated cells in the total population as robust as possible to variations in cell division parameters, and to minimize the time for tissue recovery in a non-oscillatory manner.