OneHealthStudio platform for mechanistic systems modeling - showcase PBPK/PD of siRNA

1 Introduction

We introduce OneHealthStudio, a platform designed to provide scalable and reliable mechanistic computational modeling, simulation, and analytics through scientific machine learning.

As the mechanistic models take physiology into account, ordinary differential equations are often utilized to represent the pathophysiological mechanisms. These ODE models can include over 100–1000 state variables and exhibit significant nonlinearity, which further challenges parameter estimation routine to ensure model credibility (Kapfer et al., 2019).

Here, we demonstrate how it integrates with PBPK/PD models built using OSPSuite (Eissing et al., 2011), and its merits in modeling and simulation in a few large-scale models. In the workflow, users first create PBPK/PD/QSP models using graphical user interface tools such as PK-Sim, Mobi as well as other general modeling standard like systems biology markup language (SBML). The initial project or model can then be saved in a shared format such as pkml. After that, we parsed and extracted the set of ordinary differential equations (ODE), which is encoded in pkml format.

As a scalable solution for large-scale ODE model, OneHealthStudio utilizes high-performance ODE solvers from SUNDIALS/CVODES (Serban & Hindmarsh, 2005). We benchmarked simulation performance of OneHealthStudio with OSPSuite-R and httk. In addition, we enhanced OneHealthStudio estimation engine with the scientific machine learning approach like automatic differentiation (AD). AD ensures more accurate sensitivity and Jacobians, which have been shown to yield more precise gradients in gradient-based optimization of systems biology models (Raue et al., 2013). We evaluated the performance of AD approach in our OneHealthStudio against finite-difference (FD) methods, which are currently implemented in parameter identification in neither PK-Sim/Mobi (Eissing et al., 2011) nor the R package OSPSuite-PI.

The highlights of OneHealthStudio are

• 7-10 fold faster simulation than OSPSuite-R (demonstrated in our case study of PBPK/PD modeling of siRNA treatment),
• 100 fold faster simulation than httk from US-EPA (see our separate note),
• less iterations in SciML-enabled parameter estimation and more reliable and accurate sensitivity,
• profile likelihood approach can reveal more realistic parameter uncertainty than conventional Fisher information matrix.

In the following, we exemplified OneHealthStudio in a case study of PBPK/PD model of siRNA.

2 Case study: PBPK/PD modeling of siRNA

We adapted a recently published PBPK/PD model of N-Acetylgalactosamine (GalNAc)-Conjugated siRNAs by Salim et al. (2025). In our study, we specifically took the PBPK model of ALN-AT3 (fitusiran) with subcutaneous injection to mice. As a note to size of ODEs, this PBPK model was mathematically represented by a set of ordinary differential equations with 434 states. Our OneHealthStudio extracted, parsed and annotated the entire set of ODEs.

In evaluating and benchmarking OneHealthStudio, we modified two aspects: first assuming dosing levels at 1, 2, and 5 [mg/kg] and then conducting parameter estimation. In parameter estimation, we first generated the data based on the reported parameter values in the original publication (Salim et al., 2025), see Table 1.

In our simulation, we took three types of synthetic observable:

• RNA-induced silencing complex(RISC)-siRNA: concentration of RISC-siRNA at container Liver|Periportal|Intracellular,
• mRNA: pharmacological biomarker,
• protein: pharmacological biomarker.

Then in our estimation, we applied multiple-starting value optimization by using 20 Latin hypercube sampling (Loh, 1996) in a multidimensional space (\(k_{deg}, k_{DR}, SC_{50}, S_{max}\)). We performed optimization using a trust region reflective (TRF) approach (Branch et al., 1999).

library(knitr)

# Create sample data
set.seed(789)
data <- data.frame(
  parameter = c("kdeg", "kDR", "SC50", "Smax"),
  description = c("Degradation rate constant of ASGPR on hepatocyte", "Degradation rate constant of RISC complex", "Maximum stimulation of mRNA degradation", "RISC-loaded siRNA at half-maximal stimulation"),
  value = c(1.52, 0.0033, 13.1, 3.52),
  unit = c("1/h", "1/h", "-", "nmol/L"),
  range_parameter = c("(0, inf)", "(0, inf)", "(0, inf)", "(0, inf)")
  #range_parameter_profiling = c("(, )", "(, )", "(, )", "(, )")
)


# Create the table
kable(data)

Table 1: Parameters, their description and their settting in estimation percedure

parameter	description	value	unit	range_parameter
kdeg	Degradation rate constant of ASGPR on hepatocyte	1.5200	1/h	(0, inf)
kDR	Degradation rate constant of RISC complex	0.0033	1/h	(0, inf)
SC50	Maximum stimulation of mRNA degradation	13.1000	-	(0, inf)
Smax	RISC-loaded siRNA at half-maximal stimulation	3.5200	nmol/L	(0, inf)

