In vitro PK/PD modeling of tyrosine kinase inhibitors in non‐small cell lung cancer cell lines

Abstract Tyrosine kinase inhibitors (TKIs) are routinely prescribed for the treatment of non‐small cell lung cancer (NSCLC). As with all medications, patients can experience adverse events due to TKIs. Unfortunately, the relationship between many TKIs and the occurrence of certain adverse events remains unclear. There are limited in vivo studies which focus on TKIs and their effects on different regulation pathways. Many in vitro studies, however, that investigate the effects of TKIs observe additional changes, such as changes in gene activations or protein expressions. These studies could potentially help to gain greater understanding of the mechanisms for TKI induced adverse events. However, in order to utilize these pathways in a pharmacokinetic/pharmacodynamic (PK/PD) framework, an in vitro PK/PD model needs to be developed, in order to characterize the effects of TKIs in NSCLC cell lines. Through the use of ordinary differential equations, cell viability data and nonlinear mixed effects modeling, an in vitro TKI PK/PD model was developed with estimated PK and PD parameter values for the TKIs alectinib, crizotinib, erlotinib, and gefitinib. The relative standard errors for the population parameters are all less than 25%. The inclusion of random effects enabled the model to predict individual parameter values which provided a closer fit to the observed response. It is hoped that this model can be extended to include in vitro data of certain pathways that may potentially be linked with adverse events and provide a better understanding of TKI‐induced adverse events.


INTRODUCTION
Targeted therapy in non-small cell lung cancer (NSCLC), has provided many patients with NSCLC with improved prognosis and survivability. 1,2Tyrosine kinase inhibitors (TKIs) are one of the main classes of drugs that are used as a targeted therapy for NSCLC treatment. 24][5] Common examples of TKIs include erlotinib and gefitinib.
As with all medications, there is a chance that patients may experience adverse events with a TKI therapy.Serious adverse events may become life-threatening, as the adverse event begins to interfere with their physiological system. 6Suspension of NSCLC treatment may also occur with serious adverse event diagnosis. 7This can provide an opportunity for the tumor to grow and potentially spread to other regions.Unfortunately, patients may also die due to a serious adverse event. 6or many adverse events, the mechanisms are unclear. 8,9][12] There are limited in vivo studies that investigate the effects of TKIs on different protein expressions or the activation of certain genes.Most of the in vivo studies focus on tumor response to the TKIs.4][15] Although these studies are useful, it is often difficult to interpret how these factors could potentially be related to the adverse event.
7][18][19] These studies could provide potential insights into the mechanisms for different adverse events.
Analysis on the phosphorylation pathways of RTKs can enable a model to be developed to simulate the activation or deactivation of certain growth factors which have been shown to be related to pulmonary adverse events. 20he developed in vitro model can then be transformed by considering a tumor compartment to generate an in vivo model to simulate the possible mechanism for pulmonary adverse events due to TKIs.In order to utilize these studies, it is imperative to capture the transition from when the drug is introduced to the medium and all the way to any effects the drug could have on the potential targets.The pharmacokinetic (PK) relationship among TKIs, NSCLC cells, and the cell medium remains limited.
Nonlinear mixed effects (NLME) modeling enables PK and pharmacodynamic (PD) parameter values to be estimated from an observed population. 21The use of in vitro data and NLME modeling will allow an in vitro PK/PD model to be developed for TKIs in NSCLC cell lines.This model can then potentially be expanded to simulate potential pathways for adverse events.
The aim of this paper is to develop an in vitro PK/PD model for TKIs in NSCLC cell lines with the use of in vitro data and NLME modeling.The Methods section will cover the data collection, model development (including the structural identifiability and sensitivity analysis), and the incorporation of NLME into the structural model.In the Results section, the sensitivity analysis plot is displayed along with the parameter estimates and model simulations.These results are then summarized and discussed in the Discussion section.

In vitro data collection
In vitro data were extracted from cell viability studies that observed the changes in tumor cells' viability in response are many in vitro studies which do observe the changes in different gene expressions due to TKIs.

