Step-by-step comparison of ordinary differential equation and agent-based approaches to pharmacokinetic-pharmacodynamic models.

Mathematical models in oncology aid in the design of drugs and understanding of their mechanisms of action by simulation of drug biodistribution, drug effects, and interaction between tumor and healthy cells. The traditional approach in pharmacometrics is to develop and validate ordinary differential equation models to quantify trends at the population level. In this approach, time-course of biological measurements is modeled continuously, assuming a homogenous population. Another approach, agent-based models, focuses on the behavior and fate of biological entities at the individual level, which subsequently could be summarized to reflect the population level. Heterogeneous cell populations and discrete events are simulated, and spatial distribution can be incorporated. In this tutorial, an agent-based model is presented and compared to an ordinary differential equation model for a tumor efficacy model inhibiting the pERK pathway. We highlight strengths, weaknesses, and opportunities of each approach.

INTRODUCTION

Understanding the dose–concentration–effect relationship is fundamental in drug development and mathematical models are instrumental to aid data interpretation. Pharmacokinetic‐pharmacodynamic (PKPD) models can support the prediction of the right dose. The PK part describes the body’s impact on the molecular drug concentration by investigating the relationship between dose and plasma concentrations to gain information on drug absorption, distribution, metabolism, and excretion. Active drugs and their metabolites are measured in an accessible biologic fluid, such as blood, plasma, or urine. From these assays, concentration‐time curves of the plasma and tissue drug concentration are obtained, and PK measures and parameters, such as area under the curve, maximum concentration, clearance, volume of distribution, and elimination half‐life, are determined. PK studies also consider intrinsic and extrinsic factors that influence individual and population systemic exposure. Critical intrinsic factors include gender, ethnicity, genetic polymorphism, and renal or hepatic organ dysfunction, and important extrinsic factors are drug‐drug interactions and concomitant administration of drug products with food. , The PD part describes the impact of the drug molecule on the body (i.e., receptor binding, receptor sensitivity, post‐receptor effects, and chemical interactions). It includes all pharmacological and pathophysiological responses, and therapeutical effect following administration of a drug or placebo. PD studies are designed to inform a drug’s mechanism of action and the dose‐response relationship. Response can be expressed as a direct or indirect measure of efficacy and/or safety of the drug using either biomarkers, surrogate endpoints, or clinical endpoints to quantify the PD effect.

Mathematical PKPD models can be generated and approached in different ways. In this work, we apply agent‐based models (ABMs) and compare its properties with the more traditional ordinary differential equation (ODE) approach.

ODEs are often used to model the change of drug concentration and effect over time. As a method where behavior for each entity of the system is regulated on a population level, ODEs are suitable to model phenomena that are centrally coordinated, for example, by a set of mass transfer binding reactions. Examples are concentration gradients (e.g., nutrient, oxygen gradients, or PK and PD models). Data for this modeling approach may come from the periphery, the tissue, or at the whole organism level. Solutions of differential equations typically describe the mean behavior of the system over time. Variability can be introduced in one or more parameters, as used in the nonlinear mixed‐effects approach. Complicated models that include not only changes in time but also space, for example, spatial differences in drug effect or concentration, use partial differential equations (PDEs) that depend on several variables and contain partial derivatives. ,

Agent‐based modeling

Biological phenomena have heterogeneous characteristics (e.g., cancer with regions of hypoxia, necrosis, quiescence, and proliferation). In addition, cancer growth depends on the changing micro‐environment (e.g., oxygen, glucose, and pH gradients). Drug resistance can manifest as a consequence of the limited diffusion of an effective amount of the drugs far into the tumor core due to irregular vascular structure. In such settings, time, space but also the characteristics of each individual cell or group of cells and their interaction with the surrounding environment are important for modeling drug effects.

ABMs simulate heterogeneity with one or a set of different agents that have attributes and act autonomously in an environment according to certain rules. To apply this concept to PKPD modeling, cells or a group of cells can be represented as agents, whereas the drug amount is simulated with Boolean values, ODEs, PDEs, or as agents. Attributes of the cell agents could be location, mutation rate, growth rate, drug resistance, or antigenicity. Agents can interact in an environment which can consist of different conditions (i.e., spatial location, concentration of drugs, nutrients, receptors, or other agents). Interactions with the environment and neighboring agents are based on intracellular decision‐making rules that, for example, describe the behavior of tumor cells when interacting with immune cells or the tumor‐microenvironment. An agent undergoes growth, proliferation, quiescence, apoptosis, or necrosis as a response to surrounding environmental conditions or interaction with other agents. As such, ABM is considered as a method coordinated on the cellular level with an emergent behavior from cell‐cell and cell‐moieties interactions that is not centrally coordinated by a set of mass transfer binding reactions. It can handle complex biological phenomena and local phenomena driven by discrete decision. Agents act only upon local information on the state of the system, rather than being affected by the global system state. In this way, ABMs are well‐suited to represent the transition of emergent behavior between one scale of organization (i.e., microscopic scale) to behavior observed at another (i.e., macroscopic scale). , Simple activities at the microscopic scale can cause complex behavior on the macroscopic scale. , , , , For example, mutation at the cellular level can create patterns of growth and extent of metastasis at the tumoral and stromal level. Unexplained patterns can be discovered, such as circumstances where the tumor cells are completely eliminated by the immune system, without the need for any cancer therapies. Adaption of the agents to changing circumstances is possible upon changing environmental conditions. Agents can have memory where information about the environment is stored and future behavior can be modified according to past stages. This can be used in modeling priming of T cells or memory T cells. Heterogeneity is generated by the variability in the attributes of each agent (i.e., mutation, resistance, antigenicity, or number of DNA breaks after radiation), the decision rules, and the stochastic process governing those.

