PhysiPKPD: A pharmacokinetics and pharmacodynamics module for PhysiCell

Pharmacokinetics and pharmacodynamics (PKPD) are key considerations in any study of molecular therapies. It is thus imperative to factor their effects into any in silico model of biological tissue involving such therapies. Furthermore, creating a standardized and flexible framework will benefit the community by increasing access to such modules and enhancing their communicability. PhysiCell is an open-source physics-based cell simulator, i.e., a platform for modeling biological tissue, that is quickly being adopted and utilized by the mathematical biology community. We present here PhysiPKPD, an open-source PhysiCell-based package that allows users to include PKPD in PhysiCell models. Availability & Implementation The source code for PhysiPKPD is located here: https://github.com/drbergman/PhysiPKPD.

the effect of a given drug are commonly broken down into pharmacokinetics, or PK, and pharmacodynamics, or PD. Together, they are often called PKPD. Pharmacokinetics describes how a drug is transported throughout the body to get to the cells it acts on, i.e., the drug exposure. Pharmacodynamics describes how the drug then acts on those cells, i.e., the drug response.
The dynamics involved in PKPD are generalizable to many different drugs acting in different ways. This makes it an appealing target for a single module that can flexibly handle these processes on a platform such as PhysiCell. Features on both the PK side -such as dosing schedules, loading doses, elimination, and distribution rates -and the PD sidesuch as mechanisms of action (MOA), effect, the half maximal effective concentration or EC50, and hill coefficients -are amenable to this level of abstraction.
We present here PhysiPKPD, a standardized framework to incorporate these PKPD processes in PhysiCell. We also provide several examples demonstrating PhysiPKPD and providing users with two template projects to aid in the model-building process. In this first version, we provide three options for modeling pharmacokinetics: (1) one-compartment models, (2) two-compartment models, and (3) the ability to supply a Systems Biology Markup Language (SBML)-defined model. For options (1) and (2), the equations must be linear. On the PD side, substrates cause damage to cells based on the internalized substrate.
As this damage accumulates within a cell, the MOA-associated rate parameter(s) tend towards a user-defined saturation rate. While we expect most users to add both PK and PD for a given substrate, they can be added independently if desired. For substrates that do have PD, these are set for each cell type in the model.

IMPLEMENTATION
We provide two ways to create and run PhysiPKPD models, shown in Figure 1. First, the sample projects that come with PhysiPKPD demonstrate the MOAs implemented in PhysiPKPD ( Figure 1A). Second, the user can use the template projects to jumpstart the model-building process with the full range of possible parameters and code present throughout the PhysiCell repository ( Figure 1B). We now explain how PhysiPKPD achieves these dynamics. In what follows, S will stand for the name of a substrate, C for the name of a cell type, and X for an MOA.

Pharmacokinetics
For pharmacokinetics, the user must first specify which substrates in the model follow a PK model. This is done by adding these substrate names to the comma-separated string user parameter PKPD_pk_substrate_names.

PK models
For each of these substrates, the user should specify a PK model (see Table 1). If this is not specified, PhysiPKPD will attempt to use a two-compartment model. In addition, the user should specify a Biot number, which PhysiPKPD uses to set the ratio of perivascular substrate concentration to that in the blood vessel [13]. PhysiPKPD will default to a value of 1, indicating equal concentrations inside and outside the capillary walls.  Table 1. PK model specifications. The specification column indicates the string value to be used in the configuration file for the parameter S_pk_model.

Model
Description Specification One-compartment Circulation compartment with linear elimination 1C Two-compartment Circulation and peripheral compartments with linear clearance rates 2C

SBML-defined
Any SBML-defined model with one species named circulation_concentration SBML For example, add the following to User Parameters in the configuration file: <S_pk_model type="string">2C</S_pk_model>.
One-compartment models. The PhysiPKPD one-compartment model obeys the differential where C represents the circulation concentration of the substrate and is the elimination rate (see Table 2). PhysiPKPD uses the value of C to update the Dirichlet conditions in the PhysiCell microenvironment after multiplying by the Biot number. Two-compartment models. The PhysiPKPD two-compartment model obeys the differential Here, C and are as above. P is the periphery concentration. The parameters k 12 and k 21 are the intercompartmental clearance rates and R is the ratio of the volumes of the central and peripheral compartments (see Table 2). PhysiPKPD uses the value of C to update the Dirichlet conditions in the PhysiCell microenvironment after multiplying by the Biot number. The inclusion of the periphery compartment allows for biphasic elimination in the central compartment.
SBML-defined models. If the above two models are inadequate for a user's purposes, an SBML file can be used to specify a PK model. We have used Copasi (RRID:SCR_014260) [14] to build such models, but any program that outputs an SBML file will work, e.g., sbmlutils [15].
PhysiPKPD will use the first state variable of the model to update the Dirichlet conditions in the PhysiCell microenvironment after multiplying by the Biot number.

Dosing schedules
For the one-and two-compartment models, all parameters and dosing events must be specified within the configuration file (see Table 3). For any missing parameters, PhysiPKPD will issue warnings and use default values where it can, and it will throw errors where it must. All substrates are assumed to be given intravenously so that the concentration in the central compartment has a one-time increase upon dosing, S_central_increase_on_dose.
These doses are given at regular intervals, S_dose_interval, until either the simulation ends or until the maximum number of doses has been administered. The first dose can be given at a fixed time or based on the confluence in 2D in the entire rectangular microenvironment. A loading dose can also be set with a fixed number of doses. Future releases can include other methods of administration, e.g., oral, and finer-grained control, such as more sophisticated timings for doses, e.g., M-F dosing.
For SBML-defined models, the user must include dosing events in the SBML file itself. In Copasi, this can be done by creating Events that increase the concentration of a compartment(s) at certain times. In the future, we hope to allow the user to specify a CSV file with the dosing times and amounts for a substrate along with an SBML-defined system