2.1 Evaluation metrics

We evaluated the estimation routine using automatic differentiation approach (Bücker, 2006). Jointly, we checked three following metrics,

• number of iterations until the same optimization criteria met (AD-approach vs FD-approach) using TRF optimization;
• sensitivity accuracy (OneHealthStudio-AD vs OSPSuite-R);
• simulation runtime (OneHealthStudio vs OSPSuite-R).

2.2 Model simulation

We benchmarked model simulation with 1 [mg/kg] siRNA of ALN-AT3 (fitusiran). We further provided an interactive interface, which further enable a repeated dosing regimen of siRNA, see onehealthstudio.shinyapps.io

We compared the runtime and checked simulation alignment between the two simulators OneHealthStudio and OSPSuite-R with the same tolerance in the ODE solver (AbsTol=1e-10 and RelTol=1e-10). Figure 1 showed a 8-fold speedup of OneHealthStudio compared to OSPSuiteR.

Figure 1

2.3 Parameter estimation

OneHealthStudio employed automatic differentiation to ensure accurate sensitivity and robust in estimation as well as less iteration needed.

We compared the sensitivity term from automatic differentiation and finite difference, where we exemplified the sensitivity \(\frac{\partial [RISC]}{\partial \theta}\). In FD approach, we used \(\delta_{\theta} = 1e-6\). For parameters \(k_{deg}\) and \(k_{DR}\), FD approximation worked well as AD approach. However, for parameters \(SC_{50}\) and \(S_{max}\), the numerical error in the FD approach can be reflected into a non-zeros derivative, see subplots related to \(SC_{50}\) and \(S_{max}\) in Figure 2. Note that in this PBPK/PD model, both parameters only impact PD markers mRNA and protein but not PK part, given no feedback from PD back to PK.

Figure 2: time course sensitivities of [RISC-siRNA] w.r.t. four parameters in PBPKPD model

Hence, automatic differentiation approach in OneHealthStudio show more robust sensitivity analysis than conventional finite difference approach. Furthermore, in another PBPK modeling for a small molecule, Appendix recapped comparison of AD and FD approach in computing sensitivities. We compared the number of iteration in either AD or FD approaches, see Figure 3.

Figure 3: OneHealthStudio employing AD in optimization (20 average iterations) need ~5 times less iterations in estimation procedure compared to conventional optimization with the finite difference approach (109 average iterations).

2.4 Parameter uncertainty with Fisher and profile likelihood approaches

We profiled the parameters and checked \(\chi^2\)-based threshold to extract 95% confidence interval of parameter (Raue et al., 2009; Yang et al., 2016). We also benchmarked the PL approach with Fisher matrix approach.

The PL function is expressed as \[PL(\theta_i) = \min_{\theta_{-i}}log(likelihood(\theta_{-i})) \tag{1}\] where \(\theta \in (k_{deg}, k_{DR}, SC_{50}, S_{max})\), \(i \in (1,2,3,4)\), and \(\theta_{-i}\) denotes all parameters except the parameter under profiling (\(\theta_i\)).

We plotted the profile likelihood (PL) of all four parameters \(k_{deg}\), \(k_{DR}\), \(SC_{50}\), and \(S_{max}\), with further details listed in Table 1.

2.4.1 Parameter profile likelihood

With OneHealthStudio, we profiled likelihood of four systems parameters and checked confidence interval approximation from Fisher information matrix (FIM). As a validation check, we explored the simulations with parameters at the bounds of confidence intervals for each of two approaches. Figure 4 plotted the profile likelihood for four parameters.

Figure 4: Parameter uncertainty assessment using profile likelihood or fisher information matrix.

For \(k_{DR}\), both PL and FIM return identical confidence interval curves for cost functions, see the second left subplot in Figure 4. However, FIM of \(k_{deg}\) underestimates the parameter uncertainty as the FIM-based constructed parabola is above the profile likelihood curve, check zoom-in plot in Figure 5. Let us have a detailed check in Figure 5 For \(k_{deg}\), PL can construct more realistic parameter uncertainty bounds (shown in log10 scale in Figure 5), while FIM underestimate the 95% confidence interval.

Figure 5: PL approach realistically assesses uncertainty of \(k_{deg}\), while FIM underestimated 95% CI of \(k_{deg}\),