Every agent may have a spatial distinct position. Lattice‐based models restrict positions to a fixed grid; off‐lattice models do without this restriction. To create a lattice‐based model, regular structured grids (Cartesian [2D or 3D], dodecahedral [3D]), or unstructured grids are used. , Lattice‐based methods can be further classified by their spatial resolution, meaning how many lattice sites are used for creating cells. Cellular automaton models contain at most a single cell in one lattice site. Discrete lattice‐based rules are used to update a single cell at every time step. To reduce grid artifacts, lattice sites are updated in random order. A cell can remain, move to an unoccupied adjacent lattice site (randomly or by a directional stimulus such as chemo‐ or haptotaxis), die and vacate a lattice site, or divide and allocate a progenitor cell in a neighboring site. The number of neighbor cells depend on the shape of the grid. In a square lattice, cells have four or eight neighbors (von Neumann or Moore neighborhood) or cells can be surrounded by six symmetrically located cells in a hexagonal grid. , Lattice gas cellular automaton (LGCA) models contain multiple cells in one lattice. Instead of tracking every single cell, LGCA traces a group of cells that moves through channels from individual lattice sites. Cellular Potts models utilize multiple lattice sites to simulate one cell.

Off‐lattice models can be categorized into center‐based models that focus on cell volumes and boundary‐based methods that model cell boundaries. Center‐based methods can be differentiated into three approaches: center‐based models, subcellular element models, and clusters. Center‐based models simulate each cell’s center of mass or volume by using one agent per cell. Cell positions are updated after exchanges of adhesive, repulsive, locomotive, and drag‐like forces between cell centers. Subcellular element models display cell morphology in greater detail by having multiple agents for subcellular elements of each cell. A different option is clusters, where an agent simulates a cluster of cells or functional units (e.g., breast glands or colon crypts). Boundary‐based methods include vertex‐based models or front‐tracking methods. In vertex‐based models, forces on the vertices are computed, which is useful in modeling confluent tissues. Front‐tracking methods, such as the immersed boundary method, can be used for greater spatial resolution.

Programming ABMs

An environment suitable for beginners is provided by the software Netlogo, in which models can be created and executed, and the results visualized. ABMs can be implemented in object‐oriented programming languages, such as Python, Java, C, or C++. Differences between those languages are user friendliness, speed, and available libraries. Python is beginner friendly and has a rich environment of standard libraries whereas Java, C, and C++ are often faster. An advantage of coding the model is that the researcher has control over all aspects of the model, and additional features can be easily implemented or changed. Disadvantages are the time spent programming and the need for programming skills. Because many models use the same or similar building blocks with small variations, libraries for different programming languages are available which provide frameworks and templates that allow users to design a customized model. Examples are MASON for Java, SWARM for Java and Objected C, Repast for Java and C++, Chaste as a source code in C++, and MESA for Python. A benefit is the reduced model construction time. Disadvantages could be the time spend to learn the functions of the library and that the needed feature might not be directly available. , ,

EXAMPLE: PKPD‐ABM AND PKPD‐ODE MODEL FOR SIMULATING ANTI‐CANCER TREATMENT WITH COBIMETINIB

This tutorial aims to introduce ABM, which is widely used in social sciences and systems biology. , It gives a simple example of a PKPD‐ABM of anti‐cancer treatment with the MEK inhibitor cobimetinib, targeting the RAF/MEK/ERK pathway, which plays an important role in cell proliferation and survival. This signaling pathway is initiated by binding of growth factors, cytokines, and extracellular mitogens, which activates receptor tyrosine kinases and leads to a reaction cascade resulting in activation of MEK and phosphorylation of ERK. As a consequence, cellular responses are cell proliferation and survival. The RAF/MEK/ERK pathway is frequently mutated in cancer cells leading to increased proliferation and survival. Cobimetinib aims to counteract this by decreasing phosphorylation and as a consequence decrease cell proliferation and cell survival. ,