Damage accumulation
For each cell type affected by a particular substrate, a damage-accumulation model should be specified. If this is not specified, PhysiPKPD will default to the concentration-based model. The two options are the concentration-based model and the amount-based model (see Table 4). Both of these models use proportional to the internalized substrate. This damage is repaired at the linear rate r 1 and the constant rate r 0 . As damage is an abstract quantity, the proportionality constant for the internalized substrate causing damage is set to 1. All three of these parameters (Table 5) must be included in the custom data for C, the cell type affected by substrate S. The differential equations are then given by  Set all in Custom Data of the cell definition of C. Table 6. Cell effect parameters.

S_moa_is_X
Identifies whether the MOA of S on C is X None (Values > 0.5 trigger the MOA)

S_X_saturation_rate
Limiting rate of X as damage from S grows towards infinity min −1

S_X_EC50
Damage from S at which the rate of X is halfway between the base and saturation rates damage

Hill power None
Replace X with one of the following: prolif, apop, necrosis, or motility.

Cell effect parameters
To apply the desired MOA of S on C, the damage variable S_damage is used as the input to a Hill-type function. The user must specify three parameters for this Hill-type function in addition to identifying this MOA (see Table 6). The four MOAs currently implemented in PhysiPKPD are proliferation (prolif), apoptosis (apop), necrosis (necrosis), and motility (motility). Replace X in Table 6 with the parenthetical name of the desired MOA. By default, all MOAs are assumed off, so including S_moa_is_X in the custom data is only necessary if set to true. Note that custom data in PhysiCell must be of type double, so "true" means a value > 0.5.
Multiply-targeted MOAs. In the case that two or more substrates have the same MOA on a given cell type, we compute the drug effects as multipliers or factors. That is, the algorithm computes the saturation factor, f sat , based on the user-supplied saturation rate and uses the damage to compute a factor, f, between 1 (no change) and the saturation factor (max change): When for example, multiple substrates affect proliferation, the factors for each substrate are multiplied to give the final factor for the cell. In the case of necrosis, the base rate for necrosis is often set to 0, which causes problems in computing f sat . Thus, we instead add the effects of multiple substrates targeting necrosis.

Example results
We provide several examples of implementations of PhysiPKPD: one for each MOA, one combination therapy with one anti-proliferative drug and one pro-apoptotic drug ( Figure 1A), and one using a confluence condition to start dosing. Follow the README.md file instructions to get PhysiPKPD set up in the root PhysiCell directory. The output will go to the output folder. There are many parameters that can be changed to explore various behaviors, even restricted to just one MOA. See Tables 2, 3

DISCUSSION
A PKPD model that is ready to use "out of the box" with PhysiCell will greatly benefit those seeking to use PhysiCell as a platform for agent-based modeling. It will also create a standard that can be used by any modeler looking to incorporate PKPD into their PhysiCell models.
PhysiPKPD provides exactly this functionality. We designed PhysiPKPD so that the user can largely forego editing the C++ code, instead selecting PK/PD models and setting parameter values in the PhysiCell configuration file. By integrating with the native substrate class in PhysiCell, we can provide all these dynamics without introducing additional dependencies. When the PK models are the standard one-or two-compartment models, we use the analytic solution to provide fast and accurate updates to the pharmacokinetics. More importantly, the PD model also yields an analytic solution that we use to efficiently solve these dynamics, particularly compared to using numerical methods.
The flexibility of PhysiPKPD allows the user to readily add as many drugs as desired and have them affect the various cell types in unique ways. Each drug can follow its own PK model, and if the one-and two-compartment models are not sufficient, an SBML-defined PK model can take their place. Each drug can affect each cell type in a unique way, combining any number of PhysiPKPD's MOAs. By using the damage_coloring function of PhysiPKPD, a user can easily observe the extent to which up to two drugs are affecting a given cell type.
This helps in narrowing down a parameter regime that results in observed and/or desired drug exposure and response.
We have also aimed to make PhysiPKPD easy for others to pick up and use. The samples provided are a great place to start learning these capabilities of PhysiPKPD. For those looking to build a new project with PKPD included from the beginning, we provide two sample projects to help the user get started: one with and one without SBML support for a PK model. Finally, adding PK and/or PD to a substrate in a current PhysiCell project is simply a matter of adding the correct header files, updating the Makefile, and adding the parameters to the configuration file.
There are many improvements and additions we hope to make to PhysiPKPD in the near future. With regards to PK, we hope to make it easier to include dosing events in an SBML file by allowing the user to supply a PK model and a CSV with the dosing information. With regards to PD, there are many features that can be added to allow for greater flexibility in terms of mechanisms of action, effect models, and integration with intracellular signaling.
We have so far assumed that the AUC of a substrate inside a given cell, what we called damage, can be used to determine the effect of a drug on that cell. This is similar to irreversible effects models that use the AUC of the drug concentration to determine the response [16]. Many other effects models have been used, including simple direct effect, indirect response, and signal transduction models [17]. Recent work has used Bayesian inference to determine the distribution of delay times in response to therapeutic agents [18], making stochastic effects models an appealing next step as well. This could be added directly into previously-studied intracellular models or be used as a means to coarse-grain such complex models while still allowing for heterogeneity within cell types.
We are also working with the makers of PhysiCell to further simplify and standardize the integration of PhysiPKPD into PhysiCell. A key goal of this partnership will be to allow the user to include PKPD while only needing to specify the relevant parameters and without directly interacting with the C++ code. We look forward to continuing to develop this tool as the community uses it and seeks new features.

DATA AVAILABILITY
Snapshots of the code and underlying data is available from the GigaScience GigaDB repository [19].