WHAT QUESTION DID THIS STUDY ADDRESS?
The development of an in vitro TKI model based on publicly available data and the application of nonlinear mixed effects modeling.

WHAT DOES THIS STUDY ADD TO OUR KNOWLEDGE?
With the aid of cell viability data and nonlinear mixed effects modeling, parameter estimates were obtained for an in vitro pharmacokinetic/pharmacodynamic model for the TKIs alectinib, crizotinib, erlotinib, and gefitinib.

HOW MIGHT THIS CHANGE CLINICAL PHARMACOLOGY OR TRANSLATIONAL SCIENCE?
This model can be further extended to model other potential pathways for adverse events using key observations made in in vitro studies.
to TKIs over time.Studies were included if the following details were present in the study: TKI concentrations versus time, duration of the treatment, number of cells that were placed in each well, and the name of NSCLC cell lines that were used.TKI groups that had less than 10 viability analyses were removed and the treatment duration was also limited to 72 h only.The decisions for these additional criteria to be included were made to provide adequate populations for the NLME modeling and 72 h was the most common treatment duration.As a result, cell viability for the following TKIs were collected; alectinib (120 cell-viability observations), crizotinib (129 cell-viability observations), erlotinib (61 cell-viability observations), and gefitinib (82 cell-viability observations).WebPlotDigitizer was used to extrapolate the tumor cell viability response, the minimum concentration to produce 50% inhibition (IC 50 ), the maximum inhibition response (I MAX ), and in cases where there were extra TKI dosages, the TKI concentration. 22lectron micrographs of each of the cell lines that were included were used to estimate the initial cell volume which was based on sphere volume calculations.The image software tool ImageJ was used to estimate the cell volumes. 23The in vitro tumor cell viability data as well as the cell volumes that were estimated, along with their reference sources, can be found in Tables S1-S5 in the Supplementary Materials.

Structural model in vitro TKI PK/ PD model
In an in vitro setting, the main state variables for the PK process are the concentrations of the drug in the medium and the cell.Figure 1 provides a simplified schematic diagram of the PK process in an in vitro setting.Although the TKIs will produce their effects in cells, the TKIs are first introduced into the medium, thus the absorption process from the medium into the cell has to be captured.In Figure 1, the only parameters that are known are the initial concentration of the TKI in the medium and the initial volume of the cell.Note that, in this case, the volume of the cell refers to all the cells that are grouped together in a well.
The following ordinary differential equations (ODEs) characterize the changes in concentration of TKI in the medium and in the cell: The definitions for the parameters for all of the equations as well as their respective units can be found in Table 1.Note that the time "(t)" has been dropped for all the state variables for the sake of brevity.
One of the main effects of TKIs is the reduction of the tumor size, which would reduce the NSCLC cell volume.This reduction in volume would affect the concentration of the TKI in the cell.Therefore, the volume of the cell at time t, V C(t) , can be described as: The PDs of TKIs in NSCLC cell lines can also be represented by a first order ODE. 24A first order indirect PD model was chosen to take into account the impact of TKI concentration on the cell viability and also to generate an appropriate dose response curve.The main effect of TKIs on NSCLC cell lines that is observed in cell viability studies is a reduction in cell viability which reduces the cell volume, due to inhibition of phosphorylated RTKs. 25 Equation 4 was derived to characterize the change in tumor (1) Pharmacokinetic process in an in vitro setting.TKI C and TKI M refers to the concentration of TKI in the cell and medium respectively.K DRUGIN and K DRUGOUT refers to the absorption rate of TKI form the medium into the cell and from the cell back into the medium respectively.V C and V M refer to the volume of the cell and medium, respectively.TKI, tyrosine kinase inhibitor.
cell viability as a result of TKI administration, through the use of an indirect effect compartment model. 26te that for this system of ODEs, only the parameters V C(0) , I MAX , and IC 50 are known.Moreover, TumorCellViability is the only state that is observed.Although the volume of the medium is not known, the volume ratio, V R , between the cell volume and cell medium, V M , is known (between 1 and 8). 27ence, V M can be represented as:

Structural identifiability analysis
A structural identifiability analysis assesses whether the unknown parameters of a model can be identified uniquely or otherwise based on the structural model and the states of the model that are observed.Consider the following as a generalized ODE model formulation: where dx i (t)∕ dt denotes the rate of change of the states of the ODE.The function f (x(t), , u(t)) encompasses the ODE dynamics for the vector with elements (x(t)), (which represents a vector of unobserved and unknown parameters), and u(t), the input vector.The y(t) represents the output and is considered as a smooth function with the same elements as the system states.Last, x(0) represents the vector of initial conditions for the states (x INT ).
A specific parameter, i , is said to be globally identifiable if for all of t, the following condition is met: In other words, if i can be uniquely determined from y(t), then i is globally identifiable.Therefore, i can still be classed as locally identifiable if for i there exists a neighborhood, V ( * ) where Equation 7 holds.][30][31] For a model to be classed as at-least locally identifiable, none of the parameters can be deemed as unidentifiable.Subsequent numerical estimations of unidentifiable parameters will not be reliable as there are infinitely many values that an unidentifiable parameter can take for the given observation.
In this study, the input-output relation and Taylor series approaches were used to assess the structural identifiability of the system.The Taylor series method utilizes the Taylor series coefficients of the observed output, to assess for structural identifiability. 30,31This study utilizes the Lie derivatives of the observed output to generate the inputoutput relation.The coefficients of the monomials that were generated from the input-output relation approach are then used to assess for the structural identifiability of the model.These coefficients (Taylor series and inputoutput relations) include some (or in some cases all) of the parameters that are in the ODE system.These coefficients are also assumed to be equal to unique solutions, which can be used to identify whether a parameter is identifiable.
T A B L E 1 Definitions for all of the parameters and state variables.
Both methods were used to overcome the issues with computational memory limitations and model complexity. 32ased on the model ODEs and knowledge that the parameters V C(0) , I MAX , and IC 50 are known, the model is deemed as locally identifiable.A detailed analysis of the structural identifiability of this model in an NLME context can be found in the supplementary materials.The Taylor series approach was performed in MATHEMATICA and the input-output approach was performed in MAPLE. 33,34

Sensitivity analysis
Sensitivity analysis effectively analyses the impact that changes in the model parameters will have on the model output. 35This is useful for validating whether the model is truly characterizing the system that the model aims to represent.Moreover, for models that have many parameters, it can aid in the identification of parameters which could potentially be removed from the model, reducing its parameter complexity which may additionally improve the structural identifiability of the model.Sensitivity analysis was performed on the model presented using the Sobol sensitivity indices method. 36This global sensitivity analysis method is based on analyzing the variance of the output state based on different parameter values.The variance is then decomposed into different contribution scores for each of the parameters and their impact on the output's variance.Any parameter which had a Sobol value of 0 would indicate that this particular parameter has no impact to the output of a model.[39]

Nonlinear mixed effects modeling
The NLME modeling comprises of both the population effects and the random effects.The population effects refer to the population parameter estimates, whereas, in this study, the random effects refer to the between subject variability (or interindividual variability) as well as the residual unknown variability (RUV). 21The general NLME predicted output form is given by: where i represents the individuals in the population and j represents the number of observations.The x ij refers to the states in the model ODEs that describe the state for individual i and observation j.The i refers to the individual estimated parameter values.The g x ij , i , × ij represents the residual unknown variability.In this study, MONOLIX was used to perform the NLME modeling. 40ONOLIX utilizes the stochastic approximation expectation maximization algorithm in order to estimate population parameter values. 41In short, individual parameter values are drawn from a conditional probability which contains population parameter estimates that are set by the user.These individual parameter values are then assessed for the likelihood of being able to produce the observed output.New population estimates are then formed based on these individual parameter values and the process repeats until the estimated parameter values begin to converge.