The ABM simulates the behavior of tumor cells as agents during treatment with cobimetinib. The environment is given by the time course and effect of the administered drug. To keep this model simple, no cell location is introduced. Cell death and division occur on the microscopic scale and a summary of cell behavior provides the total number of tumor cells on the macroscopic scale. The model is implemented in the Python programming language. The code can be found at https://github.com/VanThuyTruong/Tutorial. A PKPD model previously published by Wong et al. implemented in ODE, was used to simulate the cobimetinib concentration time course and the phosphorylation of ERK. The model is used here as an illustrative example and implemented as described in the original reference, for the sole purpose of numerical simulation. To account for the PK, a hybrid approach is considered, with PK being governed by an ODE system and combined with a tumor ABM to determine the interactions of cells in the presence of drug concentrations. To compare the behavior of the hybrid PKPD‐ABM tumor model, a PKPD‐ODE tumor model is implemented to highlight the relative strengths and limitations of each method.

PK model

The PKs of the drug was characterized by an ODE model that links dose taken orally to a plasma compartment and a tumor micro‐environment compartment (see part a in Figure 1). The PKs in the tumor compartment induces a reduction of phosphorylation of the ERK pathway (part b) that is used as the driving force to implement the PKPD model, either as a full ODE model (part c1) or as a hybrid ODE‐ABM (part c2).

FIGURE 1

Model structure. (a) The PK model describes tumor disposition of cobimetinib. (b) The PKPD model links the percentage pERK decrease, d(t), to the amount of cobimetinib inside the tumor compartment. Phosphorylated ERK plays an important role in cell division. Therefore, pERK could be seen as a biomarker of tumor growth and a decrease of pERK causes a decrease in cell division. (c1) the ODE model simulates the effect of the pERK decrease on the number of tumor cells in the population. k(t) is the pERK value inside one population. The mean number of cells increases with birth rate λ and decreases with death rate µ. A high pERK value favors birth, a low value favors decay over growth. (c2) An ABM, where tumor cell death or division is driven by the percentage pERK decrease caused by the amount of cobimetinib inside the tumor compartment, d(t), and the individual pERK value of each tumor cell, k(t). (a–c1) together result in a PKPD‐ODE and (a–c2) together result in a PKPD‐ABM. ABM, agent‐based model; IC50, half‐maximal inhibitory concentration; ODE, ordinary differential equation; PKPD, pharmacokinetic pharmacodynamic

We use the model by Wong et al. to describe kinetics of drug concentration in oral, plasma, and tumor compartments after a weight‐based oral dose of cobimetinib. The ODEs are linear, with absorption rate k 0 from the oral to the plasma compartment, clearance rate β from the plasma, and inter‐compartmental rate γ from the plasma to the tumor compartment. Transfer from the tumor to the plasma is represented by the inter‐compartmental rate α. Because the drug is a small molecule, it is assumed that drug excess is high and loss due to target binding, catabolism, or elimination in the tumor compartment is negligible. Therefore, no elimination rate is introduced in the tumor compartment. The time courses follow the following ODEs:

The initial condition X 0(0) is the given dose (1, 3, or 10 mg/kg of body weight) converted to µmoles/kg. (The molecular weight of cobimetinib is 531.3 g/mol.) The initial drug amount in the plasma compartment X(0) and in the tumor compartment X(0) is zero. Our parameter values are based on those of Wong et al.

PD model

The PD model describes the effect of cobimetinib by the following relationship between X(t) and the percentage pERK decrease, d(t): where IC50 is the cobimetinib amount in the tumor compartment at which the percentage pERK decrease is half of I max, the maximum percentage pERK decrease, and h is the Hill coefficient.

Depending on the initial pERK value, k(0), and the percentage pERK decrease caused by the drug d(t) at time t, the pERK value is:

In Figure 1, we illustrate the two types of tumor model, based on Equation 5, that we will consider: agent‐based (c2) and population‐based (c1). In the ABM, we assign a different initial value k(0) to each tumor cell, drawn uniformly in the range (0, k max). Values of k(0) greater than 100 represent cells that have activating mutations in B‐RAF resulting in higher activation of MEK. Individual values of k(t) evolve according to Equation 5. In the population‐based model, a single value k(t) represents the average pERK status of the population of tumor cells. That is, cell‐to‐cell variability is introduced in the agent‐based model; the total number of cells is an integer that depends on time because individual cells may die or divide according to pERK levels. The population‐based model, on the other hand, is a set of ODEs where one value of k(t) represents the influence of the drug on the population of tumor cells in the ODE. The population size is also governed by an ODE, which we now introduce.

PKPD‐ODE

The population‐based model is the PKPD‐ODE model shown as c1 in Figure 1. A single pERK value k(t) ∈ (0, k max) characterizes the tumor cell population and is updated according to Equation 5.

The size of the tumor cell population at time t, T(t), obeys:

The constant λ > 0 models the mean birth rate. Death is modelled as a decreasing function of k(t). We choose µ = 2λ, which has the effect that birth and death balance when and facilitates the comparison with ABM. We use the value λ = 0.0828 day−1 of Wong et al.

PKPD‐ABM

