The analysis was performed using pooled data of patients treated with avelumab across drug development phases with different oncology indications from multiple clinical trials. An overview of the analyzed patient population is given in Table 1. Data from 1892 patients with evaluable PK data were used in the analysis and were split into two subsets stratified by ADA incidence (ever vs never ADA+): a subset for model development (1513 patients, ~ 80% of data) and a validation subset (379 patients, ~ 20% of data). Previously, a population PK model developed on data from 3 of the 7 clinical trials was published [11].

Avelumab concentrations were quantified using an immunoassay sandwich method (interrun coefficient of variation [CV] ≤ 16.1%, bias [absolute value] ≤ 15.5%, total interrun error ≤ 19.1%). The lower limit of quantification was 0.2 μg/mL [11]. ADA testing was conducted using a tiered assay approach, whereby samples that were screened and confirmed positive for ADA were subsequently analyzed to determine the titer using the homogeneous bridging electro chemiluminescent immuno-assay [12]. The following numbers of data points/participants were used for the model development and evaluation: PK development, 12,707/1513; PK evaluation, 3362/377; CTMM development, 14,874 /1495; CTMM evaluation, 3925/375; DTMM development, 14,247/1513; and DTMM evaluation, 3766/378. The average number of ADA samples and collection times are provided in Table 1. Missing ADA data were ignored in the analyses. There was no minimum ADA data requirement for a patient to be considered evaluable. Only 7 patients with neutralizing antibody were reported, and for 76 data points, the neutralizing antibody information was missing.

A previously-developed population PK model for avelumab formed the basis for PK-ADA model development [11]. This model was 2-compartmental, with time-dependent clearance (CL; Fig. 4). Decrease in CL was modeled as a sigmoid maximal inhibitory response process [13]:

$$C_=TVCL\cdot }\left(\frac_\cdot _^}_^+_^}\right)\cdot }(_)$$

Here, \(_\) is CL in individual \(i\) at time \(t\), \(TVCL\) is the typical value of CL in the population, \(_\) is the maximal possible change in CL relative to baseline for individual \(i\), \(_\) is the time after first dose in individual \(i\), \(_\) is the time at which 50% of \(_\) is reached, \(\gamma\) describes the shape of the relationship, and \(_\) is interindividual variability (IIV) in CL for individual \(i\), defined as being normally distributed with a mean of 0 and variance of \(_^\). In addition to CL, IIV was included on central volume of distribution (V1) and \(_\) (additively for the latter). The residual error was described by a combined additive and proportional error model. For the current analysis all covariate effects were removed except body weight as a covariate on CL and V1 using standard allometry coefficients [14].

Two approaches to model ADA data were investigated: (1) ADA as a covariate in the population PK model, and (2) Markov models of ADA status. Finally, a PK-ADA model linked the Markov model for ADA status to the population PK model to assess the bidirectional effect of ADA status on PK and drug exposure on ADA status simultaneously. The modeling strategy started with model development for models with lower complexity/computational burden and progressed to models with higher complexity/computational burden. ADA status was evaluated for effect on avelumab PK as a categorical covariate on avelumab clearance using different categorizations of ADA status: ADAEVER (time-invariant, ever vs never ADA+), ADAONCE (time-varying, ADA+ at first appearance of ADA and thereafter), and ADALOCF (time-varying, ADA+ vs ADA−) as graphically demonstrated in Fig. 1.

Categories of time-variant ADA, using a representative patient

The ADA part of the joint model was considered be Markovian in nature, since it is likely to depend on state transitions over time: ADA− to ADA+ , and back again, based on the previous state, and on time. ADA status was modeled using two-state discrete-time (DTMM) and two-state continuous-time (CTMM) Markov models. The Markov ADA-models were evaluated using VPCs for both the probability of ADA status as well as the probability to change ADA state or to remain in the same ADA state. Covariate relationships identified for the DTMM were included in the CTMM. The joint PK-ADA model was constructed by linking the population PK model to the CTMM to assess the effects of ADA on avelumab clearance and the effect of avelumab exposure on ADA status.

