In JGWojtyniak/ospsreports: ospsreports Bookdown Workflow

General information {#chapgeneral}

This chapter gives short introductions into all processes and mechanisms relevant for the presented work. In addition, it provides general information regarding the system-dependent parameters, which were the same for all model simulations to ensure a comparable model environment. Furthermore, it specifies a standardized procedure for generating and evaluating virtual populations, as a necessary feature for implementing variability into \acrshort{pbpk}\index{PBPK} models. Moreover, this chapter describes the methods used for model development as well as for model evaluation.

Mathematical implementation of \acrfullpl{ddi} and \acrfullpl{dgi}

All required equations for the mathematical implementation of \acrshortpl{ddi}\index{DDI} and \acrshortpl{dgi}\index{DGI} are predefined in the \acrshort{osps} [@osps]. Hence, the purpose of this chapter is only to give a quick introduction into the mathematical background. For a more detailed description please have a look at the \acrshort{osps} documentation [@osps].

Competitive inhibition {#chapcompinhibition}

Competitive inhibition is the reversible binding of an inhibitor to the active site of an enzyme or transporter protein. Due to this, the inhibitor prevents that the native substrate can be metabolized or transportated by the enzyme or transporter, respectively. High concentrations of the native substrate are able to overcome this process (concentration-dependency). A special feature of competitive inhibition is that the \gls{vmax} remains unaffected, while the \gls{km} is apparently increased by the inhibition (equation \@ref(eq:eqcompetitivekmapp)). The (reduced) reaction velocity \gls{v} during co-administration of substrate and competitive inhibitor is described by equation \@ref(eq:eqcompetitivev):

\begin{loequation}{Competitive Inhibition} \begin{equation} \mathrm{KM_{app}=KM*(1+\frac{[I]}{Ki})} (#eq:eqcompetitivekmapp) \end{equation}

\begin{equation} \mathrm{V=\frac{(V_{max}[S])}{(KM_{app}+[S])}=\frac{\left(kcat [E]*[S]\right)}{(KM_{app}+[S])}} (#eq:eqcompetitivev) \end{equation} \end{loequation}

with \gls{vmax}=\glsdesc{vmax}, \gls{km}=\glsdesc{km}, \gls{kmapp}=\glsdesc{kmapp}, \gls{i}=\glsdesc{i}, \gls{s}=\glsdesc{s}, \gls{e}=\glsdesc{e}, \gls{ki}=\glsdesc{ki} and \gls{kcat}=\glsdesc{kcat}.

Mixed inhibition {#chapmixedinhibition}

In the case of mixed inhibition, the inhibitor can bind reversibly to the enzyme in a competitive manner or to the enzyme-substrate complex in an uncompetitive manner. Substrate and inhibitor have different binding sites on the enzyme. Equations \@ref(eq:eqmixedkm) and \@ref(eq:eqmixedvmax) describe the changes in \gls{kmapp} and \gls{vmax} in the presence of a mixed inhibitor:

\begin{loequation}{Mixed Inhibition} \begin{equation} \mathrm{KM_{app}=KM*\frac{1+\frac{[I]}{KI_{{c}}}}{1+\frac{[I]}{KI_{{u}}}}} (#eq:eqmixedkm) \end{equation}

\begin{equation} \mathrm{V_{max,app}=\frac{V_{max}}{(1+\frac{[I]}{KI_{{u}}})}} (#eq:eqmixedvmax)] \end{equation} \end{loequation}

with \gls{vmaxapp}=\glsdesc{vmaxapp}.

\acrfull{mbi} {#chapmechanismbased}

If an enzyme is irreversibly inhibited, de-novo synthesis is required to restore the initial situation regarding enzyme activity (time-dependency). A special case is the \acrshort{mbi}\index{MBI}, where an enzyme is irreversibly inhibited by a reactive species formed during the catalytic process. The change in enzyme turnover in the case of \acrshort{mbi} can be described by equation \@ref(eq:eqmixedinhibition):

\begin{loequation}{Mechanism-Based Inactivation} \begin{equation} \mathrm{\frac{d[E]{cat}(t)}{dt}=k{deg}E_0-(k_{deg}+\frac{k_{inact}I}{(K_i+I)}*[E]_{cat}(t))} (#eq:eqmixedinhibition) \end{equation} \end{loequation}

with \gls{ecatt}=\glsdesc{ecatt}, \gls{kdeg}=\glsdesc{kdeg}, \gls{e0}=\glsdesc{e0}.

Induction {#chapinduction}

Enzyme and transporter induction is caused by a drug-induced increase of gene expression which can be described by \@ref(eq:eqinduction):

\begin{loequation}{Induction} \begin{equation} \mathrm{\frac{d[E]{cat}(t)}{dt}=k{deg}E_0(\frac{1+(E_{max}*[Ind])}{EC_{50}+[Ind]})}) (#eq:eqinduction) \end{equation} \end{loequation}

with \gls{emax}=\glsdesc{emax}, \gls{ind}=\glsdesc{ind}, \gls{ec50}=\glsdesc{ec50}.

\acrshortpl{dgi}

For all \acrshortpl{dgi} implemented a change in the transporter or enzyme \gls{kcat} was assumed compared to wildtype. Although, also a change in \gls{km} would be conceivable for the polymorphisms included in this network, no evidence for this situation could be found in literature. Besides, \acrshortpl{dgi} were implemented in a stepwise procedure. First, the wildtype \gls{kcat} was optimized based on homozygous wildtype individuals. Afterwards, based on homozygous deficient individuals the enzyme or transporter \gls{kcat} for the polymorphism was estimated. Finally, the implementation was evaluated using heterozygous individuals. For this purpose the \gls{kcat} for individuals with heterozygous genotypes were calculated according to equation \@ref(eq:dgis) assuming an additive relationship between the wildtype and the deficient pharmacogene.

\begin{loequation}{\acrshort{dgi} calculation for heterozygous individuals} \begin{equation} k_{cat,heterocygous}=\mathrm{0.5k_{cat,wildtype}+0.5k_{cat,deficient}} (#eq:dgis) \end{equation} \end{loequation}

General modelling workflow overview

According to the recently released \acrfull{ema} guideline: "Guideline on the reporting of physiologcially based pharmacokinetic (PBPK) modelling and simulation." [@EuropeanMedicines2018] the use of an analysis plan is recommended. For this reason, a modelling workflow plan was developed which describes the whole \acrshort{ddi} and \acrshort{dgi} network establishment (\@ref(fig:figgeneralworkflow)). In addition the R package drake was used for the establishment of an reproducible statistical analysis plan.

loadd(general_workflow_plot)
general_workflow_plot

Metabolizing enzymes {#chapmetenzymes}

The following chapters are a brief summary of the model relevant biological background information as well as relevant enzyme and transportation processes:

\acrfull{cyp} {#chapcyp}

\acrshortpl{cyp} accounting for round about 75% of the total drug metabolism and following are the major enzymes involved in drug metabolism [@Guengerich2008]. Figure \@ref(fig:figfractioncyp) shows the fractions of clinically used drugs metabolized by the corresponding \acrshort{cyp} isoforms according to @Zanger2013. In the presented work especially the isoforms \acrfull{cyp3a4}, \acrfull{cyp3a5}, \acrfull{cyp2c8} and \acrfull{cyp2c9} are of relevance.

loadd(cyp_expression_plot)
cyp_expression_plot

Drug dissolution {#chapdissolution}

For drug dissolution one of the predefined formulations were used, depending on which template could best describe the observed data. Possible predefined formulations includes:

• Dissolved drug
• Weibull dissolution
• Particle dissolution
• Dissolution table

For a full overview see the \acrshort{osps}-manual [@osps].

Meal time events {#chapmealtime}

Unless otherwise described in the substance-specific chapters for meal time events the standard implementation as described in the \acrshort{osps}-manual was used [@osps]. If information regarding the meal composition was available they were included in the meal time events.

System-dependent parameters {#chapsystemparams}

All system-dependent parameters like reference concentrations, protein half-lives as well as tissue expression profiles of transporters and relevant metabolizing enzymes are listed in Table \@ref(tab:tabsystemparam). The relative expression in different tissues is visualized in \@ref(fig:)

loadd(system_dependent_parameter_table)
system_dependent_parameter_table

loadd(relative_expressions_plots)
subcaption="Enzyme or Transporter"
split_sub_fig(plots=relative_expressions_plots, chunk_caption="RelativeExpression", caption = "Relative enzyme or transporter expression in various tissues", maxlength = 1,subcaption=subcaption,outwidth="'100%'",ncol=1,figwidth=7,figasp=1.618,dev="'CairoPDF'",dev_args=list(NULL))

In all individuals and populations, the \acrfull{ehc}\index{EHC} was enabled (\acrshort{ehc} continuous fraction set to one) by assuming a continuous flow from bile to duodenum. Furthermore, if no information about the renal state were specified a glomerular filtration rate of one was assumed which is aequivalent to passive filtration without active secretion or reabsorption.

Virtual populations {#chapvirtualpop}

To cover the variability of simulated concentration-time profiles, virtual populations containing 1000 individuals with demographic properties (gender, age, weight, height, race) adapted to the range of each respective study were created based on the \acrfull{icrp}\index{ICRP} database. Furthermore, if \glsdesc{mean} and \acrfullpl{sd} were stated for relevant study demographics, the area intersected by the simulated and observed probability density functions were calculated as a simplified check for equality of the simulated and observed populations. If no \acrshort{sd} but the minimum and maxium values were stated, \acrshort{sd} were roughly estimated after \@ref(eq:popdist).

\begin{loequation}{Calculate \acrshort{sd} based on parameter range} \begin{equation} \mathrm{SD=\frac{(min+max)}{4}} (#eq:popdist) \end{equation} \end{loequation}

Besides, if no information on the gender of the volunteers was available, a 100% male population was assumed. To model studies that did not report the age range of their individuals, a population between 20 and 50 years of age was created. The reference concentrations of proteins for which ontogeny information were available in PK-Sim were distributed according to the variabilities provided in the PK-Sim database. If no information was available in the modeling platform, they were set to be log-normally distributed according to literature reports or otherwise with a moderate geometric \acrshort{sd} of 1.4 (35% \acrfull{cv}).

Model development {#chapmodeldevelopment}

For model development, concentration-time profiles of published clinical studies covering the full reported dosing range were used. If available this included data of intravenous and oral administration after single or multiple doses. If information about the amount of unchanged drug excreted to urine or feces were available, they were also utilized to inform the model optimization. If parameters had to be identified the Monte-Carlo algorithm implemented in PK-Sim were used. Parameter identifications were performed if either no literature values were on hand, or multiple values with a broad range were available or the parameter seemed to have a major influence on the drugs pharmacokinetic profile. Besides, model relevant information regarding all clinical studies available are listed in the corresponding chapters (\@ref(chap)) including information about the assignment of each study to the training (model building) or test data set (model evaluation). To validate, that these information are equally distributed within the training and test datasets, \acrfullpl{anova}\index{ANOVA} or chi-square goodness-of-fit tests were performed for continuous and categorical data, respectively. For model development for each study in the training dataset a mean prediction was used based on the mean demographic properties (gender, age, weight, height, race) of the respective study. If values were missing, a mean value for age (30 years), weight (73 kg) and height (176 cm) were used. Parameters of the final models are given in Tables \@ref(tab:parameterdataset). Previously developed PBPK models used for the DDI network development were extended as described in the chapters \@ref(chapmodelextensions). Parameter sets of all model relevant parameters of extended models are given in \@ref(chapmodelextensions).

Model evaluation {#chapevaluation}

For all observed concentration-time profiles population simulations using the virtual populations as described in chapter \@ref(chapvirtualpop) were generated. The predicted profiles were subsequently used for the statistical as well as the graphical model evaluation. For this purpose, simulated \acrshortpl{mean} and \acrshortpl{sd} were calculated as the observed data were mostly reported by this descriptive measures. Following, various \acrfullpl{nca}\index{NCA} parameter like the \acrfull{auc}\index{AUC} and the \gls{cmax}\index{$C_{max}$} were calculated and used for further model evaluation.

Statistical model evaluation {#chapstatevaluation}

For statistical model evaluation the following accuracy measures were calculatet. For all \acrshort{mean} and \acrshort{sd} concentration-time values available the \acrfull{mrd}\index{MRD} as well as the \acrfull{msa}\index{MSA} according to equation \@ref(eq:mrd) and \@ref(eq:msa) were calculated.

\begin{loequation}{\acrfull{mrd} calculation method} \begin{align} \mathrm{x=\sqrt{\frac{\sum_{i=1}^{N}(\log_{10}c_{obs}-\log_{10}c_{pred})^{2}}{N}}} \notag \ \mathrm{MRD=10^{x}} (#eq:mrd) \end{align} \end{loequation}

\begin{loequation}{\acrfull{msa} calculation method} \begin{align} \mathrm{Q=\frac{c_{pred}}{c_{obs}}}\notag \ \mathrm{\zeta=100*(\exp(\operatorname{median}(|\log_{e}(Q_{i}|))-1)} (#eq:msa) \end{align} \end{loequation}

Hereby a \acrshort{mrd}\index{MRD} value $≤$ 2 characterizes an adequate prediction whereas the \acrfull{msa}\index{MSA} gives a easy interpretable impression of the relative deviation of the model predictions compared to the observed values. Additionally, \acrshort{nca} ratios (equation \@ref(eq:ncaratio)) were calculated and compared.

\begin{loequation}{\acrshort{nca} ratios observed versus predicted} \begin{equation} \mathrm{NCA_{ratio}=\frac{NCA,value_{pred}}{NCA,value_{obs}}} (#eq:ncaratio) \end{equation} \end{loequation}

Again a value $≤$ 2 is considered as sufficient. Further, for each \acrshort{ddi} and \acrshort{dgi} observed and predicted \acrshort{nca} effect ratios were estimated (see equation \@ref(eq:effectratio)).

\begin{loequation}{\acrshort{ddi} and \acrshort{dgi} effect ratios} \begin{equation} \mathrm{Effect_{ratio}=\frac{\frac{NCA,value ^{DDI,DGI}{pred}}{NCA,value^{placebo}{pred}}}{\frac{NCA,value^{DDI,DGI}{obs}}{NCA,value^{placebo}{obs}}}} (#eq:effectratio) \end{equation} \end{loequation}

Finally, as a quantitative measure of the prediction accuracy for each DDI interaction as well as for all placebo concentration-time profiles the \acrfull{gmfe}\index{GMFE} were calculated for all observed \acrshort{nca} values according to \@ref(eq:gmfe):

\begin{loequation}{\acrfull{gmfe} calculation method} \begin{align} \mathrm{GMFE=10^{x}} \notag \ \mathrm{x=\frac{\sum|\log_{10}(\frac{NCA,value_{pred}}{NCA,value_{obs}})|}{n}} (#eq:gmfe) \end{align} \end{loequation}

Graphical model evaluation {#chapgraphevaluation}

For graphical model evaluation a set of different figures were generated. For all concentration-time profiles \acrfullpl{vpc} were created using the virtual populations as described in chapter \@ref(chapvirtualpop) and the subsequently calculated \acrshort{mean} and \acrshort{sd} profiles. Furthermore, for all data included in the training dataset \acrfull{gof} plots like observed values versus predicted values, time versus residuals or predicted values versus residuals were generated (see \@ref(fig:gof)). Therefor, the mean profiles as described in chap \@ref(chapmodeldevelopment) were used. Moreover, predicted versus observed \acrshort{nca} ratios as well as \acrshort{ddi} and \acrshort{dgi} effect ratios as calculated in chap \@ref(chapstatevaluation) were evaluated using the twofold limits, halffold limits and the limits proposed by Guest et al [@guest].

Parameter correlation, covariance and \acrshort{nca} parameter sensitivity

For all parameters which had to be optimized the partial derivatives calculated locally at the found optimal parameter values were used to calculate a correlation and covariance matrix which give some information about local sensitivity of the parameters included in the \acrshort{pi} [@osps]. Afterwards these information were used to visualize the correlation and covariance matrices and to estimate confidence intervals for the identification parameters [@osps].

Additionally, \acrshort{nca} parameter sensitivities of the final PBPK models were calculated measured as relative changes of the \acrshort{auc}, \acrshort{cmax} and \acrshort{tmax} of one dosing interval in steady-state conditions for simulations of the highest recommended doses. For this purpose, a mean individual as described in chapter \@ref(chapmodeldevelopment) was used. Parameters were included into the analysis if they had to been optimized, if they might have a strong influence due to calculation methods used in the model (fraction unbound) or if they had significant impact in former models (solubility, blood/plasma ratio, glomerular filtration rate fraction).

The sensitivity for the \acrshort{nca} parameter to the input parameter of interest was then calculated as the ratio of the relative change of the \acrshort{nca} parameter and the relative variation of the input parameter (see equation). For reasons of numerical stability, sensitivities were calculated as the average of several sensitivities based on different variations (see equation). The relative variations are defined by multiplication of the value in the simulation with variation factors. These variation factors are defined by equation. Hereby, the variation range was choosen in accordance with the recommendation from the current \acrshort{ema} guideline [@EuropeanMedicines2018].

\begin{loequation}{\acrshort{nca} parameter sensitivity} \begin{equation} \mathrm{S_{i,j}=\frac{\Delta PK_j}{\Delta p_i}\frac{p_i}{PK_j}} (#eq:sensitivity1) \end{equation} \begin{equation} \mathrm{S_{i,j}=\frac{\sum_{k=1}^n\frac{\Delta_k PK_j}{\Delta_k p_i}\frac{p_i}{PK_j}}{n}} (#eq:sensitivity2) \end{equation} \begin{align} \mathrm{(1+\alpha\frac{k}{n}}\notag \ \mathrm{(\frac{1}{1+\alpha\frac{k}{n}})} (#eq:sensitivity3) \end{align} \end{loequation}

JGWojtyniak/ospsreports documentation built on Nov. 14, 2019, 12:11 p.m.