In the ABM, each tumor cell is impacted by the pERK decrease individually. One agent is a tumor cell i with a value k(t) as an attribute, representing its degree of intracellular phosphorylation of ERK, determined by the amount of drug in the tumor compartment X(t) according to Equation 3 from the PK model. As the drug concentration lowers, the cell’s pERK value returns toward its initial value k(0), unless the cell dies according to the decision‐making algorithm.

Decision‐making algorithm (microscopic scale)

The fate of a cell is governed by a set of rules depending on its ERK phosphorylation level (see Figure 1 part c2). pERK thresholds determine the behavior of the cells. A cell whose pERK value is above the division threshold is said to be in the division pool; a pERK value below the death threshold consigns the cell to the death pool. Inside the division or death pool, a division or death event happens with a certain probability per unit time. In between these thresholds, a cell remains quiescent.

The ABM is a Markov process that can be summarized as follows. At time t, the number of tumor cells is N(t), each cell with its pERK value. The probability that cell i divides in the interval (t,t + ∆t) is λ(t)∆t where

Similarly, the probability that cell i dies in the interval (t,t + ∆t) is µ(t)∆t where

We use the Gillespie algorithm, which is a fair method for designating which of these 2N(t) possible events is the first, after t, to occur. We construct the following sums:

We decide which type of event occurs, then decide which cell it happens to. With probability Λ(t)/(M(t) + Λ(t)), the first event after time t is a division event. In that case, the probability that, of the N(t) possible cells, it is cell i that divides, is λ(t)/Λ(t). With probability M(t)/(M(t) + Λ(t)), the first event after time t is a death event. In that case, the probability that cell i dies is µ(t)/M(t). The initial pERK value of the parent cell at the time of division is inherited by the daughter cells.

An example of this decision‐making algorithm is given in the Supplementary Material.

Simulation outcomes

Different simulation scenarios were explored. The drug dose and the distribution of initial pERK were varied for each simulation. The ODE model shows the behavior at the population level with the overall cell number (see for example Figure 2.1d where each blue line represents one tumor population) and the homogeneous pERK value of each population (see for example Figure 2.1b where each purple line describes the pERK value inside one population). The ABM gives further insight in each cell of the population with the pERK value of each cell over time (see for example Figure 2.1c where each color represents one cell) and the number of cells in the division, quiescence, and death pool (see for example Figure 2.2b). It is possible to track the pERK value of each cell inside the population over time in Figure 2.1c. A cell can die and cause the line to end or give birth to an offspring which creates a new line starting from the time‐point of birth (marked with a cross). The crosses describe the discrete time points for events caused by the Gillespie algorithm. The scatter plots show each cell inside the population with its pERK value at specific time points (see Figure 2.2c–f) and videos at https://github.com/VanThuyTruong/Tutorial/tree/main/videos show the all scatter‐plots with cell birth and death over the simulation period.

FIGURE 2.1

Simulation of the ODE model and ABM with a single oral dose of 3 mg/kg. pERK values are initially uniformly distributed. The parameter values are given in Table S1. Figure a shows the PKPD model, common to the ODE model and ABM. Figure b and d show multiple trajectories of the PKPD‐ODE model, each with the same initial cell population size but a different initial pERK value. Figure c shows individual pERK values in one realization of the ABM Here, each cell has a different initial pERK value, chosen from the uniform distribution in (0, 200). Figure e shows the total cell numbers in 100 such realizations.

FIGURE 2.2

One realization of the ABM. In Figure 2.2c– f , one dot represents the pERK value of one cell at one timepoint. The ABM may provide a more realistic model because it captures heterogeneity, different scales, and emergent behavior. On the other hand, ODE is suitable for modeling well‐mixed compartments with mass transfer and simple interactions at one scale level. A video of the scatter plots can be found at https://github.com/VanThuyTruong/Tutorial/blob/main/videos/3mgkg%20single%20dose%20uniform%20distribution.mp4 ABM, agent‐based model; ODE, ordinary differential equation; PKPD, pharmacokinetic pharmacodynamic

Single dose treatment with 3 mg/kg

First, a treatment with a single dose is simulated with 3 mg/kg per body weight (see Figures 2.1 and 2.2). The initial pERK value for the ODE model and ABM are sampled from a uniform distribution in the range 0−200. The ODE model simulates 100 tumor cell populations with different initial pERK (see Figure 2.1b), whereas the ABM follows the fate of 100 tumor cells for a given patient (see Figure 2.1c). In both models exposure‐driven pERK reduction induces tumor shrinkage followed by regrowth as the drug level falls (see Figure 2.1b–e). A difference is that the ABM is able to track cell death and division of each cell (Figure 2.1c).

Simulating 100 times with the ABM (Figure 2.1e) yields 100 different population histories because the initial pERK distribution in each trial is different due to the stochasticity of the Gillespie algorithm. We observe that the total number of cells in a population may be higher in comparison to the ODE model, because the cells with a higher pERK values are more likely to survive and divide (see Figure 2.1d and 2.1e).