Discrete-time Markov models (DTMMs) assume time to be a discrete variable, the recorded observation times in a clinical study. DTMM equations describe the transition probability from state \(_\) to state \(_\), or more simply, the probability of each state \(_\) at the current occasion \(j\) given that state \(_\) was observed at the previous occasion \(j\)−1 [15].

The DTMM (Fig. 2) may be summarized as follows:

Structure of the discrete-time Markov model. \(_:\) probability to transition from ADA+ to ADA−; \(}}_\): probability to transition from ADA− to ADA+

$$_=P(_=\left.1\right|_=0)$$

$$_=P\left(_=\left.0\right|_=0\right)=1-_$$

$$_=P(_=\left.0\right|_=1)$$

$$_=P(_=\left.1\right|_=1)=1-_$$

Here, \(_\) is the probability of a transition from ADA− (state 0) to ADA+ (state 1) in individual \(i\) at observation occasion \(j\), \(_\) is the probability of a transition from ADA+ to ADA− in individual \(i\) at observation occasion \(j\), \(_\) is the probability that ADA remains negative in individual \(i\) at observation occasion \(j\), and \(_\) is the probability that ADA remains positive in individual \(i\) at observation occasion \(j\). \(_\) is the observed value of ADA (+ or −) in individual \(i\) at observation occasion \(j\), while \(_\) is the observed value of ADA (+ or −) in individual \(i\) at preceding observation occasion \(j\)−1.

The probabilities were estimated using a logit transformation:

$$P_} = \frac} + \eta_}} } \right)} \right)}} } \right)}}$$

Here, \(_\) is the individual realization of patient \(i\)’s probability to transition from ADA− to ADA+ where \(_\) is the typical value of the probability in the population, and \(__}\) is the IIV defined as being normally distributed with a mean of 0 and variance of \(__}^\).

The continuous-time Markov model (CTMM) is more appropriate for non-uniform observation times, although for two-state data, the DTMM and CTMM methods may produce similar results for naturally discretized data (e.g. prespecified time points based on clinical trial designs) [10]. In this method, the influence of the previous state on the probability of the current state typically decreases with increasing time between observations [16,17,18].

$$P_} = \frac} + \lambda _} \cdot e^} + \lambda _} } \right)\;t_} }} }}} + \lambda _} }}$$

$$P_} = \frac} + ~\lambda _} \cdot e^} + \lambda _} } \right)\;t_} }} }}} + \lambda _} }}$$

Here, \(_\) is the transfer rate constant from state 0 to state 1 in individual \(i\), \(_\) is the transfer rate constant from state 1 to state 0 in individual \(i\), and \(_\) is the time difference between the current and previous observations in individual \(i\) at observation occasion \(j\). Other parameters are as previously defined. Treatment and covariate effects can be applied to \(_\) and \(_\).

In the current case, the CTMM was implemented as ordinary differential equations (ODEs) representing the probabilities of ADA being negative or positive (Fig. 3):

Structure of the continuous-time Markov model (CTMM). \(_\): rate constant for ADA+ to ADA− transition; \(_\): rate constant for ADA− to ADA+ transition; \(A1\): probability for ADA+ ; \(A0\): probability for ADA−

$$\frac_}=-_\cdot _+_\cdot _$$

$$\frac_}=_\cdot _-_\cdot _$$

\(_\) and \(_\) represent the probabilities of ADA being negative and positive, respectively. Using this approach, modifications to the data file were necessary to facilitate analysis with NONMEM (version 7.4.3) [19]; compartment amounts and associated bioavailability values were reset after each observation [18]. Aside from ADA status itself, covariates on ADA state transition probabilities were considered, particularly tumor type. The transition rates were estimated using a log transformation:

\(_\) and \(_\) represent the probabilities of ADA being negative and positive, respectively. Using this approach, modifications to the data file were necessary to facilitate analysis with NONMEM (version 7.4.3) [19]; compartment amounts and associated bioavailability values were reset after each observation [18]. Aside from ADA status itself, covariates on ADA state transition probabilities were considered, particularly tumor type. The transition rates were estimated using a log transformation:

Here, \(_\) is the individual realization of patient \(i\)’s transition rate from ADA− to ADA+ where \(__}\) is the typical value of the transition rate in the population and \(__}\) is the IIV defined as being normally distributed with a mean of 0 and variance of \(__}^\).However, the data did not support the estimation of IIV in \(_\) or \(_.\)

The joint PK-ADA model combined the base population PK model (without ADA effect or IOV) with the final CTMM (including covariates). The models were linked by (1) estimating the correlation between avelumab clearance and \(_\), (2) effect of probability of ADA+ status (\(p1\)) on clearance (mono-directional joint PK-ADA model), and (3) time-course of the effect of avelumab exposure on the \(_\) and the effect of \(p1\) on clearance (bi-directional joint PK-ADA model). The typical values of avelumab clearance, \(_\), and the associated IIV parameters were estimated. All other parameters were fixed to the final typical values of the estimates of the base population PK model and the final CTMM.

The bi-directional effects are demonstrated in Fig. 4. The effect of \(p1\) on avelumab CL (\(_\)) was modeled as a linear or as a non-linear relationship between \(p1\) (independent variable) and CL (dependent variable):

Structure of the bidirectional joint PK-ADA model. CL: clearance; Q: intercompartmental clearance; V1: central volume of distribution; V2: peripheral volume of distribution; \(P_\): effect of probability of ADA+ on CL; \(C_\): effect of concentration in central compartment on probability of ADA+ ; \(_\): effect of time on CL; \(_\): maximum change in CL relative to baseline; \(_\): time for half of maximum effect; \(\gamma :\) shape of time effect curve; \(_\_c1}\): maximum change in effect of exposure on rate of transition from ADA− to ADA+ ; \(_\): concentration for half of maximum exposure effect; \(_\): shape of exposure effect curve; \(_\_p1}\): maximum change in CL due to ADA; \(_\): probability of ADA+ for half of maximum ADA effect; \(_\): shape of the ADA effect curve; \(_\): rate constant for ADA+ to ADA− transition; \(_\): rate constant for ADA− to ADA+ transition; \(p1\): probability for ADA+ ; \(p0\): probability for ADA−

$$_=\frac_}_}\cdot _\cdot ^_+_)}$$

$$_=1+sl\cdot _ (\text)$$

$$_=1+\frac_\_P1 }\cdot _^_}}_^_}+_^_}} (\text)$$

Here, \(_\), \(_\), and \(_\) are the avelumab elimination rate constant from the central compartment, avelumab CL, and central volume of distribution, respectively, for patient \(i\). In the linear model the intercept is set to 1 (no change) and \(sl\) is the slope of change in \(p1.\) The nonlinear model is a sigmoidal Emax model, and \(_\_P1}\) and \(_\) are the maximum effect and 50% of the maximum effect of \(p1\) on clearance. \(^\) describes the shape of the relationship.

The link model for the relationship between exposure and rate of ADA transition (\(_\)) is a sigmoidal Emax model:

$$_=1+\frac_\_C1\bullet _^_}}}_^_}+_^_}}$$

Here, \(_\) is avelumab concentration in the central compartment in individual \(i\) at time \(t\), \(_\_C1\bullet }\) is the typical value of the maximal possible change in \(_\) relative to baseline, \(_\) is the typical value of the concentration at which 50% of \(_\_C1\bullet }\) is reached, and \(^\) describes the shape of the relationship. The direction of the maximum change (i.e. an increase or a decrease) was not specified but was limited to a range of − 100% to 500%.