Parameter estimation
For each of the analyses, a constant error model was chosen as the best model to characterize the RUV.Equation 8can therefore be represented as: where a is a constant value (see Table 2) and ij often represents a normal distribution with a mean 0 and variance of 1.Other error models were assessed, such as the proportional and combined error model.However, the constant error model obtained the lowest standard errors and relative standard errors (RSEs).
From Equation 5, V R can be used to estimate V M .All of the estimates for V R lie in the 1-8 range (Table 2).For all of the TKIs, the estimated values for K DRUGIN are higher than the estimated value for K DRUGOUT .The model for the alectinib group has the lowest K OUT value (0.021), whereas that for crizotinib has the highest parameter estimate value for V R , whereas the erlotinib estimates yield the highest K IN value, as well as lower (8)

Prediction versus observation
Figures 3 and 4 showcase the model predictions versus the observed data using the population estimates only (Figure 3) and the inclusion on random effects (Figure 4).
In Figure 3c, the model parameter estimates for the erlotinib group only using the population estimates only has the widest range for the predicted values between 0 and 90%.The other TKI groups in Figure 3 show a restriction between 0 and 70% or, in, the case of the gefitinib, only 0 and 40% (Figure 3d).For all of the plots, the observed tumor cell viability can reach above 100%.
Figure 4a, shows that the alectinib individual parameter estimates were still constrained to lie within 40 to 90% of tumor cell viability.The other TKI groups (Figure 4bd), have predicted cell tumor viability ranges that were the same as for the observed tumor cell viability ranges.Moreover, the crizotinib, erlotinib, and gefitinib groups each achieved a correlation value of 1.The specific cell lines were chosen as these cell lines had the most varied TKI concentrations and the most data associated with them.Similar to the previous plots, the alectinib simulation response lies between 50% and 80% of tumor cell viability (Figure 5a).The crizotinib observed response is the only plot which shows a gradual decrease in cell growth due to the increase in crizotinib concentration (Figure 5b).Both the erlotinib and gefitinib plots show that the average observed values for the NSCLC cell lines HCC827 and PC9 (respectively) can fluctuate as the concentration increases.The simulated dose response curve using only the population parameter estimates only can be found in Figure S1 in the Supplementary Materials.The sensitivity analysis for the parameter TKI M (Figure 2a) shows that parameters K DRUGIN and K DRUGOUT are the most influential parameters.In Equation 1, both the volumes enter into the denominators for the latter part of the ODE, which implies that, unless the volumes are extremely small, the numerator (including the parameter K DRUGOUT ) will have more of an impact.As for the parameter K DRUGIN , it is the only parameter that is multiplied by TKI M in Equation 1. K OUT , K IN , and I MAX are also seen as important parameters for TKI M (Figure 2a).These parameters are also influential for TKI C and tumor cell viability (Figure 2b,c).The state TKI C also influences TKI M .The concentration of TKI in the cell is also greatly influenced by the volume of the cell which is affected by the tumor cell viability, where the parameters K OUT , K IN , and I MAX have significant influence.The IC 50 parameter, however, does not seem to influence any of the states to a notable effect.This may be due to the fact that, based on Equation 4, the I MAX value always impacts on the tumor cell viability more than IC 50 .The parameters estimated are not dependent on the NSCLC cell line mutation status.This decision was made to reflect the fact that the IC 50 and I MAX values would reflect any mutation that was present.Although many PK/PD models often fix the I MAX value at 1, 42 the T A B L E 2 Parameter estimates for the in vitro PK/PD models.to include the I MAX values was to aid the model to enhance characterization of the data that were observed.