Figure 2.1c, 2.2a–f provide a closer look into one population. We see that the pERK value of a majority of cells falls under the death threshold (see Figure 2.2d), those cells go into the death pool (see Figure 2.2b) and therefore can die, which decreases the number of cells (see Figure 2.2a). When the drug level falls, the cells recover and the pERK value increases (see Figure 2.2e). After ~ 70 hours, the drug concentration is low enough for the cells to gain a pERK value above the division threshold and enter the division pool, those cells cause a population growth (see Figures 2.2a, 2.2b, 2.2f). This simulation shows that cells with a low pERK level. which are assumed not being mutated, are most affected by the drug while drug treatment leads to an evolution of mutated cells with a high pERK value.

Bimodal distribution of pERK values

The next example aims to recapitulate a tumor where cells have a distribution of pERK levels based on two given phenotypes. We construct bimodal distributions from two normal distributions where 50 cells in the ABM or 50 populations in the ODE model have a pERK value with a mean of 60% pERK, whereas 50 cells in the ABM or 50 populations in the ODE model have a pERK value with a mean of 190% pERK. The standard deviations are 10% in both distributions. A single dose of 3 mg/kg was chosen, hence the PKPD of the examples in Figures 2 and 3 are identical. The ODE model shows distinct behavior for populations with different pERK values (see Figure 3.1c), although the ABM does not display the bimodality on the population level (Figure 3.1d). In the ABM, the heterogeneity is shown on the microscopic level (see additional Figure S3.2c in the Supplementary Material) and the cell number is summed up for the macroscopic level (population level). Therefore, the bimodality is only seen at the microscopic level. Taking a closer look into the example, we see that for the ABM, the number of tumor cells initially falls (Figure 3.1d). This is caused by the cells of the population with the mean pERK of 60%, which are either in the death or quiescence pool (see Figures S3.2c, S3.2d in the Supplementary Material). They die out or remain in the quiescence pool and therefore are not able to contribute to the population growth (see in the supplementary Figures S3.1b, S3.2c‐S3.2f ). The cells with high pERK value around 190% are resilient to the drug treatment, because the decrease of active pERK caused by the drug is not enough to bring those cells under the death threshold. Then, the drug amount in the tumor compartment decreases and these tumor cells leave the quiescence pool and enter the division pool where they can divide again and cause an increase of the total cell number (see Figure S3.2b, Figure S3.2e). The tumor cells start dividing and pass their pERK value at the division time‐point to the daughter cells. Looking at Figure S3.2f, we see that the offsprings have initially lower pERK values than the parent cell, because the pERK decrease caused by the drug immediately affects them. During the time course of the simulation, their pERK value increases until they regain their inherited pERK values (Figures 3.1b, S3.2f). Because of the heritage of the pERK value, a typical pattern evolves and in the time course of the simulation there will be a tumor population with a high metabolic phosphorylation pool (compare Figures S3.2c and S3.2f). The overall population growth is driven by the cells with a high pERK value (see Figures 3.1b and 3.1d). Therefore, comparing the overall cell number of 100 simulations with a uniform pERK distribution (Figure 2.1e) and the bimodal distribution (Figure 3.1d) we see similar behavior with a constant population growth after a slight decrease until ~ 50 h. This simulation emphasizes that tumor growth is driven by mutated cells with a high pERK value. This can be seen due to the heterogeneity.

FIGURE 3.1

Simulation of a single dose of 3 mg/kg. The pERK values are initially bimodally distributed. Figures a and c show the PKPD‐ODE, Figures b and d show the PKPD‐ABM as comparison. Additional simulation of the behavior of one population in the PKPD‐ABM can be found in the supplementary. A video of the scatter plots can be found at https://github.com/VanThuyTruong/Tutorial/blob/main/videos/bimodal%20pERK.mp4 ABM, agent‐based model; ODE, ordinary differential equation; PKPD, pharmacokinetic pharmacodynamic

Higher division and death threshold

The effect of changing birth and death thresholds is explored in Figure 4. The parameters are the same as in Figure 2 except for the division threshold of 150% pERK instead of 100% pERK, and a death threshold of 100% pERK instead of 50% pERK. This is compared to an ODE model (Figures 2.1b and 2.1d) and an ABM with the death and division threshold of 50% and 100% pERK (see Figures 2.1c and 2.1e). A treatment of a single dose of 3 mg/kg is given.

FIGURE 4

The effect of changing birth and death thresholds is explored with a simulation of a single dose of 3 mg/kg. The pERK values are initially uniformly distributed. Figures a and b show the PKPD‐ABM with a death threshold of 100% pERK and a division threshold of 150% pERK, as comparison. Figures c–d show the behavior of one population in the PKPD‐ABM. Scatter plots for this simulation are in the supplementary (Figures S4e–h). A video of the scatter plots can be found at https://github.com/VanThuyTruong/Tutorial/blob/main/videos/150%20division%20threshold.mp4. ABM, agent‐based model; PKPD, pharmacokinetic pharmacodynamic

The higher division and death thresholds lead to a longer residence of the tumor cells in the quiescence and death pools (see Figure 4d, and additional simulation in the Figures S4e‐S4h). Despite the drug washout over time, the majority of cells are not able to reach the division pool (Figure S4g), which lowers the growth of the tumor population.