Potential covariate relationships were explored graphically by plotting the potential covariates versus the parameters of interest. Graphical exploration procedures were only relied upon if the degree of η-shrinkage in the parameters was reasonably low (< 30%) [20].

Covariate analysis for the population PK model was restricted to the relationship between the different categorizations of ADA status (ADAEVER, ADAONE, ADALOCF, etc.) and avelumab CL. The ADA-CL relationship was implemented using a linear function:

$$_=TVCL\cdot \left(1+_\right)$$

Here, \(_\) is CL in individual \(i\), \(TVCL\) is the typical value of CL in the population, and \(_\) is the estimated of change in CL when ADA+ . For the DTMM and CTMM, covariate relationships included demographics (age, sex, and race), disease-related status (serum albumin, C-reactive protein [CRP], Eastern Cooperative Oncology Group performance score [ECOG] status, tumor type, and baseline tumor burden), previous and current concomitant therapies (previous use of biologics or PD-L1 inhibitors), and baseline laboratory results and organ function (aspartate transaminase [AST], alanine transaminase [ALT], creatinine, creatinine clearance, bilirubin, estimated glomerular filtration rate [eGFR]). Categorical covariates were tested using a linear function and continuous covariates were tested using power functions (see below).

Covariate model development was performed stepwise (using the stepwise covariate modeling tool in Perl-speaks-NONMEM [PsN], version 5.3.0) [21, 22]. Each candidate covariate was tested on each of the parameters of interest, one at a time. The parameter-covariate relationship producing the largest change in the NONMEM objective function value (OFV) was retained. This process was repeated as a series of forward model-building steps until no further parameter-covariate relationships were present that met the forward inclusion criterion (a change in OFV of − 3.843, corresponding to a nominal significance level of p = 0.05). A backward elimination process was then undertaken, in which each relationship was removed one at a time. At each backward step, the parameter-covariate relationship with the lowest change in OFV and not meeting the backward elimination criterion (a change in OFV of + 10.83, corresponding to a nominal significance level of p = 0.001) was removed. The process was concluded when no further parameter-covariate relationships could be removed.

The testing of categorical covariates was implemented using a linear function, as follows:

$$_=_\cdot \left(1+_\right)$$

where \(_\) is the parameter value for individual \(i\), \(_\) is the typical value of the parameter in the population, and \(_\) is an estimated parameter corresponding to the unique value of the categorical covariate in individual \(i\). For the largest or reference category, \(_\) was defined as 0. Covariate categories containing less than 20 patients were not separately tested (except for tumor type) but instead lumped with the reference case (typically the category with highest frequency in the population).

Testing of continuous covariates was performed using a power function, as follows:

$$_=_\cdot _} }\right)}^_}$$

where \(_\) and \(_\) are as previously defined, \(_\) is the value of the covariate in individual \(i\), \(\overline\) is the median value of the covariate in the population, and \(_\) is a parameter describing the shape of the relationship of the covariate to the parameter.

Stratification was used when appropriate to ensure that the model could be evaluated adequately across important subgroups of the data. The adequacy of the models was evaluated using visual predictive checks (VPC) [23]. The population PK and PK-ADA models were used to simulate 500 replicates (model development) or 1000 replicates (external validation) of the analysis data set. Statistics of interest were calculated from the simulated and observed data for comparison: the 2.5th, 50th (median), and 97.5th percentiles of the distributions of the simulated concentration at each sampling time bin were calculated. These percentiles of the simulated data were plotted versus time, with the original observed data set and/or percentiles based on the observed data overlaid to visually assess concordance between the model-based simulated data and the observed data.

Final models were used to predict dependent variables (avelumab concentrations and/or ADA data) based on the patients in the data set who were set aside for external validation (20% of the data set) at the individual and population level (the latter using VPCs). These were compared with observations to provide an assessment of the model’s predictive ability.