Parameters
All of the models predicted higher K IN rates compared to the K OUT rates, which are in line with what is known about tumors in their actively proliferating state.In the majority of cases, the K DRUGOUT and V R estimates had the highest RSE values.With the exception of the erlotinib parameter estimates, the shrinkage for the estimates for K DRUGIN , K DRUGOUT , and V R are above 99%, which suggests that the individual parameter estimates will likely be similar to the population parameter estimates.The shrinkage values for the K OUT estimates were lower for the crizotinib, erlotinib, and gefitinib estimates when compared to the other parameter shrinkage values.This indicates that there is less uncertainty with respect to the individual parameter estimates for K OUT .Given that K OUT has significant influence on the tumor cell viability (Figure 2c), the individual estimates for K OUT have an impact on the prediction versus observation response when individual estimates are included.
For the predictions versus observations plots (Figures 3  and 4), the alectinib plots have restricted ranges which could possibly be due to the restrictions placed on the I MAX values, as discussed previously.In Figure 4, the crizotinib, erlotinib, and gefitinib plots all achieved a correlation value of 1 suggesting that the inclusion of random effects enables the model to capture the observed data efficiently.The shrinkage values for the alectinib parameter estimates were above 90% for all of the parameter estimates except for K IN , which suggests that the individual estimates were similar to the population parameter estimates.This could explain why, even with the inclusion of individual estimates, the prediction versus observation plot for alectinib (Figure 4a), remained similar to that for the population estimates only (Figure 3a).Moreover, this suggests that perhaps there are not enough data for the model to be able to fully predict the PK and PD processes for alectinib.Nonetheless, all of the plots had high positive correlation values of at least 0.9 which suggests that the model's assumed relationship between concentration and tumor viability inhibition is similar to the observed data.
In Figure 5, the average simulated dose response curves for most of the plots show that the inclusion of the random effects provided a closer fit for the model when compared to the data.In the alectinib plot (Figure 5a), the simulated response remains at 50% of tumor cell viability whereas the observed response actually decreases to less than 20%.This outcome may be due to the low K OUT value compared to other TKI estimates (Table 2).Although the IN VITRO MODELING TKIS NSCLC CELL LINES model also has a similar predicted K OUT value (0.026), the variance is a lot higher than compared to the alectinib variance for K OUT , which could suggest that other individual K OUT values were higher, enabling a more accurate prediction.For the alectinib variance for K OUT (0.082), this leads to the assumption that individual K OUT estimates will most likely remain close to their population value.This could impede the ability of the predictions to be lower for the higher concentrations.Moreover, even though the I MAX values for the H3122 NSCLC cell line were on average 0.9, Figure 2 shows that the most influential parameter for tumor cell viability response is actually K OUT .For Figure 5c,d This in vitro model shares many similarities with its in vivo model counterpart.4][45][46] The tumor cell volume is also considered, which for many models also incorporates this as part of the effect compartment. 47he next step for this model is to include the phosphorylation cascade in order to incorporate target growth factors that have been linked to pulmonary adverse events, such as interstitial lung disease (ILD).Common approaches for predicting ILD involve the use of artificial intelligence (AI). 48,49However, although these approaches can achieve high prediction scores, there is limited information on why certain risk factors may lead to an increased likelihood of ILD.This unfortunately further confounds understanding of the mechanism behind ILD. 50The extended in vitro PK/ PD model, includes the potential mechanistic pathway of ILD (based on the phosphorylation data from pulmonary cells).When the PD elements of the in vitro model are incorporated into a patient model, the different TKI medications as well as dosing regimens can be implemented and the likelihood of ILD onset will be mechanistically based on the TKI and dosing regimen for each patient.
The main limitations in this analysis to note are that the observations and predictions are based on data that were extracted digitally and thus are only an estimate and that the data obtained are based on different in vitro studies.cases where the NSCLC lines were used, there may still be differences in the experimental and growth conditions and proliferation rates.These differences will also impact on the I MAX and IC 50 values as well as the model parameter estimates.The estimates obtained through NLME modeling are dependent on the data used.Due to these confounding factors and the sparsity of the available data, the inclusion of more data would be necessary for the alectinib, crizotinib, and gefitinib models to be utilized in a robust predictive capacity.However, the erlotinib model demonstrates that it is feasible to capture the PK/PD process based on the population estimates alone (Figure 3c).This specific model could potentially be further explored to incorporate the RTK pathway and simulate the possible mechanism(s) for adverse events.
In addition to this, the in vitro model does not include other factors, such as protein binding, the number of receptors available, active TKI metabolites, the diffusion rate from the cell membrane to the cell nucleus, etc.The main reason for these exclusions was to reduce the complexity of the model and for the model to be structurally identifiable.To conclude this study, an in vitro TKI PK/PD model was developed with the use of NLME modeling.It is believed that this is the first instance where estimates for these parameters for these TKIs in an in vitro setting have been obtained without prior experimentation for data collection.It is hoped that this model can be furthered developed in order to model other potential pathways for adverse events.
M (0) TKI concentration in the cell medium Dose (μmol/L) TKI C (0) TKI concentration in the cell 0 (μmol/L) TumorCellViability (0) Tumor cell viability 100 (%) K DRUGIN (h −1 ) TKI absorption rate from medium into cell K DRUGOUT (h −1 ) TKI absorption rate from cell back into the medium V M (L) Volume of the medium V C(0) (L) Initial cell volume K IN (h −1 ) Tumor cell proliferation rate I MAX (%) Maximum tumor cell viability inhibition IC 50 (μmol/L) TKI concentration required to produce 50% tumor cell viability inhibition