Multiple treatment cycles

A treatment of 1 mg/kg administered every 24 h with 10 cycles was simulated in Figure 5.1 with the ODE and ABM to derive pERK and tumor cell dynamic time‐course. The pERK‐time course in the ABM has less regular onset‐offset granularity due to the time steps of the Gillespie algorithm (Figure 5.1c). The time points are sparse because the time steps in the Gillespie algorithm depend on the number of cells in the death and division pools and most cells are trapped in the quiescence pool (see Figure S5.2b). Overall, the ODE and ABM show similar behavior. As seen in the Supplementary Material, cells oscillate between division and quiescence pools (Figures S5.2d,e). This results in a decrease of cell number. There is a peak every 24 h in the division pool, followed by a peak in the death pool (see time points every 24 h in Figure S5.2b). Upon chronic treatment, the pERK level within tumor cells slowly drops, causing the total cell number to decrease and tumor shrinkage (Figures S5.2b and S5.2a). After ~ 200 hours, the cells in the death pool have all died, whereas some cells with initially high pERK values remain dormant in the quiescence pool (see Figure S5.2b). This resembles drug resistance of mutated cells with a high pERK value.

FIGURE 5.1

Simulation of a multiple treatment cycles of 1 mg/kg every 24 h. The pERK values are initially uniformly distributed. Figure a shows the PKPD model, Figures b and d show the PKPD‐ODE model, Figures c and e show the PKPD‐ABM as comparison. Figure S5.2 in the supplementary show the behavior of one population in the PKPD‐ABM (model c2 in Figure 1). A video of the scatter plots can be found at https://github.com/VanThuyTruong/Tutorial/blob/main/videos/multiple%20cycles%201mgkg.mp4. ABM, agent‐based model; ODE, ordinary differential equation; PKPD, pharmacokinetic pharmacodynamic

DISCUSSION AND CONCLUSIONS

Table 1 compares ABMs and ODE models, using properties (scale, dynamics of interactions, population, space, memory, and stochastic model), strengths and limitations, and implementation (model building and qualification, communication, and applicability). A key feature of ABMs is heterogeneity at the cellular/microscopic level: every agent has an individual set of attributes. Emergent behavior may arise, without central coordination, from local phenomena driven by discrete decisions and interactions between individual agents. Adaptation of the agents to environmental changes is possible with rules set by the modeler, enabling memory to be stored to modify future behavior. Agents may also have a geographical location. , , , Due to their properties, ABMs are naturally suited to the diverse tumor microenvironment, where interactions between different cell types play critical roles in cancer development, progression, and control. The agents in our anticancer treatment example are tumor cells, each with a different level of pERK. Depending on a cell’s current pERK value, it is susceptible to either death or division, governed by stochastic processes. Tumor cells with high pERK values divide and pass their pERK values onto their offspring, whereas the cells with a lower pERK value die or stay quiescent. Thus, the average pERK value in a cell population may increase.

TABLE 1

Comparison between ODE and ABM features and implementation principles for the model‐based approach used to characterize PK and PKPD properties

Comparator	ODE	ABM
Scale	Macroscopic Mean behavior at system level	Microscopic Individual behaviors at cellular level driving emergent behavior at the system level
Dynamics of interaction	Mass transfer dictated by stoichiometric equilibrium and compartmentalization of the system	Rule‐based individual agent interaction and stochasticity
Population	Homogenous	Heterogenous
Space	Not typically implemented	Typically implemented
Memory	Not typically implemented	Typically implemented
Stochastic model	Between subject variability and residual unexplained variability Implemented at the population level by assuming parameter probabilistic distributions	Implicit feature of single agent and governed by stochasticity No distributional assumption
Model building	Rigorous statistical framework for model selection Data‐driven	Hypothesis generation and hypothesis‐testing iterative learning Simulation‐based and/or data calibration depending on empirical evidence
Model qualification	Simulation‐based diagnostic (visual predictive check) Data‐based model qualification (goodness of fit)	Model calibration based on single cell data (in vivo/in vitro/ex vivo experiments) or any data source of relevance for the biological system of interest (micro‐ or macro‐level)
Limitations	Oversimplification Structural rigidity (e.g., compartmentalization) Scalability	Overparameterization Model discrimination Uncertainty in outcomes
Strengths	Well established modeling framework Simple implementation	Emergent behavior to find plausible mechanisms for unforeseen outcomes (e.g., resistance or necrosis) Easier to scale
Computational resources	Typically not a limitation unless large number of differential equations required	Complex ABMs demand high computational power and cause long running times
Model comparison	Straightforward due to similar model structure and model discrimination criteria	More complicated than ODE due to multiplicity of rules, choice of attributes and stochasticity driving emergent behaviors
Communication	Challenging due to theoretical concept with mass transfer and binding kinetics (mathematical knowledge required) Familiar concept in PKPD with large example pool Well‐established with regulatory agencies	Easier to communicate than ODE with non‐modeler audience as more biologically interpretable Not extensively used by PKPD modelers and regulators More challenging to defend as less data‐driven than ODE due to paucity of data at the subscale level
Applicability	Modeling of well‐mixed compartments with mass transfer and simple interactions at one scale level PK models, PD models, and traditional pharmacometric models of exposure‐response Quantitative system pharmacology models with known or observable macro‐level outcomes	Simulation of complex biological systems with subscale components (atomic, molecular, cellular, tissue, organism) System biology models and quantitative system pharmacology models with limited empirical data and most relevant to elucidate unexpected behaviors (micro‐known to produce macro‐unknown)