FIMs overestimates uncertainties of both \(SC_{50}\) and \(SC_{max}\), as the constructed paraloa is below the profile likelihood at the left side (towards 0, i.e. \(\log_{10}(SC_{50})\) and \(\log_{10}(S_{max})\) approaching to \(\infty\)), see \(SC_{50}\) in Figure 6(a) and case of \(SC_{max}\) in Figure 6(b).

Furthermore, we presented parameter uncertainty of \(SC_{50}\) and \(S_{max}\) in Figure 6,

# Read and display the first image
knitr::include_graphics("./png/siRNA_v20250315/benchmark_PL_FIM_2.png")

# Read and display the second image
knitr::include_graphics("./png/siRNA_v20250315/benchmark_PL_FIM_3.png")

(a) profiling \(SC_{50}\) in log10 scale

(b) profiling \(SC_{max}\) in log10 scale

Figure 6: PL approach realistically assesses uncertainty of \(SC_{50}\) and \(S_{max}\), while FIMs overestimate 95% CIs of \(SC_{50}\) and \(S_{max}\).

2.4.2 Simulation at profiled bounds with 95%

In profiling four parameters, \(k_{deg}\) and \(k_{DR}\) had well-bounded 95% CIs, see Figure 4. We showed the simulation at the profile bounds in Figure 7 and Figure 8, respectively. These simulation results confirmed the well-bounded confidence intervals of parameters \(k_{deg}\) and \(k_{DR}\).

Figure 7: simulation at the two-side profile bounds of \(k_{deg}\).

Figure 8: simulation at the two-side profile bounds of \(k_{DR}\).

Furthermore, we verified that FIM overestimated uncertainty of estimate of \(k_{deg}\), as simulations at PL bounds in Figure 7.

In two parameters \(SC_{50}\) and \(S_{max}\) with half-bounded PL shown Figure 4, would simulation at the bounds (@) also verified the cases.

Figure 9: simulation at the profile bound of \(SC_{50}\), which has only lower bound.

Figure 10: simulation at the profile bound of \(S_{max}\), which has only lower bound.

2.4.3 Model-based experimental design

From above section, PL-based uncertainty provided more realistic reflection of uncertainty around the parameters and showing. It was resulted from the fact that PL approaches takes nonlinearity of the biological mechanism into account. Section 2.4.1 with parameter uncertainty already show the capacity of PL approach to reveal realistic parameter uncertainty where FIM approaches failed.

In the context of experimental design, one can expect the impact of combination of readouts and dosing levels on profile likelihood curves. This could offer a way to better perform experimental design by taking into account with relatively small set of experimental data. Previous examples in systems biology applications have shown the applicability with profile likelihood to better balance the data and complexity of model structures (Litwin et al., 2022; Yang et al., 2016, 2021). With high-performance OneHealthStudio modeling and simulation capacity, we would experience PL-approach for other large-scale models.

3 Conclusion and Perspectives

Our study highlights OneHealthStudio with higher simulation efficiency and better estimation performance:

• 8-fold speedup in simulation runtime compared to OSPSuite-R,
• more accurate sensitivity available from automatic differentiation compared to finite difference approach,
• 5-fold fewer iterations in parameter estimation,
• more realistic assessment of parameter uncertainty in large scale model from profile likelihood-based approach compared to conventional Fisher information matrix.

In near future, we plan to extend use cases of OneHealthStudio to model viability in a mixed-effect modeling fashion. Given the scalability of automatic differentiation in OneHealthStudio, we will benchmark with Hamiltonian Monte Carlo Bayesian inference (Elmokadem et al., 2023) and nonlinear mixed-effect modeling (Fidler et al., 2019). One the technical side, OneHeathStduio currently only supports computation on CPU, while we are developing and scaling OneHealthStudio with jax framework (Bradbury et al., 2018), which would lead to utilize more power computational units like GPU and TPU. Together, we aim to deliver more advanced and reliable computational analysis workflows for mechanistic modeling applications like model informed drug development and next generation risk assessment approach (NGRA).

Appendices - PBPK of a microsomal triglyceride transfer protein (MTP) inhibitor (see details also in another note)

Recap of PBPK model of MTPi

We used a PBPK model for the small molecule lomitapide, a microsomal triglyceride transfer protein (MTP) inhibitor used in homozygous familial hypercholesterolemia (HoFH) patients. The PBPK model for lomitapide was constructed using Open Systems Pharmacology Suite (OSPSuite). Model building for rats used the small molecule workflow in PK-Sim, then model simulation and parameter estimation and visualization were conducted in OneHealthStudio. We summarized experimental protocol and data as the following

• four-days oral administration followed by another two-days PK measurements,
• six dosing levels: 0.625, 1.25, 2.5, 5.0, 10.0, and 20.0 [mg/kg/day],
• data extract from NDA.