Figure 2
Figure 2 visualizes the Sobol indices for each of the parameters in the model.In Figure 2a, K DRUGIN has the highest Sobol indices at 0.8.Other parameters that have Sobol indices greater than 0 include K OUT , K IN , and I MAX .For TKI C (Figure 2b), V C , K IN , I MAX , K OUT , and K DRUGOUT have higher Sobol indices compared to the other parameters.Last, for the observed output state (tumor cell viability), K OUT has the highest Sobol indices of ~1.Other parameters with high Sobol indices are K IN and I MAX (Figure 2c).
values for the random effects compared to the other TKIs.The gefitinib estimates have the lowest population estimate for the parameter K IN .Overall, the erlotinib parameter estimates have the lowest shrinkage compared to those for the other TKIs.

Figure 5
Figure 5 visualizes the simulated dose response curves based on the parameter estimates.The specific cell lines were chosen as these cell lines had the most varied TKI concentrations and the most data associated with them.Similar to the previous plots, the alectinib simulation response lies between 50% and 80% of tumor cell viability (Figure5a).The crizotinib observed response is the only plot which shows a gradual decrease in cell growth due to the increase in crizotinib concentration (Figure5b).Both the erlotinib and gefitinib plots show that the average observed values for the NSCLC cell lines HCC827 and PC9 (respectively) can fluctuate as the concentration increases.The simulated dose response curve using only the population parameter estimates only can be found in FigureS1in the Supplementary Materials.

F I G U R E 2
Sobol indices for the parameters in the model.(a) Sobol indices with respect to TKI M .(b) Sobol indices with respect to TKI C .(c) Sobol indices with respect to TumorCellViability.IC 50 , half-maximal inhibitory concentration; I max , maximum inhibition; TKI, tyrosine kinase inhibitor.
, the tumor cell viability fluctuates as the concentration of the TKIs increases.Whereas the plots are based on the one specific cell line and specific concentrations, the data originate from different studies which have different cell viability response values at the same TKI concentrations.The model was able to capture the fluctuations as the individual estimates are based on the individual studies, from which the mean was taken to produce the simulated response.

F I G U R E 4
Predictions versus observations for tumor cell viability using individual parameter estimates.(a) Alectinib model.(b) Crizotinib model.(c) Erlotinib model.(d) Gefitinib model.The p values for the alectinib, crizotinib, erlotinib, and gefitinib treatments are, respectively, as follows: 1.068 × 10 −85 , 0*, 5.897 × 10 −169 , and 1.478 × 10 −184.*Note that a value of 0 does not show that a value is zero but rather that the value is smaller than the minimum value allowed to be computed in R.