Abbreviations: ABM, agent‐based model; ODE, ordinary differential equation; PKPD, pharmacokinetic pharmacodynamic.

Other examples of ABMs include the model of Cockrell and Axelrod, simulating the heterogeneous cell types and their proliferation kinetics inside human colon crypts, where quiescent stem cells are at the bottom of the crypt, proliferating cells are close to the bottom third, and differentiated cells are placed in the top two thirds. Proliferating cells are killed by cytotoxic drugs, whereas quiescent stem cells are resistant to cytotoxic drugs due to their low probability of dividing. The on‐lattice ABMs of Kather et al. were used to predict survival, and guide new strategies for immunotherapy, in a colorectal cancer patient cohort. Effects of chemotherapy, immunotherapies, and cell migration inhibitors alone and in combination were simulated.

In ODE models, individual agents and individual interactions are not explicitly considered. Instead, macro‐level outputs are driven by mass transfer dictated by stoichiometric equilibrium and compartmentalization of the system (e.g., blood and tumor compartments). In our ODE model, all the cells have the same pERK value and constitutes a homogenous population; death or division of individual cells is not recorded. ODEs may be used to represent a subset of outcomes following therapeutic interventions at the macro‐level. , Elements of stochasticity may be introduced by assuming that parameters are drawn from probability distributions (to model between‐subject variability) and that measurements are subject to noise. In ODE models, adaptation to changing circumstances may be prespecified in the parameter settings and initial conditions. In most circumstances, ODE models do not have inherent memory features and do not possess a trace of the input variables within the system unless PDE are implemented to localize the interactions. ,

Although ABMs take many different forms, all ODE models have the same mathematical structure, common to many areas of science. Algorithms for numerical solutions of ODEs are widely available and understood, as are methods of data fitting and calibration against observed data. , ODE models can have a simply descriptive purpose, or both a predictive and a descriptive purpose. The simplest models, with the fewest parameters, provide an approximation whose accuracy relies on stochasticity not being preponderant.

In principle, ODE models are based on population‐level assumptions and data, whereas ABMs are constructed from understanding at the level of individual agents. For example, in vitro assays or immunohistochemical data may inform a model. When data for macroscopic phenomena are not available, an ABM can be built as a hypothesis‐testing device, where assumptions lead to observable collective behavior via simulation. Of course, the more complicated a model, the more parameters it contains, the greater is the risk of overfitting and nonidentifiability. Another example is the tumor microenvironment modeled by Kather et al., which can be calibrated in a patient‐specific way. In the models we have introduced, the ABM requires more input than the ODE: the initial distribution of pERK values, and the death and birth thresholds. Analytical methods are required when parameters do not represent directly observable quantities. In that case, statistical estimation techniques, such as maximum likelihood estimation or the method of moments, are applied to a given dataset to select appropriate values. If analytical methods are not suitable for a given ABM, methods involving the generation of simulated data need to be considered. Those methods can be classified in frequentist approaches and Bayesian approaches. Frequentist approaches are distance‐based or likelihood‐based (e.g., the simulated minimum distance method or the methods of simulated moments). One example for calibration of ABMs using machine learning would be the work from Lamperti et al. Cockrell and An simulate systemic inflammatory response syndrome, with an ABM used to compare treatment options with different dosing regimens and drug combinations.

ABMs can be computationally challenging, depending on their complexity and the chosen model paradigm. In cellular Potts models, changes in cell shape and direct cell‐cell interaction are directed by Monte Carlo simulations and energy minimalization. Off‐lattice methods bring the benefit of a more realistic simulation because cells can have various positions with respect to each other and freedom to move in any direction instead of being ordered on a grid. However, this comes with the disadvantage of higher computational cost, because special algorithms are necessary to efficiently handle cell‐cell neighborhoods. During movement or division of cells, cell collisions with nearby cells need to be considered, which can be challenging in densely packed areas or populations. Improved approximations of cell biomechanics invariably come with computational cost. , ODEs, in contrast, seldom require large computational resources.

A strength of ABMs is the use of biological rules which makes communication of the model easier and more intuitive for a non‐modeling audience. Components of two ABMs can be combined in a modular fashion to create meta‐models. Each ABM can have different agents and outcomes, making comparison between different models challenging. On the other hand, advantages of ODE based models commonly used in the PKPD community are the simplicity of implementation, and relative ease in fitting experimental data in a statistically robust manner with clear decision rules for model selection and well‐established tools for model qualification. Communication with non‐modeler audiences relies on understanding of mass transfer and binding kinetics.

