Experimentally calibrated multiscale model predicts schedule dependent drug combination effects.

Skip to main content An official website of the United States government Here's how you know Here's how you know Official websites use .gov A .gov website belongs to an official government organization in the United States. Secure .gov websites use HTTPS A lock ( Lock Locked padlock icon ) or https:// means you've safely connected to the .gov website. Share sensitive information only on official, secure websites. Search Log in Dashboard Publications Account settings Log out Search… Search NCBI Primary site navigation Search Logged in as: Dashboard Publications Account settings Log in Search PMC Full-Text Archive Search in PMC Journal List User Guide PERMALINK Copy As a library, NLM provides access to scientific literature. Inclusion in an NLM database does not imply endorsement of, or agreement with, the contents by NLM or the National Institutes of Health. Learn more: PMC Disclaimer | PMC Copyright Notice NPJ Syst Biol Appl . 2026 Mar 3;12:52. doi: 10.1038/s41540-026-00669-4 Search in PMC Search in PubMed View in NLM Catalog Add to search Experimentally calibrated multiscale model predicts schedule dependent drug combination effects Othmane Hayoun-Mya Othmane Hayoun-Mya 1 Barcelona Supercomputing Center (BSC-CNS), 1-3 Plaça Eusebi Güell, 08034, Barcelona, Spain 2 Universitat Politècnica de Catalunya - BarcelonaTech (UPC), Barcelona, Spain Find articles by Othmane Hayoun-Mya 1, 2 , Arnau Montagud Arnau Montagud 1 Barcelona Supercomputing Center (BSC-CNS), 1-3 Plaça Eusebi Güell, 08034, Barcelona, Spain 3 Institute for Integrative Systems Biology (I2SysBio), CSIC-UV, Valencia, Spain Find articles by Arnau Montagud 1, 3 , Alfonso Valencia Alfonso Valencia 1 Barcelona Supercomputing Center (BSC-CNS), 1-3 Plaça Eusebi Güell, 08034, Barcelona, Spain 4 Catalan Institute for Research and Advanced Studies (ICREA), 23 Passeig Lluís Companys, 08010, Barcelona, Spain Find articles by Alfonso Valencia 1, 4, # , Miguel Ponce-de-Leon Miguel Ponce-de-Leon 1 Barcelona Supercomputing Center (BSC-CNS), 1-3 Plaça Eusebi Güell, 08034, Barcelona, Spain Find articles by Miguel Ponce-de-Leon 1, ✉, # Author information Article notes Copyright and License information 1 Barcelona Supercomputing Center (BSC-CNS), 1-3 Plaça Eusebi Güell, 08034, Barcelona, Spain 2 Universitat Politècnica de Catalunya - BarcelonaTech (UPC), Barcelona, Spain 3 Institute for Integrative Systems Biology (I2SysBio), CSIC-UV, Valencia, Spain 4 Catalan Institute for Research and Advanced Studies (ICREA), 23 Passeig Lluís Companys, 08010, Barcelona, Spain ✉ Corresponding author. # Contributed equally. Received 2025 Sep 9; Accepted 2026 Feb 7; Collection date 2026. © The Author(s) 2026 Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/ . PMC Copyright notice PMCID: PMC13065843  PMID: 41776194 Abstract Therapeutic synergy emerges from interactions between molecular drug action, intracellular signaling, and tissue-level transport dynamics. We developed a multiscale model integrating these scales to predict schedule-dependent drug combination effect in the AGS cell line. Calibrated solely on single-drug growth curves, the model accurately predicted population-level outcomes of drug combinations without combination-specific training. This demonstrates the model’s capacity to suggest mechanistic multiscale insights into the logic of drug combinations within the AGS cell line, establishing a computational platform for the systematic in silico exploration of virtual multiscale experiments on drug diffusion and dosing schedules. Cross-scale analysis revealed that combination therapy efficacy flows across scales: population-level pharmacokinetics dictate the sequence of molecular target engagement within individual cells, determining collective cell-fate decisions. Our simulations predict that inhibiting the PI3K/AKT axis before MEK is more effective than the reverse order, disabling a pro-survival rebound and locking cells into an apoptotic state. Ultimately, by capturing phenomena that single-scale approaches cannot, this framework generates translationally relevant hypotheses, providing a versatile platform for optimizing drug scheduling, formulation, and combination strategies where matched molecular Boolean models and phenotypic data are available. Subject terms: Biophysics, Cancer, Computational biology and bioinformatics, Drug discovery, Systems biology Introduction A central challenge in oncology is predicting how cancer cells will respond to therapeutic interventions, particularly to combinatorial therapies designed to overcome drug resistance and improve patient outcomes 1 – 3 . While single-target targeted therapies can show initial efficacy, tumors frequently develop resistance through compensatory pathway activation and clonal selection 4 , leading to treatment failure 5 . Combination therapies have emerged as a promising strategy to address this challenge, targeting multiple pathways simultaneously to prevent or delay resistance 2 , 6 . However, therapeutic response is an inherently multiscale process, where drug-target interactions at the molecular level trigger complex signaling cascades that determine cell-fate decisions, ultimately shaping population-level dynamics and tissue-wide outcomes 7 . The synergistic effects of drug combinations often emerge from this intricate interplay across scales, making them difficult to predict from single-scale experiments alone 8 . Computational systems biology provides a powerful and increasingly essential toolkit for studying and understanding the effect and mechanistic underpinnings of synergistic drug combinations 9 . To capture the inherently spatial nature of tumor biology, agent-based models (ABMs) have become a widely-used approach 10 . In an ABM, individual cells are modeled as autonomous agents that interact with each other and their local microenvironment according to a set of prescribed rules. This bottom-up strategy allows for the simulation of emergent phenomena such as clonal selection, resource competition, and necrotic core formation that cannot be captured by traditional, non-spatial population models like ordinary differential equations (ODEs) 11 – 13 . To ground agent behaviors in established molecular biology, they can be coupled with intracellular network models 14 . Here, boolean models offer a particularly well-suited framework 15 , 16 . By representing signaling components as nodes that can be either ON or OFF, they capture the logical structure and causal relationships within signaling pathways that drive critical cell-fate decisions, such as proliferation and apoptosis 17 . The Boolean qualitative approach is invaluable when detailed, quantitative kinetic parameters are unknown or difficult to obtain, as is often the case 18 , 19 . Integrating qualitative and quantitative formalisms into hybrid multiscale models represents a frontier in computational systems biology of cancer, promising a mechanistic bridge from molecular perturbation to tissue-level response 10 , 20 , 21 . However, their transition from research tools into quantitatively predictive engines has been impeded by a set of fundamental challenges 9 . Firstly, the intertwined problems of model calibration and inter-model coupling. Then, multiscale models are often under-constrained, containing dozens to hundreds of parameters that span biological scales 22 . Moreover, estimating these parameters to reproduce experimental data is a major computational challenge requiring high-performance computing strategies to explore the parameter space effectively 23 – 26 . On top of this is the challenge of designing the cross-scale interface, which involves the set of transfer functions that translate discrete Boolean states into continuous biophysical parameters 21 , 27 . A poorly designed or calibrated interface between models can sever the mechanistic link between scales, leading to biologically unrealistic simulations and undermining the model’s predictive credibility 27 . Even if a model is rigorously calibrated and mechanistically robust, it faces a final, and perhaps most critical, translational gap: ensuring that its predictions can be meaningfully applied in real-world scenarios 28 . Translational potential is a particularly relevant issue for models of drug synergy. Nearly all synergy measurements are established in vitro, using well-mixed, two-dimensional monolayer cultures where cells experience uniform and constant drug concentrations 29 . This idealized scenario bears little resemblance to the reality of a solid tumor growing within the human body. In vivo, the complex 3D architecture of the tumor microenvironment, including a dense extracellular matrix and a poor, heterogeneous vasculature, creates significant transport barriers for effective drug delivery 30 . Consequently, drugs with different physicochemical properties will penetrate tumor tissue at different rates, creating complex and time-varying spatiotemporal gradients 31 . The optimal synergistic ratio identified in vitro may never be achieved simultaneously across the entire tumor 6 . This pharmacokinetic discrepancy is a primary suspect in the frequent clinical failure of combinatorial drug treatments that show immense promise in preclinical assays 32 , 33 . Addressing this specific and clinically relevant challenge is paramount for the field. In this work, we present a comprehensive multiscale modeling framework that systematically addresses this set of challenges. Our approach integrates a previously validated Boolean model of the AGS gastric cancer cell line 34 with a 3D agent-based model using the PhysiBoSS tool 21 , using the AGS Boolean network as a validated regulatory core for drug synergy prediction and applying the multiscale coupling to explore how these predictions translate to tissue-level behavior. To surmount the calibration and interface challenges, we implement a massively parallel pipeline using the EMEWS framework 25 , 26 to perform evolutionary optimization, systematically identifying the parameters that govern the crucial cross-scale interface. A key strength of our approach is that, although the resulting multiscale model is calibrated only on experimentally measured single-drug response time-courses from the AGS line, the model is validated by its ability to predict the temporal dynamics of drug combination outcomes at the cell population level that correspond to experimentally defined synergies. This enables us to predict the treatment response and analyze the temporal evolution of the drug interaction at the molecular, cellular, and population levels for different drug pairs, suggesting context-specific pathway dependencies in the two synergistic combinations tested (PI3K+MEK and AKT+MEK inhibitors). Building on this, we leverage the multiscale AGS model as an in silico laboratory to systematically probe the translational gap. This framework enables performing complex spatiotemporal simulations that single-scale Boolean approaches cannot execute, such as sweeping drug diffusion parameters, testing staggered dosing schedules, and quantifying asymmetric efficacy windows. Our simulations suggest that both the timing and sequence of synergistic drug administration substantially influence therapeutic outcomes in the AGS context. This approach underscores the necessity of accounting for pharmacokinetic variability and scheduling in the rational design of combination therapies. It further suggests that a calibrated multiscale framework may help generate quantitative, testable hypotheses that require experimental validation to inform more effective clinical strategies. Moreover, this PhysiBoSS and EMEWS multiscale model calibration and simulation methodology is potentially applicable to other cell lines or cancer models when supported by the experimental data needed to calibrate the model and a cell-line or cancer-type specific Boolean model. Results Development of a multiscale agent-based model for drug response in a gastric cancer cell line We developed a comprehensive, multiscale agent-based model that integrates molecular signaling networks, cellular behaviors, and population dynamics (Fig. 1 ). Using this model, we investigated the mechanisms of drug action and synergy in a gastric cancer cell line, as well as the effects of spatial and temporal heterogeneity on drug effectiveness. The model simulates a growing culture of AGS cells within a microenvironment where cells interact mechanically and where diffusible substances, such as nutrients and drugs, interact with individual cells (Fig. 1 A). Each cell is represented as an individual agent with defined physical properties (e.g., position, volume) and biological processes, including: i) proliferation and death (via apoptosis or necrosis), governed by an intracellular signaling network based on the AGS Boolean model (Fig. 1 A); ii) a pressure-dependent contact inhibition mechanism that modulates proliferation as cells approach confluence, essential for capturing realistic in vitro growth dynamics (Fig. 1 B); iii) a model of drug transport and binding kinetics (Fig. 1 C). The AGS Boolean model governs cell fate decisions by regulating proliferation and death signals in response to perturbations, including exposure to different kinase inhibitors. The network includes the known targets of the drugs used in this study: PI103 (PI3Ki), PD0325901 (MEKi), and AKT Inhibitor VIII (AKTi) (Supplementary Table 2 ). Fig. 1. Multi-scale hybrid model of AGS cancer cell line drug response experiments. Open in a new tab Our multiscale model integrates cellular, molecular, and microenvironmental scales to model drug responses and synergistic interactions. A Left: Simulation domain representing an approximation to an RTCA well-plate (610 × 610 × 20 μm 3 ) containing ~1500 cells in a disk of 305 μm radius with a 2.8 μm cell spacing. PI3K inhibitor (green) and MEK inhibitor (blue) were added to the microenvironment at 1280 simulation minutes and with the same I C 50 concentrations from the experimental reference. This setup allowed for replicating the in vitro setup in which the reference single and combined drug perturbations were performed. Right: Detailed cell agent model showing the embedded AGS Boolean network with interconnected signaling pathways, highlighting key drug target nodes (yellow) and model readouts categorized as antisurvival or apoptosis-promoting (red) and prosurvival or growth-promoting (green) signals that regulate apoptosis and proliferation, respectively. Transfer functions are represented with a sigmoidal curve in a box. B Chart representing the pressure-based growth inhibition model. C Drug transport and binding kinetics model showing diffusion across the cell membrane followed by reversible binding to target proteins with specific k 1 and k −1 rates, forming drug-target complexes [DT] that modulate Boolean network activity. We modeled drug transport within each agent assuming simple diffusion dynamics (Eq. ( 2 )) with molecules diffusing through the extracellular microenvironment and across cell membranes according to drug-specific permeability coefficients until reaching a steady-state with the extracellular medium (see Supplementary Material for a more detailed description). The internalized inhibitors bind to their molecular targets via a reversible kinetic model (Eq. ( 4 )), with the resulting drug-target complex concentration directly inhibiting specific target nodes in the intracellular Boolean network. To quantitatively link drug-target binding to signaling inhibition, we use a Hill equation transfer function that calculates the probability of target node inactivation based on drug-target complex concentration (Eq. ( 5 )). At each simulation time step, this probability is used in a stochastic process to determine whether the target Boolean node (PI3K, MEK, or AKT) becomes inactive, enabling coupling between drug presence and Boolean signaling responses. Further details on the binding model and the Hill equation implementation are provided in Methods section. The Boolean model regulates cell fate by generating multi-valued and not mutually exclusive readouts that can simultaneously denote proliferation and apoptosis. The proliferative signal is calculated from three key prosurvival Boolean nodes (cMYC, TCF, RSK), while the apoptotic signal is derived from three key antisurvival Boolean nodes (FOXO, Caspase8, Caspase9). To consolidate these complex signals into a single value for each process, we quantify proliferative and apoptotic signaling through weighted sums. These integrated scores, which reflect the varying contributions of each node, directly modulate cellular proliferation and apoptosis rates through two different Hill functions (see Methods for further details). Under control conditions, this model allows AGS cells to maintain their experimentally observed doubling times and basal apoptosis rates. When drugs are introduced, perturbations propagate through the Boolean network, altering the integrated scores and modulating cellular phenotypes via calibrated transfer functions. The calibrated model replicates control condition and single-drug experimental response patterns To establish a baseline for our simulations, we first calibrated the model to replicate control growth dynamics using the experimental growth curves where no drug was added. This calibration process involved optimizing five baseline growth parameters that govern contact inhibition and basal proliferation: the pressure half-max ( K p ) and Hill coefficient ( H p ), the basal proliferation rate ( r γ , b a s a l ), the initial cell population radius ( R p o p u l a t i o n ), and the intercellular spacing (Table 1 ). For this optimization, we utilized two distinct algorithms: the Covariance Matrix Adaptation Evolution Strategy (CMA-ES) and a Genetic Algorithm (GA) (see Methods section for more details), and both of them converged to a stable low RMSE value (Supplementary Figs. 1, 2 ). The resulting calibrated model successfully reproduces the experimental time course, achieving a Pearson correlation coefficient of 0.99 ( p < 0.05) and accurately capturing both the exponential growth phase and the gradual approach to carrying capacity, as observed in Fig. 3 . Table 1. Model parameters used in the multiscale PhysiBoSS simulation, grouped by functional role Parameter Description Lower bound Upper bound Units Baseline growth parameters K p Proliferation half-max 1.0 10.0 Dimensionless H p Proliferation Hill coefficient 1.0 10.0 Dimensionless r b a s a l , γ Basal proliferation rate 1 × 10 −4 1 × 10 −3 1/min R p o p u l a t i o n Initial population radius 100.0 500.0 μm Cell spacing Inter-cellular distance 1.0 5.0 μm Drug transport and binding parameters k d r u g , X Drug permeability 1 × 10 −7 1 × 10 −4 μm/min K 1, d r u g , X Drug binding rate 1 × 10 −9 1 × 10 −4 M −1 ⋅ min −1 K −1, d r u g , X Drug unbinding rate 1 × 10 −9 1 × 10 −4 min −1 Cellular response parameters H γ Hill coefficient for growth 0.1 15.0 Dimensionless K A . γ Half activation for growth 1 × 10 −5 1.0 Dimensionless r r e s p o n s e , γ Growth response rate 0.01 100.0 1/min H α Hill coefficient for apoptosis 0.1 15.0 Dimensionless K A , α Half activation for apoptosis 1 × 10 −5 1.0 Dimensionless r m a x , α Maximum apoptosis rate 1 × 10 −3 1 × 10 −2 1/min r r e s p o n s e , α Apoptosis response rate 0.01 100.0 1/min Boolean network integration parameters ω γ , c M Y C cMYC influence on growth 0.0 1.0 Dimensionless ω γ , T C F TCF influence on growth 0.0 1.0 Dimensionless ω γ , R S K RSK influence on growth 0.0 1.0 Dimensionless ω α , F O X O FOXO influence on apoptosis 0.0 1.0 Dimensionless ω α , C a s p 8 Caspase8 influence on apoptosis 0.0 1.0 Dimensionless ω α , C a s p 9 Caspase9 influence on apoptosis 0.0 1.0 Dimensionless Boolean model simulation parameters d t i n t r a c e l l u l a r Intracellular timestep 1.0 40.0 Min s c a l i n g i n t r a c e l l u l a r Boolean to continuous scaling 1.0 1000.0 Dimensionless Open in a new tab The first group governs drug transport and target binding, the second group controls cellular growth and apoptotic responses, and the third group mediates the connection between Boolean network states and continuous cellular phenotypes. Specific parameter values are based on literature sources for variables regarding basal proliferation rates 35 – 38 ,drug permeability 39 – 44 ,drug binding and unbinding rates 45 , 46 and apoptosis rates 47 – 49 . Fig. 3. Population-level dynamics in drug combinations recover experimentally-quantified synergy in a model calibrated using single drug experiments. Open in a new tab The multiscale model, using parameters derived exclusively from single-drug calibrations, successfully recapitulates the population-level dynamics effects of combination therapies. Dotted lines indicate experimental curves. Solid lines indicate simulated curves. A For the PI3Ki-MEKi combination, the simulation (solid orange line) accurately reproduces the experimental data (dotted orange line), showing a strong growth arrest leading to a plateau significantly below the levels of either single agent (blue and green) or the control (gray). B Similarly, for the AKTi-MEKi combination (solid orange line), the model predicts the profound growth inhibition observed experimentally, with the cell count stabilizing at a level markedly lower than that of the single-agent or control conditions. Shaded regions represent variability across experimental replicates and aggregated top simulated growth curves. Building on this calibrated baseline, we then calibrated the parameters for individual drug responses. Using the same CMA-ES and GA optimization algorithms, we explored an 18-parameter space for each single-drug perturbation (PI3Ki, MEKi, and AKTi) to fit experimental dose-response data. These parameters govern drug-specific internalization, target-binding, apoptosis, and proliferation modulation (Table 1 ). While both optimization algorithms started with similarly high RMSE values, CMA-ES consistently achieved better final RMSE values across all single-drug conditions. The convergence analysis confirmed that CMA-ES maintained a steeper descent toward optimal solutions, an advantage that was particularly evident in the AKT inhibition scenarios (See Supplementary Fig. 1 ). For each single-drug condition, we selected the top 10% best parameter sets, identified by the lowest root-mean-square error (RMSE) values from our optimization runs. The growth curves generated by these best-performing parameter sets are shown in Fig. 2 . With these parameters, our simulations quantitatively reproduce the experimental results for the three conditions analyzed. Specifically, the Pearson correlation coefficients between the average simulated growth curves and the experimental data are 0.984, 0.977, and 0.991 for PI3Ki, MEKi, and AKTi, respectively (all p < 0.001). Furthermore, a key finding from our independent calibrations is the significant overlap observed in the distributions of these top-performing parameter sets across the three different inhibitors. This suggests that a core set of parameters remains consistent regardless of the specific drug perturbation. A detailed statistical analysis and visualization of these parameter distributions can be found in the Supplementary Material (Supplementary Fig. 6 , Supplementary Table 6 ). Fig. 2. Single drug calibration and emergent synergistic effects in drug combinations. Open in a new tab Comparison of simulated growth curves (from the top 10% of calibrated parameter sets) with experimental data for three single-agent kinase inhibitors. A For PI3K inhibition (blue), the model (line with points) reproduces the cytostatic effect, capturing the final growth plateau of the experimental data (solid line) (Pearson’s r = 0.984). B For MEK inhibition (green), the model captures the complex, non-monotonic response where cell count peaks before declining to a stable level (Pearson’s r = 0.977). C For AKT inhibition (purple), the model accurately replicates the sustained growth suppression that does not result in a plateau (Pearson’s r = 0.991). In all panels, experimental (dashed) and simulated (solid) control growth curves are shown in gray. The vertical line indicates the time of drug administration. Shaded regions represent variability across experimental replicates and aggregated top simulated growth curves. To evaluate the calibrated model under perturbation, we ran simulations for each single-drug condition by sampling parameters from an empirical distribution constructed using the top 10% of best-fitting sets. The resulting time-course predictions are shown in Fig. 2 . For PI3K inhibition (Fig. 2 A), the experimental curve shows a sharp reduction in growth rate leading to a plateau. Our simulation quantitatively reproduces the final cell count at this plateau but does not fully capture the rapidness of the growth stabilization seen experimentally toward the end of the time course. In the case of MEK inhibition (Fig. 2 B), the experiment shows a complex, non-monotonic pattern where cell count peaks and then declines to a stable level. The model successfully captures this overall trend, including the final plateau, although there is a slight discrepancy in the timing of the growth deceleration observed between 2500 and 3000 min. Finally, for AKT inhibition (Fig. 2 C), the data demonstrates a consistent reduction in growth rate without reaching a definitive plateau. The model accurately replicates this sustained but slower growth dynamic throughout the simulation. Emergence of drug synergy from a multiscale mechanistic model After calibrating the parameters for each single-drug treatment by fitting the experimental curves, we focused on investigating whether the model could reproduce experimentally observed synergistic drug combination effects at the population level without being explicitly calibrated on drug combination data. Our model utilizes 15 shared parameters describing core cellular pathways and 3 drug-specific parameters for each inhibitor, which govern drug transport and binding (Supplementary Table 4 ). To test whether we could capture the synergy for a specific drug combination, such as PI3Ki-MEKi, we defined a dedicated parameter space for that pair. First, for the parameters shared between both drugs, we derived their distributions by comparing the top 10% of parameter sets from the single-drug calibration experiments for each individual inhibitor. In other words, the consensus distribution for each common parameter was defined as the union of its best-fitting values from both single-agent calibrations (Supplementary Table 5 ). We then supplemented these with the top 10% of each drug-specific parameter for each inhibitor in the combination. This approach, applied to each drug pair, results in a distinct parameter space for every synergy, built exclusively from single-agent data. In this way, synergy is tested as an emergent property of the model, rather than being imposed through direct calibration on combination data. A more detailed explanation of how these parameter spaces are defined is provided in the Methods section and Supplementary Material (See section “Model calibration workflow for synergy emergence”). We then used these dedicated parameter distributions to test the model’s accuracy. For each drug combination, such as PI3Ki-MEKi, we sampled 5000 parameter sets from the empirical distribution and ran simulations for the constituent single agents as well as the combination. The simulations using these parameter sets accurately reproduced the corresponding single-drug experimental curves, confirming that the consolidated parameter space preserved single-agent dynamics. The combination simulations, in turn, quantitatively reproduced the enhanced growth inhibition characteristic of the synergy, achieving strong correlations with the experimental data for both the PI3Ki-MEKi and AKTi-MEKi combinations (Pearson’s r = 0.998 and 0.948, respectively; both p < 0.001). This result is visually represented by the marked reduction in cell proliferation relative to both single-agent and control conditions (Fig. 3 , orange curves) and shows that synergistic effects emerge from the model without explicit calibrating on combination data. Furthermore, an analysis of the best-performing parameter sets from both synergy simulations revealed that their distributions were remarkably similar (see Supplementary Material , Supplementary Table 7 ), despite the fact that both combinatorial therapies use different mechanisms as we will show in following sections. The simulations qualitatively reproduced the temporal dynamics observed in the experimental synergy curves for both drug pairs (Fig. 3 ). For the PI3Ki-MEKi combination (Fig. 3 , left panel), the model successfully captured the strong synergistic interaction. The simulation correctly reproduced the initial slow growth followed by a sustained growth arrest, leading to a plateau at a normalized cell count of ~50. This outcome is substantially lower than the cell counts achieved with either PI3Ki or MEKi alone. A similarly potent synergy was reproduced for the AKTi-MEKi combination (Fig. 3 , right panel). In this case, the simulation also showed a profound growth inhibition, with the cell count stabilizing at a normalized value below 40 cells, a level markedly lower than that observed for either AKTi or MEKi as single agents. Multiscale model proposes drug-specific signaling dynamics and cell fate decisions To dissect the mechanisms of drug action, we employed our calibrated multiscale model to systematically analyze how different treatments impact intracellular signaling and cell-fate decisions. Therapeutic efficacy is largely determined by a drug’s ability to modulate cell proliferation and apoptosis. Our model provides insight into these processes by quantifying the dynamic activity of key prosurvival (cMYC, TCF, RSK) and antisurvival (FOXO, Caspase8, Caspase9) Boolean model nodes, their integrated effects on cellular growth and apoptosis rates (summarized by the metrics S p r o and S a n t i ), and the resultant changes in viable cell populations, obtaining a cross-scale overview of the single and synergistic drug perturbations. The temporal dynamics of key signaling and phenotypic variables for each treatment are presented in Fig. 4 . PI3K inhibition, for example, induced a potent cytostatic effect, characterized by a sharp, sustained reduction in the cellular growth rate and a marginal increase in apoptosis (Fig. 4 , first column). This phenotype was underpinned by a strong suppression of pro-survival signaling ( S p r o decreased to 0.15) with only a minor corresponding increase in anti-survival signaling ( S a n t i = 0.10). Consequently, PI3K inhibition alone was sufficient to arrest cell population growth but not to induce notable regression. In contrast, AKT inhibition produced a different dynamic profile (Fig. 4 , second column), causing a rapid but transient suppression of pro-survival signals ( S p r o < 0.3) and a spike in anti-survival signals ( S a n t i around 0.40). However, this effect was not sustained, as pro-survival pathways recovered, allowing for the resumption of cell population growth after the initial disruption. MEK inhibition elicited yet a different, primarily pro-apoptotic, response (Fig. 4 , third column). While the suppression of growth-related signals was more gradual than in the PI3Ki cases ( S p r o decreased to around 0.4), it was accompanied by a robust induction of anti-survival signaling ( S a n t i arouns 0̃.25). This signaling profile translated to a net increase in the apoptosis rate, resulting in a moderate but steady reduction in the viable cell population. Fig. 4. Cross-scale modeling of drug effects: from signaling to population dynamics. Open in a new tab The multiscale model captures the effects of three single kinase inhibitors (PI3Ki, MEKi, AKTi) and two drug combinations (PI3Ki-MEKi, AKTi-MEKi), shown in columns. The leftmost schematic illustrates how Boolean readout nodes are combined to determine cell fate decisions. For each treatment, the top row shows the temporal evolution of pro-survival (green) and anti-survival (orange) Boolean readout node signals across the cell population. The middle row presents the corresponding phenotypic rates: proliferation (green) and apoptosis (orange), derived from these signals via transfer functions. The bottom row displays the resulting population dynamics, with live cell count (green) and apoptotic cell count (orange) over time. Drug administration occurs at t = 1200 min (vertical line). Single-drug treatments reveal distinct signaling and phenotypic profiles, while drug combinations produce enhanced suppression of growth and increased apoptosis, leading to strong population-level inhibition. Shaded bands indicate standard deviation from simulation replicates. Analysis of the combination therapies suggested a mechanistic basis for their observed synergetic effect (Fig. 4 , fourth and fifth columns). Effective combinations successfully merged the distinct advantages of each single agent to produce a superior therapeutic outcome. The PI3Ki-MEKi combination, for example, paired the strong, sustained suppression of pro-survival signaling seen with PI3Ki with the induction of apoptosis characteristic of MEKi. This dual-pronged attack on both growth and survival pathways resulted in a significantly elevated apoptosis rate and substantial cell culture regression, an effect not achievable by either drug alone. Likewise, the AKTi-MEKi combination achieved a powerful cytotoxic effect by robustly suppressing pro-survival signals while simultaneously elevating anti-survival signaling, demonstrating a clear synergistic advantage over monotherapy. To generate hypotheses about the network-level basis for these synergistic responses, we analyzed the calibrated weights of the individual signaling nodes. This analysis suggests that the synergy is achieved not just by amplifying existing signals but by shifting the relative importance of different pathways. For example, to reproduce the PI3Ki-MEKi synergy, the model required a reduced contribution from the pro-survival node cMYC and an increased influence from the Caspase9 apoptotic pathway. Similarly, the AKTi-MEKi synergy was best represented by a state where the influence of the Caspase8 anti-survival node became more prominent (see Supplementary Material , Supplementary Figs. S9 , S10 for the node weight distributions). These findings provide specific, testable hypotheses about how signaling networks may adapt to produce a synergistic response. Sensitivity analysis reveals the mechanistic drivers of therapeutic response Having calibrated our multiscale model to capture a range of drug-specific responses, we next sought to identify the most influential parameters that govern its behavior. To achieve this, we performed a machine learning-based sensitivity analysis (using a Random Forest model with SHAP) to quantify how each component of our model influences the ability to reproduce the experimental data (Supplementary Fig. 8 ). Based on our sensitivity analysis, we hypothesized that the two synergistic combinations operate through distinct biological mechanisms. While the general regulation of apoptosis was paramount for all treatment responses (Supplementary Fig. S8A ), the specific drivers of synergy differed for each drug pair. The PI3Ki-MEKi synergy was found to globally prime cells for death; it simultaneously lowered the threshold for committing to apoptosis (Half-max for apoptotic response, K A , α ) and increased the rate of its execution ( r α , max ) (Supplementary Fig. 8B ). The AKTi-MEKi combination (Supplementary Fig. S8C ), however, achieved its effect by amplifying the FOXO-mediated cell death axis, indicating that synergy is achieved by making the network highly reliant on this specific cancer suppressor pathway. The model’s ability to resolve such distinct mechanisms was further validated by its capacity to explain drug-specific phenotypes, quantitatively accounting for the distinct transient signaling profile observed under AKTi treatment through the unique influence of the TCF pro-growth pathway. See Supplementary Material for an in-depth analysis of the SHAP distributions. Differences in scheduling time affect synergistic efficacy To investigate how treatment scheduling and biophysical properties of each drug influence therapeutic synergy, we developed a three-dimensional (3D) simulation of a solid cancer environment. Building on our 2D model calibrated with in vitro data 34 , this 3D setup allowed us to systematically test the effects of the drug addition timing and diffusion coefficients on PI3Ki-MEKi and AKTi-MEKi efficacy (Supplementary Fig. 13 ). A complete description of the model parameters can be found in the “Model Parameter Calibration” subsection in the Methods section. First, we confirmed that the efficacy of the simulated drug combination still recovers the experimentally defined synergistic drug effect in our 3D simulation setting (Supplementary Fig. 13 ). Applying the inhibitors in combination reduced the final live-cell fraction to below 10%. This was significantly more effective than either single-agent treatment ( p < 0.005), although all treatments markedly outperformed the no-drug control. Next, we investigated how treatment efficacy depends on the administration timing of the synergistic drugs. For a fixed diffusion coefficient (600 μ 2 m/min), we compared three distinct drug administration schedules: (1) PI3Ki/AKTi-first, in which PI3K or AKT inhibitor is administered prior to MEK inhibitor; (2) MEKi-first, where the MEK inhibitor is given before the PI3K or AKT inhibitor; and (3) simultaneous administration, where both inhibitors are introduced at the same time (Fig. 5 B). The PI3Ki/AKTi-first and simultaneous schedules were predicted to be highly effective, both resulting in a final cell survival of less than 15%. In opposition, reversing the order to a MEKi-first schedule increased cell survival to ~33% ( p < 0.001). A more detailed simulation sweep of drug timing confirmed this dependency, predicting that efficacy begins to decline when MEKi precedes the second drug by approximately six hours, whereas the PI3Ki/AKTi-first schedule was consistently predicted as effective as, or marginally better than, simultaneous administration (Supplementary Figs. 14, 15 ). Fig. 5. Timing of drug administration critically determines therapeutic synergy efficacy in a 3D AGS cell population model. Open in a new tab A Schematic of the 3D simulation environment. Drugs are administered at the top surface and diffuse into a layer of AGS cells, whose individual states are tracked throughout the simulation. B Final live cell count (% of no-drug control) for AKTi+MEKi (left) and PI3Ki+MEKi (right) combinations under three administration schedules. The results demonstrate that the MEKi-first schedule is significantly less effective than administering the PI3K/AKT inhibitor first or simultaneously. Bars represent mean ± s.d. ( p < 0.001, p < 0.01). C Multi-scale time-course analysis for the six scenarios presented in ( B ). Each column corresponds to a specific drug combination and timing schedule. Rows display dynamics at different biological scales, including pro- and anti-survival node activation, target node activity, intracellular signal levels, cellular rates (growth, apoptosis), and population counts (alive, apoptotic). The MEKi-first schedules (third and sixth columns) uniquely show a rebound in pro-survival signaling and target node activity. This leads to attenuated apoptosis and a higher final cell count, providing a direct mechanistic explanation for the reduced efficacy observed in ( B ). To understand the mechanistic basis for this timing-dependent synergy, we analyzed the dynamics of key signaling nodes from the pathway and cell-fate network (Fig. 5 C). In the highly effective PI3Ki/AKTi-first and simultaneous schedules, the treatment predicted a rapid and sustained inhibition of the target nodes’ activity. This initial perturbation cascaded through the network, causing a shutdown of growth-promoting signals (c-MYC, RSK and TCF activation levels) and, concurrently, a lasting activation of apoptosis-increasing signals (FOXO, Caspase8 and Caspase9 activation levels). This sustained pro-apoptotic drive agrees with the high treatment efficacy, with cell survival plateauing below 15%. In contrast, the MEKi-first schedule failed to fully inhibit prosurvival pathways. Initial MEK inhibition allowed PI3K/AKT activity to persist at ~50%. This partial suppression was insufficient to fully inhibit downstream pro-survival signals, which only dipped transiently, and it resulted in a delayed and attenuated activation of apoptosis-promoting cascades. The addition of PI3Ki/AKTi predicted a secondary, but much weaker and delayed apoptotic signal. This incomplete response ultimately permitted cell recovery and is consistent with the significantly higher final survival in this scenario, our multiscale model predicts. These results indicated that the timing effect emerges from the way the signaling network integrates drug-induced perturbations through the interconnected PI3K/AKT and MAPK pathways and the regulatory feedback loops between them. To investigate this, we performed Boolean network simulations by perturbing the target node at different time points in a similar way as the in silico experiment. These simulations allowed us to identify a critical negative feedback loop centered on the IRS1 protein (Supplementary Fig. 16 ), an activator of the PI3K/AKT pathway that is suppressed downstream of MEK. These results explained the scheduling effect: administering MEKi first removes this suppressive feedback, permitting a transient pro-survival rebound in PI3K/AKT signaling that compromises overall efficacy. Conversely, targeting the PI3K/AKT pathway first or simultaneously is more effective because it preemptively disables this escape route. The relative mobility of the drugs further shaped the optimal schedule. By simulating a wide range of diffusion coefficients for both inhibitors, we identified a general principle for effective synergy in our model: PI3K/AKT signaling must be suppressed throughout the cell population before ERK activity (targeted by MEKi) significantly decreases (See Supplementary Figs. S14 , S15 ). Consequently, if the PI3K/AKT inhibitor diffused as fast or faster than the MEK inhibitor, simultaneous administration was effective. However, when the MEK inhibitor was the faster agent, the PI3K/AKT inhibitor had to be administered first to compensate for its slower penetration and preempt a pro-survival rebound. This principle proved robust; even in a simulated “worst-case" scenario with a large diffusion mismatch, a staggered schedule that ensured the correct drug arrival order at the cellular level successfully restored therapeutic efficacy. Discussion In this work, we present a model development protocol implemented within the PhysiBoSS framework to investigate drug synergy in a gastric cancer cell line. Our protocol integrates a Boolean network model of key intracellular signaling pathways with an agent-based model of cell behaviors. A central methodological advance of this study is the improved coupling between molecular and cellular scales, using a previously validated gastric adenocarcinoma Boolean network as the foundation and focusing on how faithfully the multiscale extension translates known intracellular predictions into tissue-scale outcomes. Specifically, we employ Hill functions to translate extracellular drug concentrations into the inhibition for intracellular targets such as PI3K, MEK, and AKT. The resulting changes in network state are then mapped to cell fate decisions using a weighted, multi-node readout: the collective activity of pro-apoptotic nodes (FOXO, Caspase8, Caspase9) determines the apoptosis rate, while pro-survival nodes (TCF, RSK, cMYC) determine the proliferation rate 12 , 50 . By calibrating the weights assigned to these signaling outputs, this model development protocol enables a more nuanced and experimentally faithful mapping between molecular states and cellular phenotypes, compared to single-node readouts 20 , 51 . Building on this methodological advance, we employed our multiscale model to accurately predict the temporal dynamics of drug synergies at the cell population level. Crucially, the model was calibrated solely on single-drug data, making these population-level combination predictions truly predictive rather than fitted to combination effects, with the resulting outcomes matching the synergistic classifications established by Loewe analysis in ref. 34 . The model reproduced the experimental synergy curves with high fidelity and also provided cross-scale insight into the underlying mechanisms. In the following sections, we discuss our contributions to model calibration, drug synergy analysis, and temporal dynamics of drug administration in the AGS system. These results demonstrate the effectiveness of our model development protocol, combining multiscale calibration with the ability to predict population-level combination outcomes that correspond to experimentally defined synergies. Our multiscale model, calibrated against experimental data for single kinase inhibitors (MEKi, PI3Ki, AKTi) selected from the 34 panel based on the availability of RTCA time-series data and confirmed Boolean-level phenotypic responses, successfully reproduced their individual effects on cell population growth. A key finding is that the calibrated model also recovers the synergistic growth inhibition observed in PI3Ki-MEKi and AKTi-MEKi combination therapies, an emergent behavior that was not explicitly part of the training data (Figs. 2 and 3 ). This finding suggests that the model can hint towards fundamental mechanisms of drug action rather than simply fitting to observed data. Further analysis revealed a common, overlapping parameter space that robustly described both single-agent and combination effects (Supplementary Tables S6 , S7 ). The existence of this shared set of parameters suggests that the observed synergies arise from the same underlying biological network that governs individual drug responses, reinforcing the model’s biological plausibility 52 , 53 . The initial calibration against single-drug experimental data constrains the parameter space within empirically-informed bounds, establishing a foundation for predictive modeling. We systematically evaluate parameter robustness by validating parameter sets based on their ability to generate simulations that closely match observed single-drug responses. The predictive validation process evaluates whether models calibrated on single-drug data can predict population-level temporal dynamics of drug combination effects. While synergies are predicted by the Boolean model, the resulting population-level growth and apoptosis dynamics emerge from a parameter space calibrated solely on single-drug responses. The model’s predictive capability validates that it can point towards a more nuanced understanding of time-resolved molecular and single-cell level changes across temporal and spatial scales in the AGS model. Furthermore, the model can provide mechanistic hypotheses regarding distinct temporal signatures of growth inhibition for each drug. PI3K and MEK inhibitors induced a response plateau, characteristic of a predominantly cytostatic effect consistent with their known roles in regulating cell cycle progression 52 , 54 . In contrast, AKT inhibition produced a more linear and sustained decline in the cell population, suggesting a stronger induction of apoptosis. This distinction reflects the differential engagement of pro-survival and pro-apoptotic pathways by each compound, governed by AKT’s central regulatory role 54 , 55 . These emergent dynamics are linked to the model’s architecture, which separates cell fate decisions from their subsequent execution. The Boolean network module evaluates intracellular signals to determine a cell’s fate (e.g., commitment to irreversible apoptosis), while the agent-based module simulates the temporal progression of that fate. This structure allows the model to deconvolve population-level outcomes into the underlying competition between cytostatic and apoptotic responses. The model thus illustrates how combining the cytostatic pressure of MEKi with the pro-apoptotic pressure of PI3Ki or AKTi effectively overcomes compensatory survival signaling, leading to the potent and sustained synergistic growth suppression observed in dual-pathway inhibition strategies 52 , 53 , 56 . The choice of calibration algorithm was critical for uncovering these mechanisms. Parameter estimation for complex biological models is a well-established challenge where the performance of different optimization methods can vary significantly 57 . While both GA and CMA-ES were evaluated, CMA-ES proved superior for this complex, high-dimensional parameter space (Supplementary Fig. S1 ). Although GA showed faster initial convergence, CMA-ES ultimately identified solutions with significantly lower error by sustaining a broader exploration of the parameter landscape. This result is consistent with other findings that highlight the robustness of CMA-ES in demanding biological modeling contexts 20 , 58 – 60 . The implications of this performance difference are significant for the model’s interpretation. The wider parameter distributions identified by CMA-ES represent multiple, mechanistically distinct parameter sets that can all yield high-quality fits to the experimental data. This is a known feature of complex biological models often described as “sloppy" 22 , 61 . In contrast, the tendency of GA to rapidly exploit narrow regions of the parameter space risks premature convergence on solutions that, while predictive, may not represent the full range of plausible underlying biology 60 . Therefore, our findings underscore that the exploration-exploitation trade-off, managed by the optimization algorithm, is a critical factor that determines not only a model’s predictive accuracy but also its capacity for mechanistic discovery. Future work could benefit from exploring hybrid strategies that leverage the initial exploratory power of GA with the local refinement capabilities of CMA-ES 62 . While the potential for drug synergy may be encoded within the static topology of a signaling network 34 , the Boolean layer alone has several limitations. It operates in abstract update steps without biological time scales, uses binary perturbations that cannot capture continuous dose responses, and assumes well-mixed conditions, ignoring spatial heterogeneity in drug distribution and cell-cell mechanical interactions. Our multiscale analysis reveals that the ultimate realization of synergy is governed by temporal dynamics and spatial context, requiring a framework that bridges intracellular signaling to population-level outcomes. We show that synergy arises from the intracellular signaling response to the population-level cell-fate decisions. This cross-scale link is mediated by a calibrated interface of model parameters 50 , 51 . Specifically, our model demonstrates how effective drug combinations can reshape the cell’s temporal signaling profile, forcing a transition from transient, recoverable states into stable, apoptotic attractors, which has been experimentally shown in various cancer types 63 – 65 . By mechanistically linking a specific drug-induced signaling signature to its macroscopic outcome, our framework provides a powerful lens through which to understand and predict the population-level outcomes. Zooming in on the molecular details, our model suggests that this commitment to apoptosis is achieved through distinct, drug-specific mechanisms. The analysis indicates that different synergistic combinations modulate the network in unique ways to enforce a decisive antagonism between growth-promoting and pro-apoptotic pathways. Our model reveals that the two synergistic combinations achieve their efficacy through distinct biological routes. The PI3Ki-MEKi synergy appears to emerge from a global sensitization of the cell’s apoptotic machinery. In this case, the model’s output was primarily controlled by parameters governing the general apoptotic capacity, which our analysis linked to an increased reliance on the intrinsic, Caspase-9-mediated pathway (Fig. 4 and Supplementary Fig. S10 ). This finding is consistent with experimental studies showing that dual PI3K-MEK inhibition powerfully induces apoptosis 66 . This suggests a mechanism based on priming the cell to lower its overall threshold for apoptosis, which is then executed via Caspase-9. In contrast, the AKTi-MEKi combination appears to exploit a more specific signaling vulnerability. Our model suggests this synergy is achieved by suppressing the oncogenic driver cMYC while activating pro-apoptotic FOXO signaling, a dynamic consistent with the known mutually inhibitory relationship between these pathways 67 . Critically, the model indicates this signaling shift creates a strong dependency on the extrinsic Caspase-8 pathway to induce cell death (Supplementary Fig. 11 ). This finding aligns with the general principle that AKT inhibition can sensitize cancer cells to therapy-induced apoptosis 68 . This points to a different therapeutic logic: the dynamic creation of a synthetic lethal state where the cell’s fate is critically dependent on the Caspase-8 checkpoint. The ability of the calibrated model to uncover these subtle, yet critical, distinctions in the mechanism underscores its value as a tool for generating precise, context-specific hypotheses about the nature of drug synergy. A significant insight from our work is that therapeutic synergy is governed not only by the magnitude of intracellular signals but also by the dynamic, cross-scale communication that dictates the final phenotypic response. Uncovering this required a key methodological contribution: integrating our mechanistic model with a machine learning-based sensitivity analysis using Random Forest and SHAP. This approach revealed that parameters controlling apoptosis were the most critical drivers of the model’s behavior. This finding underscores a deeper principle in cancer therapy, which is that a cell’s intrinsic readiness to undergo apoptosis, represented here by the high sensitivity of these interface parameters, can be a more critical determinant of drug efficacy than the initial signal suppression itself 69 . Still, these findings should be interpreted within the appropriate context. The parameter hierarchies identified by our sensitivity analysis reflect the mechanisms most critical for reproducing the specific experimental conditions used for calibration, rather than providing a global characterization of the system. This approach, in line with mechanistic learning frameworks 70 , help clarify how parameters governing apoptosis and the mapping from network states to phenotypes shape population-level outcomes in our AGS model. Using the calibrated multiscale framework as an in silico laboratory, we model drug diffusion, transport, and binding dynamics to connect Boolean network readouts to phenotypic traits. This enables us to investigate spatiotemporal constraints, such as spatially resolved drug delivery, contact inhibition, and dose scheduling, which the Boolean network alone cannot address. Using our calibrated model in a 3D setting, we systematically performed a series of in silico experiments to investigate the impact of drug scheduling on therapeutic efficacy. Unlike the reference experiments, which were limited to simultaneous addition of synergistic drug combinations, our simulations enabled us to test a wide range of administration sequences, diffusion coefficients and timing intervals for PI3K/AKT and MEK inhibitors. These simulations predicted that the sequence of pathway inhibitions can determine the therapeutic outcome, with the greatest cell death observed when the PI3K/AKT inhibitor was given before or together with the MEK inhibitor. A central objective of this work is to demonstrate how such a multiscale framework can generate quantitative, mechanistic hypotheses that guide experimental prioritization. While these computational predictions provide actionable starting points for laboratory investigation, we acknowledge that direct experimental validation remains essential to confirm the predicted timing dependencies, quantify the actual efficacy windows, and assess whether the observed asymmetries translate to in vitro and ultimately in vivo settings. Our study suggests a mechanistic rationale for the critical role of drug scheduling in synergistic cancer therapy. The model pointed towards a negative feedback loop, where ERK suppresses the PI3K activator IRS1, as a potential control point for this schedule-dependence 55 . Our cross-scale analysis of the simulation results provides a computational hypothesis regarding how this single feedback event dictates the downstream cell-fate decision (Fig. 4 ). When the MEK inhibitor is administered first, the resulting PI3K/AKT rebound is sufficient to prevent the complete shutdown of pro-survival nodes (e.g., c-MYC, RSK), while simultaneously delaying the activation of key apoptotic effectors like FOXO and Caspase8. This incomplete response permits cell recovery. Conversely, inhibiting the PI3K/AKT axis first preemptively disables this escape route, triggering a rapid and sustained inhibition of pro-survival signaling and a robust activation of the apoptotic machinery, committing cells into a state of high cytotoxicity. This computational finding offers a testable hypothesis for the observed timing-dependent synergy that warrants experimental investigation. Furthermore, our multiscale framework illustrates how this intracellular rule is modulated by the physical realities of drug transport in a 3D environment. The full multiscale simulations revealed an asymmetric temporal window of efficacy and enabled us to quantify the pharmacodynamic properties of the synergistic drug pair. By systematically sweeping timing offsets in a battery of PhysiBoSS experiments, we found that efficacy begins to decline when MEKi precedes the second drug by approximately six hours, whereas the PI3Ki/AKTi-first schedule remains consistently effective. We then explored the underlying molecular mechanism with simplified synchronic Boolean simulations, which pointed towards the role of the ERK-IRS1 feedback loop, but operate in an abstracted system without biological time or spatial context. This quantitative drug profiling, which provides biologically coherent and potentially translational insights into scheduling delays and differential drug diffusion, benefits from the multiscale model integration of multiple biological scales. This is particularly relevant given that MEK inhibitors and PI3K/AKT inhibitors exhibit different molecular properties that influence their tissue penetration and pharmacokinetic profiles 71 , with MEK inhibitors tending to be smaller in size, which can lead to faster penetration. Our model suggests that such a mismatch in diffusion rates could mean that even a nominally simultaneous systemic dose does not result in simultaneous target engagement within the AGS cell population. In fact, if the PI3K/AKT inhibitor reaches the cells first, this would inadvertently create the very PI3K/AKT-first sequence our model identifies as optimal. This insight provides a testable hypothesis for the clinical efficacy of simultaneous dosing: it may work precisely because it achieves the correct sequential action at the cellular level. This has direct translational implications, suggesting that efficacy might be rescued or enhanced by implementing staggered dosing schedules or by altering drug formulations to ensure the PI3K/AKT inhibitor achieves target saturation before MEK is inhibited. The model’s ability to quantitatively predict how a specific head start for a slower drug could overcome a large diffusion mismatch highlights its potential utility as part of systems designed for the prediction of novel, clinically relevant dosing strategies 9 . These predictions exemplify how the multiscale layer bridges the gap between Boolean signaling abstractions and the spatial, pharmacokinetic realities of in vitro and in vivo settings, enabling hypothesis generation that neither approach could produce in isolation. Despite these results, the present work presents some limitations on the technical and scientific sides. The current multiscale model is a simplification of a real cancer, lacking relevant features of the cancer microenvironment such as blood flow, extracellular matrix and immune components. The predictive power of the model is necessarily defined by its underlying assumptions, including the focus on a cell-line-specific Boolean network, which may overlook off-target effects or parallel resistance pathways. Moreover, the multiscale framework inherits the limitations of the Boolean model: drugs that produce experimental effects but no Boolean-level phenotypic change, such as the TAK1 inhibitor in the AGS network, cannot be calibrated within the current pipeline. Additionally, we are restricted to the specific RTCA time-series dataset obtained through personal communication, as there are currently no public databases that systematically archive this type of continuous viability data, limiting our ability to extend the calibration to additional compounds. That said, the framework is modular and could be extended to additional drugs within the same cell line as suitable calibration data become available. For example, dose-response curves from resources such as DepMap could be used to calibrate transfer functions and enable predictions across broader perturbation panels in AGS. Extending to different cell lines, however, would require constructing cell-line-specific Boolean networks and recalibrating phenotypic parameters. The current results should be interpreted within the AGS system and cannot yet be generalized beyond this context without such recalibration and network inference, though emerging tools for automated network construction may help streamline this process in future work 72 . Calibration was performed using population-level data, which may not fully capture the single-cell heterogeneity that drives differential therapeutic responses. Additionally, our study uses a cell line model rather than primary tumor tissue, which inherently implies reduced heterogeneity compared to in vivo systems. Specifically, we cannot capture the full spectrum of subclonal diversity that recent literature identifies as a primary driver of combination efficacy in patient cohorts 73 – 75 . Nevertheless, we note that synergistic interactions are well established in homogeneous cell line systems and can emerge from intracellular signaling interactions 76 . Beyond the biological simplifications, our multiscale framework introduces several methodological assumptions. First, the Boolean formalism uses discrete ON/OFF states that simplify graded interactions. While each agent runs an independent stochastic Boolean simulation introducing cell-level variability, the update schemes impose artificial temporal constraints. Multiscale calibration is therefore essential to adjust parameters like apoptosis rates so outputs match experimental dynamics despite the abstract Boolean timing. Second, we model a homogeneous AGS cell population representing an adenocarcinoma cell line, with interactions simplified to contact inhibition and computational constraints preventing a one-to-one replica of the in vitro system. Third, cross-scale coupling assumes time-scale separation between fast Boolean dynamics and slower phenotypic processes, using transfer functions to map discrete states to continuous rates. Fourth, the model depends on the Boolean network topology and logic rules, which may not capture all regulatory interactions or context-dependent modifications. Finally, to maintain tractability, we use coarse-grained representations of drug diffusion, nutrient gradients, and cell mechanics, omitting details like active transport or heterogeneous ECM that could be incorporated in future work. While these choices enable computational feasibility and interpretability, they constrain the fidelity of certain biological processes. Future work will explore alternative coupling strategies and hybrid formalisms, and sensitivity analyses to systematically assess these assumptions. Future work will advance this framework along three complementary axes. First, we will develop a dynamic network model informed by multi-timepoint data to elucidate temporal feedback loops and condition-specific attractors. Second, we will broaden the calibration beyond the AGS system by incorporating additional cell lines and drugs whenever matched molecular and growth data become available, enabling systematic cross-comparison of how different phenotypes and regulatory networks modulate synergy. Third, we plan to leverage the calibrated multiscale platform to systematically explore how spatial heterogeneity, subclonal diversity, and microenvironmental factors shape therapeutic outcomes. This includes implementing mutation models to study resistance trade-offs through variations in drug uptake kinetics or binding affinities, including experimentally-derived heterogeneity from patient or cohort data, as well as creating a mechanistic learning platform for 3D simulations that captures key cancer microenvironment features, including heterogeneous initial geometries, kinetically-realistic transport processes, specific cell-cell mechanical interactions, and spatially resolved dosing protocols. This will allow us to test hypotheses about resistance emergence, optimal scheduling under realistic drug diffusion constraints, the effect of actual heterogeneous cancer population on drug combination efficacy and the interplay between drug diffusion and cellular organization. Additionally, extending the calibration to full dose-response surfaces rather than single IC50 concentrations would enable more comprehensive validation of the framework’s predictive capacity across concentration ranges. These developments, combined with preclinical validation of scheduling predictions, will create a powerful in silico platform for discovering robust therapeutic strategies with greater clinical translation potential. Methods Drug treatment experimental data We used previously pulished experimental data 34 , in which AGS gastric cancer cells were treated with a panel of kinase inhibitors targeting key oncogenic signaling pathways. We repurposed this data to validate the outputs of our computational simulations. The inhibitors were selected in the original study for their specificity, potency, and relevance to cancer-associated signaling networks and are reported in Supplementary Table 1 . Cell proliferation in the original study was measured in real time using the xCELLigence RTCA system, which records changes in electrical impedance to derive a Cell Index used as a proxy for viable cells. Experimental dose-response curve processing The dose-response curves for the AGS cell line were obtained from the Supplementary Material of ref. 34 . For the purpose of simulation, calibration and validation, we extracted the log-scale drug response curves manually. Each curve was extracted five times to ensure consistency, and the resulting replicates were averaged. We then transformed the data to linear scale, fitted a sigmoidal curve, and computed the Hill coefficient, which was then employed for parameterizing the Hill values for the drug-specific drug perturbation Hill equation (Eq. ( 5 )). Drug response curves were fitted using a four-parameter Hill equation to quantify the concentration-dependent effects of kinase inhibitors on cell proliferation. Experimental data from dose-response experiments were preprocessed by normalizing cell proliferation indices to a 0–1 scale, where 0 represents complete growth inhibition and 1 represents untreated control growth. The drug effect was calculated as the inverse of the normalized cell index (drug effect = 1 − cell index). The Hill equation was implemented as: f ( x ) = bottom + top − bottom 1 + 1 0 EC50 log 1 0 x Hill coefficient 1 where x is the log-transformed drug concentration, EC50 log is the log-transformed concentration producing 50% of the maximum effect, Hill coefficient describes the steepness of the dose-response relationship, and bottom and top represent the minimum and maximum drug effects, respectively. Curve fitting was performed using scipy’s curve_fit function with bounded optimization to ensure biologically plausible parameter ranges: EC50 log ∈ [ x min , x max ] , Hill coefficient ∈ [0.1, 5], bottom ∈ [0, 0.5], and top ∈ [0.5, 1.1]. Initial parameter guesses were set to the mean of the concentration range for EC50 log , 1.0 for the Hill coefficient, and the observed minimum and maximum values for the bottom and top parameters, respectively. The fitted Hill coefficient provides a quantitative measure of the drug’s potency and the steepness of its concentration-response relationship, with higher values indicating a more sensitive response to concentration changes. This parameter was subsequently used to inform the value for this parameter in the Boolean target activation target function (see Eq. ( 5 )), in our multiscale model simulations. Experimental growth time-course data processing Experimental Cell Index time course data were provided through personal communication from Åsmund Flobak. The data consisted of Cell Index measurements recorded by the xCELLigence RTCA system at regular intervals over ~5.5 days, with three technical replicates per condition. We selected drugs that met two conditions: RTCA time-series data were available, and the Boolean model predicted a reduction in proliferation. Drugs such as the TAK1 inhibitor, which showed experimental effects but no Boolean-level response in the AGS Boolean network, were excluded from the calibration because the multiscale model relies on the predictive scope of the underlying molecular model. Raw data processing was necessary to account for experimental artifacts inherent to impedance-based measurements, particularly during the initial cell settling phase and the pre-treatment period. We first averaged the three technical replicates to generate a mean Cell Index value for each time point. To eliminate artifacts from initial cell settling and pre-treatment disturbances, we excluded data between 400 and 2500 min, where the Cell Index measurements typically show artefactual changes in Cell Index computation. The remaining time series data was interpolated using cubic spline interpolation with the SciPy “interp1d" function to generate evenly spaced measurements at 40-min intervals. To enable direct comparison between different experimental conditions and with the simulation data, the Cell Index values were normalized to a 0–100 scale using min-max normalization, where the minimum and maximum values were derived from the wild-type control growth curves. This normalization approach preserved the relative growth dynamics while standardizing the scale across all experimental conditions. The final processed data retained measurements up to 4200 min (70 h), capturing the complete drug response profile. The processed data was used for seven experimental conditions: control (wild-type), single drug treatments (PI103 targeting the PI3K at 0.70 μm, PD0325901 targeting the MEK proteins (MAP2K1/2) at 35.00 nM, and AKTi-1,2 (AKT Inhibitor VIII) targeting AKT (AKT1, AKT2, AKT3) and combination treatments (PD0325901 at 17.50 nM + PI103 at 0.35 μm, and AKTi-1,2 + PD0325901). The processed time course data is provided in the Supplementary Data . Definition and quantification of synergy We assessed drug interactions using Loewe additivity following ref. 34 . This method quantified interactions based on dose equivalence at a pre-specified effect level. For a chosen effect, we first identified each agent’s monotherapy dose that produced that effect. Combination doses were then expressed as fractions of their respective monotherapy-equivalent doses and summed. The summed fraction determined interaction type: exactly 1 indicated additivity, less than 1 indicated synergy, and greater than 1 indicated antagonism. Throughout this work, we use the term “synergy” to refer to drug combinations that were labeled as synergistic by the Loewe framework in ref. 34 . Specifically, we refer PI3Ki-MEKi and AKTi-MEKi combinations. Our multiscale simulations aim to reproduce the population-level outcomes of these known synergistic combinations rather than to compute new synergy metrics, and we evaluate success by comparing predicted growth curves to experimental time course measurements. Multiscale modeling and simulation protocol Agent-based multiscale simulations were performed using PhysiBoSS , an open-source C++ software designed for multi-scale mechanistic modeling of biological systems. PhysiBoSS allows incorporating signaling networks expressed as Boolean networks within each agent in order to model cellular decision-making from the molecular level. The PhysiBoSS simulation consists of two primary components: a microenvironment formed by a 3D mesh where partial differential equations solve the diffusion of various substrates, and a set of discrete entities or agents representing individual cells with encoded rules dictating their behavior and phenotype. Simulations were executed for 4200 simulation minutes, parallelized across 8 threads, using PhysiBoSS version 2.2.0. The model operates across multiple temporal scales through a hierarchical process schedule: substrate diffusion is computed at each 0.01 min, mechanical updates are computed at each 0.1 min, cell phenotype decisions are updated every 6 min. Boolean network has a default update time (intracellular_dt) of 10 min, but as it is a calibrated parameter, in our model it ranges from 1 to 40 min. The 2D simulation domain was of 600 × 600 × 50 (μm) with 20 × 20 × 20 μm voxels, replicating in vitro experimental conditions, simulated for 4200 min (Fig. 1 A). The simulation initializes ~1500 cells in a disk pattern with 300 μm radius and 2.9 μm cell spacing to match the RTCA well-plate configuration used in the in vitro experiments. Cell spacing refers to the center-to-center distance between adjacent cells, computed as a multiple of the cell radius, and defines the initial spatial arrangement of the simulation. For the 3D scenario experiments for optimizing drug dosage timing, we built a 3D cell setup of 225 × 225 × 700 μm, with 20 × 20 × 20 μm voxels and simulated them for 5000 min. The initial setup was a box of cells with dimensions 225 × 225 × 100 μm, filled with a cell spacing of 2.86 μm. Simulation output was stored at 40-min intervals, emulating experimental data acquisition intervals. The AGS PhysiBoSS model is provided in https://github.com/bsc-life/ags_synergy_paper , along with a simple tutorial for running the model, as well as the different initial configurations for replicating single-drug calibrated experiments, the synergies, and the 3D pharmacodynamic heterogeneity assays. Boolean model simulations are performed using the MaBoSS algorithm with default parameter values. MaBoSS is executed in each agent every 10 min of simulation time by default, and 50 steps are executed. Discrete time is deactivated, and the pseudorandom seed is set to 100. The scaling setting is set to 1, and inheritance of the Boolean model state is activated, which transfers the information of the Boolean simulation state from a cell to its daughter cell after division. The probability of the initial node state OFF was 1 for all nodes except for AKT, MEK, PI3K, TAK1, betacatenin, GSK3, and p38alpha, scaling and the frequency for MaBoSS execution ( i n t r a c e l l u l a r _ d t ) were calibrated for each single drug, as these change when the drug affects the Boolean network. In addition to the signaling pathways, the Boolean model includes two functional multivalued readouts that represent opposing cell fate outcomes: (a) Prosurvival, indicating growth-promoting signals, governed by the states of the TCF, RSK, and cMYC nodes; and (b) Antisurvival, indicating pro-apoptotic signaling, determined by the activation of Caspase8, Caspase9, and FOXO nodes. These readouts are multivalued with continuous ranges between 0 and 1, where 0 represents complete inhibition and 1 represents full commitment to the cellular response. Importantly, these two readouts are not mutually exclusive, allowing the model to capture intermediate phenotypic states between proliferation and apoptosis. As a result, different single or combined drug perturbations dynamically shift the balance of prosurvival and antisurvival signaling. Although the model incorporates multivalued logic and thus extends beyond classical Boolean networks, we follow established field conventions in referring to it as a ’Boolean’ model. The states of the Prosurvival nodes are coupled to the agent-based model’s growth model, while the Antisurvival nodes are linked to the agent’s death model. The specific implementation of this coupling between Boolean model readouts and the agent-based growth and death mechanisms is described in more detail in the following sections. Therefore, without any drug perturbations, all Prosurvival readouts are set to 1, and all Antisurvival readouts are set to 0, so there is maximum growth at the AGS doubling time and the baseline apoptosis rate. Simple diffusion model for drug transport We model drug internalization as a process of passive diffusion across the cell membrane, consistent with Fick’s second law 77 . The drug flux is governed by a drug-specific permeability coefficient, with the selected range of values supported by literature (Table 1 ). This approach assumes that the net flux of a drug across the cell membrane is proportional to the difference in concentration between the cell interior and the extracellular environment, modulated by the membrane permeability and the cell surface area. J = P ⋅ S ⋅ ( C int − C ext ) 2 This diffusion mechanism is implemented as an ordinary differential equation (ODE) that updates the intracellular and extracellular drug concentrations at each diffusion time step (0.01 min), allowing the system to dynamically track drug distribution across the cellular and extracellular compartments. The net flux J (in attomoles per minute, a m o l / m i n ) of a drug across the cell membrane is given by Equation ( 2 ), where P is the permeability coefficient (converted to μm/min), S is the cell surface area (μm 2 ), C int is the intracellular drug concentration (mM), and C ext is the extracellular drug concentration (mM). Permeability is typically provided in cm/s. We convert it to μ m/min using: P μ m / min = P cm / s × 60 × 1 0 4 3 Intracellular concentration is computed from the total internalized substrate (in a m o l ) divided by the cell volume (in μm 3 ), with appropriate conversion factors to yield mM. The following pseudocode summarizes the ODE-based update for each cell and drug: For each cell: For each drug: // Retrieve parameters permeability = cell-specific permeability (cm/s) cell_radius = cell radius (μm) cell_surface = 4 * π * (cell_radius) 2 (μm 2 ) cell_volume = cell volume (μm 3 ) total_internalized = total internalized drug (amol) external_concentration = extracellular drug concentration (mM) // Convert permeability to μ m/min permeability_um_min = permeability * 60 * 1e4 // Compute internal concentration (mM) internal_concentration = total_internalized / cell_volume // Ensure non-negative concentration if internal_concentration < 0: internal_concentration = 0 // Compute diffusion constant (μm 3 /min) diffusion_constant = permeability_um_min * cell_surface // Compute concentration difference (mM) concentration_difference = internal_concentration - external_concentration // Convert concentration difference to mol/μm 3 concentration_difference_mol_um3 = concentration_difference * 1e-18 // Compute flux (mol/min) flux = diffusion_constant * concentration_difference_mol_um3 // Convert flux to amol/min flux_amol_min = flux * 1e18 // Update cell's net export rate for this drug cell.net_export_rate[drug] = flux_amol_min The model assumes passive diffusion only, with no active transport or binding. The flux is defined as positive when the net movement is from inside to outside the cell, and negative when the net movement is into the cell. This ODE is integrated at each diffusion time step (0.01 min) to update the internal and external drug concentrations. Agent-level contact inhibition model AGS cell growth rate is affected by the exerted pressure from neighbor agents, allowing modeling of confluence. The pressure affects growth through a Hill equation, where contact inhibition parameters ( k A p , H p ) modulate the growth response to pressure from neighboring cells, and the basal growth rate ( r basal γ ) determines the maximum proliferation rate in the absence of pressure. calibrated on the control condition experimental growth curve. To ensure that baseline cellular proliferation dynamics accurately reproduce experimental observations, control curve calibration was performed using both CMA-ES and GA algorithms. These algorithms were employed to calibrate the 5-dimensional parameter space encompassing “AGS Baseline Behavior Parameters" (Table 1 ): initial spatial configuration ( R tumor , cell_spacing), pressure-based contact inhibition parameters ( k A p , H p ), and basal growth rate ( r basal γ ). Internalized drug to Boolean network perturbation transfer functions The model implements a multiscale approach to drug pharmacodynamics linking extracellular drug concentration to intracellular target inhibition. Drug molecules diffuse through the microenvironment according to Fick’s laws and enter cells via permeability coefficients calibrated for each inhibitor (PI3Ki, MEKi, and AKTi) based on literature 39 – 44 . Diffusion coefficients were set to 600 μm 2 /min based on literature 45 , 46 . Once internalized, drugs reversibly bind to their molecular targets following mass action kinetics. The temporal dynamics of drug-target complex formation are governed by: d [ D T ] i d t = k 1 ( [ D ] i − [ D T ] i ) ( [ T ] − [ D T ] i ) − k − 1 [ D T ] i 4 where [ D ] i represents intracellular drug concentration, [ T ] is the total target concentration, and [ D T ] i is the drug-target complex concentration, with mass conservation maintaining free concentrations. The computed drug-target complex concentration serves as input to a Hill function that converts continuous binding levels into discrete Boolean network perturbations: P i n a c t i v a t i o n = [ D T ] H K 1 2 + [ D T ] H 5 This probability determines stochastic inactivation of the corresponding target node (PI3K, MEK, or AKT) during Boolean network simulation. Drug-target mappings correspond to: PI103 (PI3K), PD0325901 (MEK), and AKTi-1,2 (AKT). Half-maximal inhibition values ( K 1 2 ) were derived from experimental IC 50 data, while Hill coefficients ( H ) were obtained from dose-response curves 34 . The model incorporates three specific kinase inhibitors with their corresponding molecular targets: PI103 (PI3K inhibitor), PD0325901 (MEK inhibitor), and AKTi-1,2 (AKT inhibitor). Half-maximal inhibition values ( K 1 2 ) for the Hill function (Eq. ( 5 )) were fixed based on experimental IC 50 values as detailed in Supplementary Table 1 . Hill coefficients ( H ) were derived from dose-response curves presented in the supplementary material of ref. 34 . At any given diffusion timestep (0.01 min), mass conservation principles ensure that the free drug concentration is [ D ] f r e e = [ D ] i − [ D T ] i and the free target concentration is [ T ] f r e e = [ T ] − [ D T ] i , where total concentrations remain constant and only the distribution between free and bound forms changes according to the binding kinetics described by Eq. ( 4 ). The stochastic inactivation process uses the probability P i n a c t i v a t i o n calculated from Eq. ( 5 ) to determine whether the corresponding target node becomes inactive during each Boolean network simulation time step, enabling direct coupling between continuous drug-target binding dynamics and discrete Boolean network states. Boolean network to cellular phenotype transfer functions The Boolean network state is translated into cellular phenotypes through transfer functions that map network readout states to growth and apoptosis rates. Proliferative signaling is computed as S γ = ∑ i ∈ A w i ⋅ x i where x i is the Boolean state of node i , w i its associated weight. The weights are constrained such that ∑ i ∈ A w i = 1 summed over the set of relevant growth-promoting nodes (cMYC, TCF, RSK), ensuring S γ lies within the interval [0, 1]. This weighted signal is then mapped to proliferation rate via a Hill function: r γ = S γ H γ K A , γ + S γ H γ 6 Similarly, pro-apoptotic signals are computed as S α from anti-survival nodes (FOXO, Caspase8, Caspase9) and mapped to apoptosis rates: r α = S α α × r α , max K A , α + S α α + r b a s a l , α 7 Proliferation rates range from baseline doubling rate (when S γ = 1) to full growth arrest, while apoptosis rates vary from basal value (when S α = 0) to maximum r α ,max (when S α = 1). Phenotypic rate changes are implemented gradually through linear convergence controlled by response parameters r response, γ and r response, α , which define the slope of the linear transition from the current phenotypic state to the target state specified by the S γ and S α parameters. This gradual implementation prevents instantaneous phenotypic transitions and ensures smooth adaptation to changing cellular conditions. All transfer function parameters were calibrated for each experimental condition using the EMEWS optimization framework. Complete parameter specifications and calibration ranges are provided in Supplementary Tables 4, 5 . Model calibration with GA and CMA-ES optimization algorithms We conducted parallel model exploration and high-throughput hypothesis testing using the EMEWS framework (Extreme-scale Model Exploration with Swift) 25 . This high-performance computing workflow, built on the Swift/T parallel scripting language 78 , enables seamless integration of model exploration algorithms with parallel simulation evaluation. Building upon our previous PhysiBoSS-EMEWS workflow 20 , we implemented three distinct parameter calibration approaches: a uniform parameter sweep for the parallel evaluation of candidate sets, and two evolutionary algorithms: a Genetic Algorithm (GA) and Covariance Matrix Adaptation Evolution Strategy (CMA-ES). We implemented both evolutionary algorithms using the Python DEAP package 79 . The GA implementation utilized the “eaSimple” function from DEAP 79 , which represents an evolutionary algorithm 80 . The algorithm iteratively performs the following steps for each generation: (1) Runs the PhysiBoSS experiment for each individual and obtains the objective metric (RMSE); (2) Selects the next generation individuals by sampling from the current population based on the metric; (3) Applies crossover and mutation operators on the best or Hall of Fame individuals to produce offspring; (4) Evaluates the new individuals; and (5) Replaces the current population with the offspring. This process continues until convergence on a minimum RMSE value is achieved. In our PhysiBoSS-EMEWS DEAP implementation, we used tournament selection with a size of 3 individuals, uniform crossover with a 0.5 probability, and a custom mutation function with a 0.2 probability. For CMA-ES, we employed the DEAP “Strategy” module, which implements the basic parameters of this algorithm 58 . All CMA-ES parameters were kept at their default values, with sigma set to 1. The CMA-ES implementation uses the “eaGenerateUpdate” function in DEAP, which iteratively performs the following steps: (1) Runs the PhysiBoSS simulation for each individual in the population; (2) Evaluates the fitness of each individual; (3) Updates the strategy parameters, dynamically adjusting the step size and covariance matrix; (4) Generates a new population; and (5) Evaluates the new fitness scores. This process continues until convergence is reached. The developed model includes many different parameters detailed in Table 1 , with some having experimentally estimated ranges from literature (maximum apoptosis rate, doubling time, drug permeability) while others required calibration through optimization-via-simulation. We employed both Genetic Algorithm (GA) and Covariance Matrix Adaptation Evolution Strategy (CMA-ES) evolutionary algorithms to calibrate unknown parameters by fitting the experimental growth curves from control and single drug treatment conditions (PI3Ki, MEKi, AKTi). The fitness function minimized a composite Root Mean Square Error (RMSE) metric between experimental and simulated growth curves, with both curves normalized between 0 and 1 using min-max scaling. To ensure robust fitting quality, the base RMSE was enhanced with three additional criteria: differences between experimental and simulated means of the last 20 timepoints, final timepoint differences, and a cell death penalty for simulations with no recorded apoptotic cells. RMSE calculations focused on the period from drug addition (1280 simulation minutes) until experiment completion (4200 simulation minutes). We first calibrated control curve parameters using a 5-dimensional parameter space comprising AGS Baseline Behavior Parameters (Table 1 ): initial spatial configuration ( R p o p u l a t i o n , cell spacing), pressure-based contact inhibition parameters ( k A , p , H p ), and basal growth rate ( r b a s a l , γ ) constrained to match AGS doubling time of 20–24 h 35 – 38 . Subsequently, for single-drug fitting, we explored an 18-dimensional parameter space comprising 3 drug-specific “Drug Transport and Binding Parameters" ( k d r u g , X , K 1, d r u g , X , K −1, d r u g , X ) and 15 common parameters including “Cellular Response Parameters" governing growth ( H γ , K A , γ , r γ , r e s p o n s e ) and apoptosis ( H α , K A , α , r m a x , α , r r e s p o n s e , α ), and “Boolean Network Integration Parameters" with node weights for growth ( ω γ , c M Y C , ω γ , T C F , ω γ , R S K ) and apoptosis ( ω α , F O X O , ω α , C a s p 8 , ω α , C a s p 9 ), plus Boolean simulation parameters ( d t i n t r a c e l l u l a r , s c a l i n g i n t r a c e l l u l a r ). The maximum apoptosis rate ( r m a x , α ) was constrained between 24-48 hours based on experimental reports 47 – 49 . Both algorithms were run for a maximum of 50 generations or until convergence, with an initial population of 300 individuals and 4 replicates per individual (See Supplementary Fig. 1 ). Synergy parameter space construction To investigate drug synergies in the 2D simulation setting, a consensus parameter space approach was developed where each single drug (PI3Ki, MEKi, AKTi) was independently calibrated using CMA-ES. The 10% of parameter sets with lowest RMSE values were identified from each calibration experiment. To identify parameters that show consistent behavior across different drug inhibition conditions, we performed a statistical analysis of the top 10% parameter distributions from independent calibrations. For each parameter, we calculated the union and intersection of parameter ranges across all conditions, the percentage of overlap between distributions (ratio of intersection to union range), and the mean coefficient of variation (CV) across conditions. A parameter was considered a consensus parameter if it showed high overlap percentage (>65%) between the three drug conditions and consistent behavior reflected by a relatively low mean CV. A consensus set of the 15 common parameters was extracted from these top-performing parameter sets to serve as the foundation for defining the parameter space of the synergy. Two distinct synergy parameter spaces were defined for PI3Ki-MEKi and AKTi-MEKi combinations by incorporating the consensus common parameters with respective drug-specific parameters from each combination. To define the parameter space for synergy experiments, we employed a union-based approach where for each parameter p , we calculate its range as: Range p = [ min ( μ 1 − 3 σ 1 , μ 2 − 3 σ 2 ) , max ( μ 1 + 3 σ 1 , μ 2 + 3 σ 2 ) ] , where μ i and σ i represent the mean and standard deviation of parameter p in each single-drug experiment i . Both common and drug-specific parameters were compared between constituent drugs, with the union of their value ranges used to define parameter bounds for each specific synergy space, ensuring mechanistic consistency between constituent drugs. Drug-specific parameters maintained their original ranges obtained from the top 10% best individuals for single-drug calibration experiments. Each synergy space encompassed a 21-dimensional parameter domain comprising 15 common parameters (Boolean-to-phenotype transfer function parameters ( H γ , K A , γ , H α , K A , α ); Boolean network readout weights ω γ , c M Y C , ω γ , T C F , ω γ , R S K and ω α , F O X O , ω α , C a s p 8 , ω α , C a s p 9 ; maximum apoptosis rate r m a x , α ; and Boolean simulation parameters) plus 6 drug-specific transport and binding parameters (3 for each drug: k d r u g , X , K 1, d r u g , X , K −1, d r u g , X from Table 1 ). Within each synergy space, uniform sweeps of 5000 parameter combinations were performed, simulating combination treatments (drugs at half their individual I C 50 concentrations) and corresponding single-drug controls for systematic synergy assessment. The top 10% best-performing individuals from each synergy space were identified for subsequent analysis. Parameter distributions between PI3K+MEK and AKT+MEK synergy conditions were compared by calculating the mean and standard deviation in each condition, effect size (Cohen’s d) between conditions, and percentage of overlap between parameter ranges. Effect sizes were classified as negligible ( d < 0.2), small (0.2 ≤ d < 0.5), medium (0.5 ≤ d < 0.8), or large ( d ≥ 0.8). Distribution overlap was categorized as high (>65%), medium (35–65%), or low (<35%). 3D drug administration experiments We constructed a 3D simulation environment (225 × 225 × 700 μm) with 20 μm³ voxels, initializing cells in a 100 μm Z-axis layer with 2.86 μm spacing, and simulated for 5000 min. For these experiments, we selected the top 100 best individuals from each synergy parameter sweep and gathered the mean value for each parameter (Fig. 3 ). The experimental design systematically explored two parameter spaces: a 2-dimensional space including diffusion coefficients of both drugs, and an extended 4-dimensional space incorporating drug addition timing for each synergistic drug pair. To explore the effects of spatiotemporal heterogeneity in drug diffusion, different timing scenarios were tested by varying the interval ( Δ t ) between drug administrations, ranging from −5000 to 5000 min (negative values indicate drug Y administered before X, positive values indicate X administered before Y). Drugs were administered at half their I C 50 concentrations (Supplementary Table 2 ) and sustained throughout the simulation. Diffusion coefficient combinations were tested at 6, 60, 600, and 6000 μm 2 /min for both drugs, generating 16 experimental conditions with 5 replicates each (Supplementary Figs. 14, 15 ). Control experiments included single-drug treatments at each diffusion coefficient and drug-free negative controls, with single-drug controls administered at time zero to establish baseline effects independent of timing variations (Supplementary Fig. 12 ). These were used as baselines to validate the synergy in our 3D simulation scenario by quantifying the efficacy of the combination (% of alive cells in comparison to the no-drug control condition). This approach enabled valid synergy comparisons across different administration schedules while maintaining a consistent reference point. Machine learning-enhanced sensitivity analysis We performed a global sensitivity analysis using Random Forest (RF) regression combined with SHAP (SHapley Additive exPlanations) value analysis to quantify parameter importance across the explored parameter space. RF regression models were trained as computationally efficient surrogate models of our multiscale PhysiBoSS model. Training datasets comprised calibrated parameter sets (Table 1 ) and their corresponding RMSE fitness scores from CMA-ES calibration experiments for each condition: control (20,920 sets), MEKi (41,899 sets), AKTi (45,000 sets), PI3Ki (45,000 sets), PI3Ki-MEKi (45,000 sets), and AKTi-MEKi (45,000 sets). Parameter sets served as input features and RMSE values as target variables. For each treatment condition, RF regressors were implemented using scikit-learn (v1.2.2) with random subsampling of 1000 parameter sets to ensure computational tractability. Models were constructed with 42 decision trees and training was parallelized across available CPU cores. SHAP values were computed using the trained RF surrogate models to quantify parameter importance across the explored parameter space. For each treatment condition, SHAP values were calculated for all parameters across the subsampled datasets (1000 parameter sets per condition). Parameters were ranked by their overall influence on RMSE output based on the distribution of SHAP values. SHAP analysis was implemented using the shap Python package (v0.41.0), with TreeExplainer used for efficient computation with the RF models. Supplementary information Supplementary Information (17.8MB, pdf) Supplementary data 1 (4.3KB, csv) Supplementary data 2 (4.3KB, csv) Supplementary data 3 (4.3KB, csv) Supplementary data 4 (2.5KB, csv) Supplementary data 5 (4.4KB, csv) Supplementary data 6 (2.5KB, csv) Acknowledgements We thank Åsmund Flobak for providing the Boolean network model, the time-course experimental data, and for his valuable input and feedback on the manuscript. This work was supported by the European Union’s Horizon Europe Research and Innovation programme projects HANAMI (grant agreement No. 101136269), by the Horizon 2020 projects CREXDATA (grant agreement No. 101092749) and PerMedCoE (grant agreement No. 951773) and by the ERA PerMed project Oncologics (ERAPERMED2020-036). O.H.M. was supported by a predoctoral “AGAUR-FI Joan Oró” fellowship from “Secretaria d’Universitats i Recerca del Departament de Recerca i Universitats de la Generaliat de Catalunya i del Fons Europeu Social Plus” (2025 FI-1 01021). Author contributions Conceptualisation: O.H.M., M.P., A.M., and A.V.; methodology: O.H.M., M.P. and A.M.; software: O.H.M., M.P. and A.M.; in silico experiments and simulations: O.H.M.; results analysis: O.H.M. and M.P.; writing original draft preparation: O.H.M. and M.P.; writing review and editing: O.H.M., M.P., A.M., and A.V.; visualisation: O.H.M. and M.P.; supervision: M.P., A.M., A.V. All authors have read and agreed to the published version of the manuscript. Data availability The complete model calibration results dataset generated in this study has been deposited in the Zenodo repository under accession number 10.5281/zenodo.17702426. Processed experimental growth timecourse curves are included in the Supplementary Data. Code availability The multiscale PhysiBoSS model was implemented using PhysiCell (version 1.12.0) integrated with PhysiBoSS addons and the MaBoSS Boolean network simulation engine. The complete, self-contained multiscale model including all configuration files, scripts and tutorials is available at https://github.com/bsc-life/ags_synergy_paper . The analysis scripts were implemented in Python and require the following third-party packages: pandas for data manipulation and analysis, numpy for numerical computing, scipy for scientific computing (including optimization, interpolation, and statistical functions), matplotlib and seaborn for data visualization, scikit-learn for machine learning algorithms (Random Forest regression, feature importance analysis, and model evaluation), shap for SHapley Additive exPlanations to quantify parameter importance, SALib for global sensitivity analysis (Sobol indices), joblib for parallel processing, statsmodels for advanced statistical modeling, pcdl (PhysiCell Data Loader) and pctk for reading PhysiBoSS simulation outputs, hillfit and neutcurve for drug response curve fitting, and vtk (Visualization Toolkit) for converting simulation outputs to ParaView-compatible formats. All scripts were developed and tested using Python 3. Competing interests The authors declare no competing interests. Footnotes Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. These authors contributed equally: Alfonso Valencia, Miguel Ponce-de-Leon. Supplementary information The online version contains supplementary material available at 10.1038/s41540-026-00669-4. References 1. Jin, H., Wang, L. & Bernards, R. Rational combinations of targeted cancer therapies: Background, advances and challenges. Nat. Rev. Drug Discov. 22 , 213–234 (2023). [ DOI ] [ PubMed ] [ Google Scholar ] 2. Fitzgerald, J. B., Schoeberl, B., Nielsen, U. B. & Sorger, P. K. Systems biology and combination therapy in the quest for clinical efficacy. Nat. Chem. Biol. 2 , 458–466 (2006). [ DOI ] [ PubMed ] [ Google Scholar ] 3. Pritchard, J. R., Lauffenburger, D. A. & Hemann, M. T. Understanding resistance to combination chemotherapy. Drug Resist. Updates 15 , 249–257 (2012). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 4. Konieczkowski, D. J., Johannessen, C. M. & Garraway, L. A. A convergence-based framework for cancer drug resistance. Cancer cell 33 , 801–815 (2018). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 5. Dhanyamraju, P. K., Schell, T. D., Amin, S. & Robertson, G. P. Drug-tolerant persister cells in cancer therapy resistance. Cancer Res. 82 , 2503–2514 (2022). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 6. Mokhtari, R. B. et al. Combination therapy in combating cancer. Oncotarget 8 , 38022–38043 (2017). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 7. Alber, M. et al. Integrating machine learning and multiscale modeling–perspectives, challenges, and opportunities in the biological, biomedical, and behavioral sciences. npj Digital Med. 2 , 115 (2019). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 8. Zhang, J. D., Sach-Peltason, L., Kramer, C., Wang, K. & Ebeling, M. Multiscale modelling of drug mechanism and safety. Drug Discov. Today 25 , 519–534 (2020). [ DOI ] [ PubMed ] [ Google Scholar ] 9. Yankeelov, T. E. et al. Multi-scale Modeling in Clinical Oncology: Opportunities and Barriers to Success. Ann. Biomed. Eng. 44 , 2626–2641 (2016). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 10. Metzcar, J., Pletz, K., Rocha, H. L. & Rozum, J. C. Translating and evaluating single-cell Boolean network interventions in the multiscale setting. arxiv http://arxiv.org/abs/2501.16052 . 2501.16052. (2025). 11. An, G. et al. Optimization and Control of Agent-Based Models in Biology: A Perspective. Bull. Math. Biol. 79 , 63 (2017). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 12. Ghaffarizadeh, A., Heiland, R., Friedman, S. H., Mumenthaler, S. M. & Macklin, P. PhysiCell: An open source physics-based cell simulator for 3-D multicellular systems. PLOS Comput. Biol. 14 , e1005991 (2018). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 13. Montagud, A., Ponce-de Leon, M. & Valencia, A. Systems biology at the giga-scale: Large multiscale models of complex, heterogeneous multicellular systems. Curr. Opin. Syst. Biol. 28 , 100385 (2021). [ Google Scholar ] 14. Rocha, H. L. et al. A hybrid three-scale model of tumor growth. Math. Models Methods Appl. Sci. 28 , 61–93 (2018). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 15. Albert, I., Thakar, J., Li, S., Zhang, R. & Albert, R. Boolean network simulations for life scientists. Source Code Biol. Med. 3 , 16 (2008). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 16. Bloomingdale, P., Nguyen, V. A., Niu, J. & Mager, D. E. Boolean Network Modeling in Systems Pharmacology. J. Pharmacokinetics Pharmacodyn. 45 , 159–180 (2018). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 17. Calzone, L., Barillot, E. & Zinovyev, A. Logical versus kinetic modeling of biological networks: Applications in cancer research. Curr. Opin. Chem. Eng. 21 , 22–31 (2018). [ Google Scholar ] 18. Abou-Jaoudé, W. et al. Logical Modeling and Dynamical Analysis of Cellular Networks. Front. Genet. 7 , 94 (2016). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 19. Kauffman, S. A. Metabolic stability and epigenesis in randomly constructed genetic nets. J. Theor. Biol. 22 , 437–467 (1969). [ DOI ] [ PubMed ] [ Google Scholar ] 20. Ponce-de Leon, M. et al. Optimizing Dosage-Specific Treatments in a Multi-Scale Model of a Tumor Growth. Front. Mol. Biosci. 9 , 836794 (2022). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 21. Ponce-de Leon, M. et al. PhysiBoSS 2.0: A sustainable integration of stochastic Boolean and agent-based modelling frameworks. npj Syst. Biol. Appl. 9 , 1–12 (2023). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 22. Gutenkunst, R. N. et al. Universally Sloppy Parameter Sensitivities in Systems Biology Models. PLOS Comput. Biol. 3 , e189 (2007). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 23. Cockrell, C. & An, G. Utilizing the Heterogeneity of Clinical Data for Model Refinement and Rule Discovery Through the Application of Genetic Algorithms to Calibrate a High-Dimensional Agent-Based Model of Systemic Inflammation. Front. Physiol. 12 , 662845 (2021). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 24. Joslyn, L. R., Kirschner, D. E. & Linderman, J. J. CaliPro: A Calibration Protocol That Utilizes Parameter Density Estimation to Explore Parameter Space and Calibrate Complex Biological Models. Cell. Mol. Bioeng. 14 , 31–47 (2020). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 25. Ozik, J. et al. High-throughput cancer hypothesis testing with an integrated PhysiCell-EMEWS workflow. BMC Bioinforma. 19 , 483 (2018). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 26. Ozik, J., Collier, N., Heiland, R., An, G. & Macklin, P. Learning-accelerated discovery of immune-tumour interactions. Mol. Syst. Des. Eng. 4 , 747–760 (2019). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 27. Ruscone, M. et al. Building multiscale models with PhysiBoSS, an agent-based modeling tool. Brief. Bioinforma. 25 , bbae509 (2024). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 28. Horvath, P. et al. Screening out irrelevant cell-based models of disease. Nat. Rev. Drug Discov. 15 , 751–769 (2016). [ DOI ] [ PubMed ] [ Google Scholar ] 29. Foucquier, J. & Guedj, M. Analysis of drug combinations: Current methodological landscape. Pharmacol. Res. Perspect. 3 , e00149 (2015). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 30. Minchinton, A. I. & Tannock, I. F. Drug penetration in solid tumours. Nat. Rev. Cancer 6 , 583–592 (2006). [ DOI ] [ PubMed ] [ Google Scholar ] 31. Trédan, O., Galmarini, C. M., Patel, K. & Tannock, I. F. Drug Resistance and the Solid Tumor Microenvironment. JNCI J. Natl. Cancer Inst. 99 , 1441–1454 (2007). [ DOI ] [ PubMed ] [ Google Scholar ] 32. Ballesta, A., Zhou, Q., Zhang, X., Lv, H. & Gallo, J. Multiscale Design of Cell-Type–Specific Pharmacokinetic/Pharmacodynamic Models for Personalized Medicine: Application to Temozolomide in Brain Tumors. CPT Pharmacomet. Syst. Pharmacol. 3 , 112 (2014). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 33. Li, X. L., Oduola, W. O., Qian, L. & Dougherty, E. R. Integrating Multiscale Modeling with Drug Effects for Cancer Treatment. Cancer Inform. 14 , 21–31 (2016). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 34. Flobak, A. et al. Discovery of Drug Synergies in Gastric Cancer Cells Predicted by Logical Modeling. PLOS Comput. Biol. 11 , e1004426 (2015). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 35. Cowley, G. S. et al. Parallel genome-scale loss of function screens in 216 cancer cell lines for the identification of context-specific genetic dependencies. Sci. Data 1 , 140035 (2014). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 36. Kim, H. J. et al. Forty-nine gastric cancer cell lines with integrative genomic profiling for development of c-MET inhibitor. Int. J. Cancer 143 , 151–159 (2018). [ DOI ] [ PubMed ] [ Google Scholar ] 37. Matozaki, T. et al. Missense mutations and a deletion of the p53 gene in human gastric cancer. Biochem. Biophys. Res. Commun. 182 , 215–223 (1992). [ DOI ] [ PubMed ] [ Google Scholar ] 38. Niapour, A. & Seyedasli, N. Acquisition of paclitaxel resistance modulates the biological traits of gastric cancer AGS cells and facilitates epithelial to mesenchymal transition and angiogenesis. Naunyn Schmiedebergs. Arch. Pharmacol. 395 , 515–533 (2022). [ DOI ] [ PubMed ] [ Google Scholar ] 39. Barnett, S. F. et al. Identification and characterization of pleckstrin-homology-domain-dependent and isoenzyme-specific Akt inhibitors. Biochem. J. 385 , 399–408 (2005). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 40. Knight, Z. A. et al. A Pharmacological Map of the PI3-K Family Defines a Role for p110 α in Insulin Signaling. Cell 125 , 733 (2006). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 41. Lindsley, C. W. et al. Allosteric Akt (PKB) inhibitors: Discovery and SAR of isozyme selective inhibitors. Bioorg. Med. Chem. Lett. 15 , 761–764 (2005). [ DOI ] [ PubMed ] [ Google Scholar ] 42. Raynaud, F. I. et al. Biological properties of potent inhibitors of class I phosphatidylinositide 3-kinases: From PI-103 through PI-540, PI-620 to the oral agent GDC-0941. Mol. Cancer Therapeutics 8 , 1725–1738 (2009). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 43. Sebolt-Leopold, J. S. et al. Blockade of the MAP kinase pathway suppresses growth of colon tumors in vivo. Nat. Med. 5 , 810–816 (1999). [ DOI ] [ PubMed ] [ Google Scholar ] 44. Thompson, N. & Lyons, J. Recent progress in targeting the Raf/MEK/ERK pathway with inhibitors in cancer drug discovery. Curr. Opin. Pharmacol. 5 , 350–356 (2005). [ DOI ] [ PubMed ] [ Google Scholar ] 45. Baish, J. W. et al. Scaling rules for diffusive drug delivery in tumor and normal tissues. Proc. Natl. Acad. Sci. 108 , 1799–1803 (2011). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 46. Tunggal, J. K., Cowan, D. S., Shaikh, H. & Tannock, I. F. Penetration of anticancer drugs through solid tissue: A factor that limits the effectiveness of chemotherapy for solid tumors. Clin. Cancer Res. 5 , 1583–1586 (1999). [ PubMed ] [ Google Scholar ] 47. Li, Y., Liu, Y., Shi, F., Cheng, L. & She, J. Knockdown of Rap1b Enhances Apoptosis and Autophagy in Gastric Cancer Cells via the PI3K/Akt/mTOR Pathway. Oncol. Res. 24 , 287–293 (2016). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] [ Retracted ] 48. Wang, T.-E., Wang, Y.-K., Jin, J., Xu, B.-L. & Chen, X.-G. A novel derivative of quinazoline, WYK431 induces G2/M phase arrest and apoptosis in human gastric cancer BGC823 cells through the PI3K/Akt pathway. Int. J. Oncol. 45 , 771–781 (2014). [ DOI ] [ PubMed ] [ Google Scholar ] 49. Will, M. et al. Rapid Induction of Apoptosis by PI3K Inhibitors Is Dependent upon Their Transient Inhibition of RAS–ERK Signaling. Cancer Discov. 4 , 334–347 (2014). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 50. Gross, A. et al. Representing dynamic biological networks with multi-scale probabilistic models. Commun. Biol. 2 , 1–12 (2019). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 51. Fumiã, H. F. & Martins, M. L. Boolean Network Model for Cancer Pathways: Predicting Carcinogenesis and Targeted Therapy Outcomes. PLOS ONE 8 , e69008 (2013). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 52. Engelman, J. A. et al. Effective use of PI3K and MEK inhibitors to treat mutant Kras G12D and PIK3CA H1047R murine lung cancers. Nat. Med. 14 , 1351–1356 (2008). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 53. Roberts, P. J. et al. Combined PI3K/mTOR and MEK Inhibition Provides Broad Antitumor Activity in Faithful Murine Cancer Models. Clin. Cancer Res. 18 , 5290–5303 (2012). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 54. Kwong, L. N. & Davies, M. A. Navigating the Therapeutic Complexity of PI3K Pathway Inhibition in Melanoma. Clin. Cancer Res. 19 , 5310–5319 (2013). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 55. Chandarlapaty, S. Negative feedback and adaptive resistance to the targeted therapy of cancer. Cancer Discov. 2 , 311–319 (2012). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 56. Shimizu, T. et al. The Clinical Effect of the Dual-Targeting Strategy Involving PI3K/AKT/mTOR and RAS/MEK/ERK Pathways in Patients with Advanced Cancer. Clin. Cancer Res. 18 , 2316–2325 (2012). [ DOI ] [ PubMed ] [ Google Scholar ] 57. Moles, C. G., Mendes, P. & Banga, J. R. Parameter estimation in biochemical pathways: a comparison of global optimization methods. Genome Res. 13 , 2467–2474 (2003). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 58. Hansen, N. & Ostermeier, A. Completely derandomized self-adaptation in evolution strategies. Evolut. Comput. 9 , 159–195 (2001). [ DOI ] [ PubMed ] [ Google Scholar ] 59. Igel, C., Hansen, N. & Roth, S. Covariance matrix adaptation for multi-objective optimization. Evolut. Comput. 15 , 1–28 (2007). [ DOI ] [ PubMed ] [ Google Scholar ] 60. Jędrzejewski-Szmek, Z., Abrahao, K. P., Jędrzejewska-Szmek, J., Lovinger, D. M. & Blackwell, K. T. Parameter Optimization Using Covariance Matrix Adaptation–Evolutionary Strategy (CMA-ES), an Approach to Investigate Differences in Channel Properties Between Neuron Subtypes. Front. Neuroinformatics 12 , 47 (2018). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 61. Walpole, J., Papin, J. A. & Peirce, S. M. Multiscale computational models of complex biological systems. Annu. Rev. Biomed. Eng. 15 , 137–154 (2013). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 62. Guo, Z., Hou, D. & He, Q. Hybrid Genetic Algorithm and CMA-ES Optimization for RNN-Based Chemical Compound Classification. Mathematics 12 , 1684 (2024). [ Google Scholar ] 63. Potter, D. S. et al. Dynamic BH3 profiling identifies pro-apoptotic drug combinations for the treatment of malignant pleural mesothelioma. Nat. Commun. 14 , 2897 (2023). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 64. Bhola, P. D. et al. High-throughput dynamic BH3 profiling may quickly and accurately predict effective therapies in solid tumors. Sci. Signal. 13 , eaay1451 (2020). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 65. Potter, D. S., Du, R., Bhola, P., Bueno, R. & Letai, A. Dynamic BH3 profiling identifies active BH3 mimetic combinations in non-small cell lung cancer. Cell Death Dis. 12 , 741 (2021). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 66. Miyoshi, S. et al. Antitumor activity of MEK and PI3K inhibitors against malignant pleural mesothelioma cells in vitro and in vivo. Int. J. Oncol. 41 , 449–456 (2012). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 67. Peck, B., Ferber, E. C. & Schulze, A. Antagonism between FOXO and MYC Regulates Cellular Powerhouse. Front. Oncol. 3 , 96 (2013). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 68. Almhanna, K. et al. MK-2206, an Akt inhibitor, enhances carboplatinum/paclitaxel efficacy in gastric cancer cell lines. Cancer Biol. Ther. 14 , 932 (2013). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 69. Pistritto, G., Trisciuoglio, D., Ceci, C., Garufi, A. & D’Orazi, G. Apoptosis as anticancer mechanism: Function and dysfunction of its modulators and targeted therapeutic strategies. Aging 8 , 603–619 (2016). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 70. Metzcar, J., Jutzeler, C. R., Macklin, P., Kohn-Luque, A. & Bruningk, S. C. A review of mechanistic learning in mathematical oncology. Front. Immunol. 15 , 1363144 (2024). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 71. Jokinen, E. & Koivunen, J. MEK and PI3K inhibition in solid tumors: Rationale and evidence to date. Therapeutic Adv. Med. Oncol. 7 , 170–180 (2015). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 72. Ruscone, M., Vazquez, M. & Valencia, A. Intelligent Tool Orchestration for Rapid Mechanistic Model Prototyping: MCP Servers as AI-Biology Interfaces. bioRxiv , https://www.biorxiv.org/content/10.1101/2025.09.10.675105v1 (2025). 73. Palmer, A. C., Chidley, C. & Sorger, P. K. A curative combination cancer therapy achieves high fractional cell killing through low cross-resistance and drug additivity. eLife 8 , e50036 (2019). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 74. Gambardella, G. et al. A single-cell analysis of breast cancer cell lines to study tumour heterogeneity and drug response. Nat. Commun. 13 , 1714 (2022). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 75. Huang, X. et al. Single-cell systems pharmacology identifies development-driven drug response and combination therapy in B cell acute lymphoblastic leukemia. Cancer Cell 42 , 552–567.e6 (2024). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 76. Menden, M. P. et al. Community assessment to advance computational prediction of cancer drug combinations in a pharmacogenomic screen. Nat. Commun. 10 , 2674 (2019). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 77. Lipinski, C. A., Lombardo, F., Dominy, B. W. & Feeney, P. J. Experimental and computational approaches to estimate solubility and permeability in drug discovery and development settings. Adv. Drug Deliv. Rev. 23 , 3–25 (1997). [ DOI ] [ PubMed ] [ Google Scholar ] 78. Wozniak, J. M. et al. Swift/T: Scalable data flow programming for many-task applications. In Proceedings of the 18th ACM SIGPLAN Symposium on Principles and Practice of Parallel Programming (PPoPP ’13) 309–310, 10.1145/2442516.2442559 (2013). 79. Fortin, F.-A., De Rainville, F.-M., Gardner, M.-A., Parizeau, M. & Gagné, C. DEAP: Evolutionary algorithms made easy. J. Mach. Learn. Res. 13 , 2171–2175 (2012). [ Google Scholar ] 80. Baeck, T., Fogel, D. B. & Michalewicz, Z. (eds.) Evolutionary Computation 1: Basic Algorithms and Operators (CRC Press, 2018). Associated Data This section collects any data citations, data availability statements, or supplementary materials included in this article. Supplementary Materials Supplementary Information (17.8MB, pdf) Supplementary data 1 (4.3KB, csv) Supplementary data 2 (4.3KB, csv) Supplementary data 3 (4.3KB, csv) Supplementary data 4 (2.5KB, csv) Supplementary data 5 (4.4KB, csv) Supplementary data 6 (2.5KB, csv) Data Availability Statement The complete model calibration results dataset generated in this study has been deposited in the Zenodo repository under accession number 10.5281/zenodo.17702426. Processed experimental growth timecourse curves are included in the Supplementary Data. The multiscale PhysiBoSS model was implemented using PhysiCell (version 1.12.0) integrated with PhysiBoSS addons and the MaBoSS Boolean network simulation engine. The complete, self-contained multiscale model including all configuration files, scripts and tutorials is available at https://github.com/bsc-life/ags_synergy_paper . The analysis scripts were implemented in Python and require the following third-party packages: pandas for data manipulation and analysis, numpy for numerical computing, scipy for scientific computing (including optimization, interpolation, and statistical functions), matplotlib and seaborn for data visualization, scikit-learn for machine learning algorithms (Random Forest regression, feature importance analysis, and model evaluation), shap for SHapley Additive exPlanations to quantify parameter importance, SALib for global sensitivity analysis (Sobol indices), joblib for parallel processing, statsmodels for advanced statistical modeling, pcdl (PhysiCell Data Loader) and pctk for reading PhysiBoSS simulation outputs, hillfit and neutcurve for drug response curve fitting, and vtk (Visualization Toolkit) for converting simulation outputs to ParaView-compatible formats. All scripts were developed and tested using Python 3. Articles from NPJ Systems Biology and Applications are provided here courtesy of Nature Publishing Group ACTIONS View on publisher site PDF (3.3 MB) Cite Collections Permalink PERMALINK Copy RESOURCES Similar articles Cited by other articles Links to NCBI Databases Cite Copy Download .nbib .nbib Format: AMA APA MLA NLM Add to Collections Create a new collection Add to an existing collection Name your collection * Choose a collection Unable to load your collection due to an error Please try again Add Cancel Follow NCBI NCBI on X (formerly known as Twitter) NCBI on Facebook NCBI on LinkedIn NCBI on GitHub NCBI RSS feed Connect with NLM NLM on X (formerly known as Twitter) NLM on Facebook NLM on YouTube National Library of Medicine 8600 Rockville Pike Bethesda, MD 20894 Web Policies FOIA HHS Vulnerability Disclosure Help Accessibility Careers NLM NIH HHS USA.gov Back to Top