fitting performance and sensitvity from AD and FD appraoches

We compared the predicted time-course PK for a range of dose levels, see Figure 11.

Figure 11: Time-course fitting of PBPK model of lomitapide to rat PK data.

We further compared the predicted PK with observed PK data in [\(\mu\)mol/L], see Figure 12.

Figure 12: comparison of PBPK model to rat PK data to repeated dosing data, dashed lines indicate 2 fold bounds.

sensitivity \(\frac{\partial[lomitapide]}{\partial\theta}\) w.r.t. parameter logP and fub

We examined the sensitivities w.r.t two parameters. OneHealthStudio provided accurate sensitivities using the AD approach, while the FD approach showed poor performance.

• lipophilicity (logP)

Figure 13: comparison of sensitivity from AD and FD approaches w.r.t lipophilicity logP

• Fraction unbound (fub)

Figure 14: comparison of sensitivity from AD and FD approaches w.r.t fraction unbound fub.

References

Bradbury, J., Frostig, R., Hawkins, P., Johnson, M. J., Leary, C., Maclaurin, D., Necula, G., Paszke, A., VanderPlas, J., Wanderman-Milne, S., & Zhang, Q. (2018). JAX: Composable transformations of Python+NumPy programs (Version 0.3.13) [Computer software]. http://github.com/jax-ml/jax

Branch, M. A., Coleman, T. F., & Li, Y. (1999). A subspace, interior, and conjugate gradient method for large-scale bound-constrained minimization problems. SIAM Journal on Scientific Computing, 21(1), 1–23.

Bücker, M. (2006). Automatic differentiation: Applications, theory, and implementations. Springer.

Eissing, T., Kuepfer, L., Becker, C., Block, M., Coboeken, K., Gaub, T., Goerlitz, L., Jaeger, J., Loosen, R., Ludewig, B., et al. (2011). A computational systems biology software platform for multiscale modeling and simulation: Integrating whole-body physiology, disease biology, and molecular reaction networks. Frontiers in Physiology, 2, 4.

Elmokadem, A., Zhang, Y., Knab, T., Jordie, E., & Gillespie, W. R. (2023). Bayesian pbpk modeling using r/stan/torsten and julia/sciml/turing. jl. CPT: Pharmacometrics & Systems Pharmacology, 12(3), 300–310.

Fidler, M., Wilkins, J. J., Hooijmaijers, R., Post, T. M., Schoemaker, R., Trame, M. N., Xiong, Y., & Wang, W. (2019). Nonlinear mixed-effects model development and simulation using nlmixr and related r open-source packages. CPT: Pharmacometrics & Systems Pharmacology, 8(9), 621–633.

Kapfer, E.-M., Stapor, P., & Hasenauer, J. (2019). Challenges in the calibration of large-scale ordinary differential equation models. IFAC-PapersOnLine, 52(26), 58–64.

Litwin, T., Timmer, J., & Kreutz, C. (2022). Optimal experimental design based on two-dimensional likelihood profiles. Frontiers in Molecular Biosciences, 9, 800856.

Loh, W.-L. (1996). On latin hypercube sampling. The Annals of Statistics, 24(5), 2058–2080.

Raue, A., Kreutz, C., Maiwald, T., Bachmann, J., Schilling, M., Klingmüller, U., & Timmer, J. (2009). Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood. Bioinformatics, 25(15), 1923–1929.

Raue, A., Schilling, M., Bachmann, J., Matteson, A., Schelke, M., Kaschek, D., Hug, S., Kreutz, C., Harms, B. D., Theis, F. J., et al. (2013). Lessons learned from quantitative dynamical modeling in systems biology. PloS One, 8(9), e74335.

Salim, E. L., Kristensen, K., & Sjögren, E. (2025). Whole-body physiologically based pharmacokinetic modeling of GalNAc-conjugated siRNAs. Pharmaceutics, 17(1), 69.

Serban, R., & Hindmarsh, A. C. (2005). CVODES: The sensitivity-enabled ODE solver in SUNDIALS. International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, 47438, 257–269.

Yang, H., Meijer, H. G., Buitenweg, J. R., & Van Gils, S. A. (2016). Estimation and identifiability of model parameters in human nociceptive processing using yes-no detection responses to electrocutaneous stimulation. Frontiers in Psychology, 7, 1884.

Yang, H., Van der Stel, W., Lee, R., Bauch, C., Bevan, S., Walker, P., Van de Water, B., Danen, E. H., & Beltman, J. B. (2021). Dynamic modeling of mitochondrial membrane potential upon exposure to mitochondrial inhibitors. Frontiers in Pharmacology, 12, 679407.