Multiscale modeling aims to include various spatio‐temporal scales from atomic to molecular, cellular, multicellular, organ, and whole body. When explicit representation of individuals is not needed, a continuous description with differential equations can be used. On the molecular scale, interactions (e.g., receptor‐ligand interactions, consumption and production of oxygen, nutrient, and cell‐cell signaling molecule concentration) can be described with ODEs. To model local conditions and environmental changes, such as availability of oxygen, nutrient, and hormones controlled by diffusion from molecularly rich regions (e.g., blood vessels and tumor edge), PDEs may be used. ODEs or PDEs may be integrated into simulations to lower computational cost. This type of hybridization of a discrete model in a continuum environment is often found in the most complete descriptions of the tumor morphology.

Comparing hybrid multiscale ABMs to a multiscale ODE model, as in the work of Milberg et al., we can see additional differences between those methods. Their physiology‐based quantitative pharmacology model predicts how the interaction of the immune system and the tumor microenvironment in a patient affects checkpoint blockade therapies administered as mono‐, combo‐, and sequential therapies. Agents may be located in compartments: lymph node, blood, tumor, lungs, gastrointestinal tract, spleen, and liver, and the periphery. Heterogeneity is represented in the percent expression of each immune checkpoint in the cancer cells as an input into the model. The model consists of 282 ODEs and 218 algebraic equations to describe interactions between cells and trafficking of cells. In contrast, an ABM‐based implementation of this model would require only six agents (with attributes in parentheses): Antigen, antigen presenting cells (resident, mature, and aCTLA4 expression), CD8+ T‐cells (naïve, primed, activated, CTLA4, and PD‐1), Tregs (CTLA4 and PD1), myeloid derived suppressive cells, and tumor cells (proliferating and PD‐L1). The behavior of these agents is then controlled by a total of 16 interaction rules, with each cell type making use of a subset of these. In principle, this could be much more efficiently implemented and easily adapted should further interactions or cell types be required.

There are various hybrid multiscale ABMs in the literature. Oduola and Li modeled cancer growth during treatment with lapatinib using a multiscale model with a stochastic hybrid system, where concentration of proteins and gene expression levels were represented with ODEs. The cellular level contains a cellular automata model on a grid. A multiscale ABM of tumor angiogenesis is provided by Olsen et al. This model includes the molecular level (VEGF and diffusion), cellular level (genetic control and space), and tissue level (cells, blood vessels, and angiogenesis). Cancer cells are modeled as agents on a three‐dimensional grid. Chaplain and Powathil developed two hybrid multiscale models that study the effects of cancer treatments, including a combination of radiation and chemotherapy. , The microenvironment contains the concentration of oxygen modeled with a PDE. At the subcellular level, the cell cycle depends on concentration of complexes that are described with ODEs. Cellular automata and Potts models simulate the cellular level. Cess and Finley created a multiscale ABM of macrophages and T cells within a tumor. Tumor cells, M1 and M2 macrophages, and T cells are modeled as agents on a two‐dimensional lattice, whereas diffusible mediators, such as IL‐4 and IFN‐γ, are simulated with PDEs. Neural networks are used to reduce the mechanistic model into a simple input/output model. The multiscale compartment model by Gong et al. describes the biological processes involved in tumor development and anti‐tumor immune response. Cytotoxic T lymphocytes and cancer cells are modeled as agents in a three‐dimensional space, with division, migration, cytotoxic killing, and immune evasion. PDEs describe the molecular scale (IL‐2 secretion and transport).

ABMs are computational models with heterogeneous agents in which the behavior of individual agents is fundamental. Simulations may be intuitive because they recapitulate biological processes. ABMs are suitable for simulating complex biological systems with subscale components (molecular, cellular, tissue, and organism) and inherent emerging behavior. On the other hand, systems of ODEs are well‐suited for simulating processes that can be approximated as homogeneous, well‐mixed systems, and would be best suited for traditional pharmacometric analyses with sufficient data (population PK and PD models, and physiology‐based PK models), or for simplistic theoretical PKPD models. For quantitative clinical pharmacology models, ODEs can also be implemented to recapitulate complex biological systems but would rely on extensive model assumptions, including parameter distributions. In summary, ABMs can provide more detailed insights into complex biological systems and are often complemented with ODEs in hybrid multiscale models. Both methodologies have their strengths and weaknesses, depending on context and purpose. With the advent of more single‐cell experiments, gene expression, and spatial transcriptomics, and other technological advances in imaging, we anticipate that the use of ABMs in discovery and drug development will increase with direct applications in PKPD and pharmacometrics to help elucidate dose scheduling and rationalize combination strategy in oncology and other therapeutic areas.

CONFLICT OF INTEREST

All authors declared no competing interests for this work.

Supplementary Material

Click here for additional data file.