Studies

Data used for model building consisted of five phase II/III trials and one phase I trial where CZP was administered subcutaneously to patients with moderate to severe RA (summarized in Table 1 [18, 31]). PK and ADA were mostly measured simultaneously, see Table 1. Overall, a total of 840 patients with 6898 CZP plasma concentration measurements and 6557 ADA measurements were available for the analysis. The proportion of patients that were reported as being ADA positive (ADA level exceeding the threshold) was 9.8% in the data used for model building. Informed consent was obtained from all study subjects. The studies were conducted in accordance with the applicable regulatory and International Council for Harmonisation–Good Clinical Practice requirements, the ethical principles that have their origin in the Declaration of Helsinki, and the local laws and regulations of the study sites.

Table 1 Data used in the analysis Observed variables

Two observed variables were considered; individually weighted PK residuals (\(__}\)) and ADA measurements (\(__}\)), to jointly incorporate the information from (i) the PK exposure of CZP as predicted from a model assuming no ADA production and (ii) information from the ADA assay, which is sensitive to the drug. EM-algorithms were used for estimation and parameters were MU-referenced.

To estimate the mode of \(__}\), the fit from a population PK model to the available CZP concentrations was used. The model was a one-compartment model including covariate effects of weight on CL/F and ethnicity (Japanese) on CL/F and apparent central volume (V/F). ADA was not part of the model. Inter-individual variability (IIV) was present on CL/F and V/F and a proportional residual error model was used. The model was initially fit (parameter estimation) to data from the first dosing occasion, thereafter the obtained individual parameter estimates predicted the remaining dosing occasions. Good model performance, i.e. unbiased residuals, is expected at the first dosing occasion, since patients were drug-naïve and unlikely to have formed ADA impacting the disposition of CZP at this time (last PK observation of first dosing occasion was at 23.6 days [average]). For future dosing occasions and stationary PK characteristics, predictions are also expected to be unbiased resulting in residuals centred around zero. However, in the presence of ADA, the model is expected to over-predict the observed CZP PK at subsequent dosing occasions, due to the absence of a covariate effect of ADA in the model, resulting in negative \(__}\)with a distribution deviating from the expected. Using the \(__}\) instead of the actual PK observations of CZP as an observed variable was motivated by two factors: \(__}\)have an expected distribution, which can be incorporated in the developed model, and, \(__}\) give an objective measure of when model performance indicates deviation from stationarity, expected when ADA are formed and not accounted for in the model.

In the studies included, ADAs were measured using a validated ELISA according to standards at the time of the bioanalysis [18, 31], since that time of analysis the ADA assay method has followed the lifecycle management and has been updated in later indications according to the state of art bioassay evolution for ADA detection. Subjects were classified as ADA clinically positive when at least one ADA sample was found to be above a set threshold of 2.4 U/mL. The threshold of 2.4 U/mL was derived based on its apparent clinical relevance where individuals in a phase II study with antibody measurements > 2.4 U/mL had lower CZP systemic exposure than those with measurements < 2.4 U/mL [18, 31]. Following this identification of clinically positive ADA, individuals were further classified as being persistent or transient ADA producers by individual subject profile inspection. The lower limits of quantification (LLOQ) of the CZP and ADA assays were 0.41 µg/mL and 0.6 U/mL, respectively, and all measurements that were below the quantification limit were set to their respective limits.

Hidden-Markov model development

A two-state MHMM related the observed variables (\(__}\) and \(__}\)), to the unobserved (hidden) underlying ADA production and the dependency of the observations on previous observations. A subject with a high measurement of ADA will likely have a high measurement at the next time point, thus a first-order Markov element was considered in this work. Production of ADA (SADA) and no production of ADA (SNOADA) were the two underlying states that were set as a representation of the immune response and the probability of transition between states were estimated. The transition probabilities sum up to 1 and only two elements of the transition probability matrix were estimated: probability of transitioning to SADA given that the previous state was SNOADA (πNOADA−ADA) and probability of transitioning to SNOADA given that the previous state was SADA(πADA−NOADA). he probabilities of remaining in the same state could be derived (πNOADA−NOADA and πADA−ADA.) All subjects were assumed to be drug naïve according to the studies criteria and therefore, the stationary distribution vector was not estimated, but rather the model was set up so that everyone started in SNOADA. The general description of an MHMM and its implementation in NONMEM has previously been described in Brekkan et al. [20]. The emission probabilities, which relate the probabilities of the observed variables to the hidden state at that time, were modeled as continuous random variables that could be correlated through a bivariate Gaussian probability density function (PDF):

$$\displaylines}}},}}}|S = s} \right) = \frac}}}^2\sigma _}}}^2(1 - \rho _s^2)} }}e \cr ^}\left( }}} - P}}}}}}}}}}}} \right)}^2} - 2\left( }}} - P}}}}}}}}}}}} \right)\left( }}} - AD}}}}}}}}}}}} \right) + }}} - AD}}}}}}}}}}} \right)}^2}} \right)} \cr}$$

(1)

where\(}_}_}}\) and \(}_}_}}\) are observed variables of interest, s is the current state (hidden), which can be either SNOADA or SADA, \(}_}\)and \(}_}_}}\) are state specific individual model predictions of the variables, \(}_}_}_}}}^\) and \(}_}_}_}}}^\) are state-specific variances of the variables (residual error) and ρs is the correlation between the variables. Two emission probability functions were required, one for each hidden state.

The model was extended to include IIV on several parameters. The individual values for the parameters were modeled assuming a normal distribution, exemplified for PKRES:

where i denotes an individual,\(__}}_\) is individual estimate of PKRES in a certain state (s), \(__}\) is the population estimate of the mode of the PKRES, and \(__}\) is a random effect assumed N(0, \(__}^\)) describing the deviation between the individual and typical values. A logit function was included on all probabilities to prevent individual probabilities exceeding 1 and then back converted to the probabilities using the logistic function, exemplified for πNOADA−ADA:

$$__}=\frac\text\text(\text\text\text\text}__}+ __})}\text\text(\text\text\text\text}__}+ __})}$$

(3)

where \(__}\) is the individual transition probability of going from SNOADA to SADA and \(__}\) is a random effect assumed N(0, \(__}^\)).

The population parameters to be estimated in the model and their expected estimates are presented in Table 2 and the model structure is presented in Fig. 1.

Table 2 Parameters in the developed bivariate hidden-Markov model Fig. 1

The general model structure of the bivariate hidden-Markov model given pharmacokinetic (PK) observations (\(__}\)) and anti-drug antibody measurements (\(__}\)). The two hidden states, SNOADA (grey) and SADA (orange) represent no production of ADA and production of ADA, respectively. The dashed vertical line separates the hidden part of the model from the observable part of the model. Observations were modeled using a bivariate Gaussian function. The probability of transitioning from SNOADA to SADA is πNOADA−ADA and the probability of transitioning from SADA to SNOADA is πADA−NOADA. The probabilities of staying in the respective states are given by πNOADA−NOADA and πADA−ADA. The dashed arrows represent the emission probabilities, i.e., the probabilities of the observations given the hidden state (for example \(P\left(__},__}|S=_\right)\) for the ADA producing state)

First, the bivariate MHMM for PKRES and ADAMES, was developed. Thereafter, the bivariate model was decoupled into two univariate models for the observed variables to determine the impact of incorporating both observed variables in the model. The individual state sequences resulting from the bivariate model and two univariate models were compared to the clinical ADA classification.

Parameter estimation and software details

Analysis of the data was performed in NONMEM and consisted of model estimation using maximum likelihood and importance sampling (IMP), followed by a post-hoc step in which the most likely state sequence in every individual was obtained using the Viterbi algorithm. The Viterbi algorithm was obtained as a downloadable subroutine for NONMEM available on the developer’s website [21, 22]. Perl-speaks-NONMEM (PsN) was used for execution and intermediate inspection of runs [23].

The likelihood of the first observation record in the model (assuming just one individual) is calculated as the probability of the observations (both at the same time in the bivariate model) in SNOADA (Eq. 1). The individual contribution of the states to the entire likelihood is kept and used for the next record. For the second record the likelihood is calculated as follows:

$$_=\frac_\bullet _+_\bullet _\right)\bullet P(_,_|S=_)}$$

(4)

$$_=\frac_\bullet _+_\bullet _\right)\bullet P(_,_|S=_)}$$

(5)

where P(PKRES, ADAMES) is the bivariate Gaussian PDF (in the bivariate model) and \(_\) and \(_\) are the transition probabilities. The individual contribution of the states to the likelihood, \(_\) and \(_\), are updated at each record and are defined as the entire likelihood divided by the likelihood of SNOADA and SADA, for the two respective states. Therefore, when estimating the model parameters, the likelihood of the data is maximized with respect to the transition probabilities and parameters describing observed variables.

Model evaluation and model predictions

Model parameter estimates were compared to their prior expectations (Table 2). Furthermore, model performance with regards to termination properties was considered, and models that terminated successfully in NONMEM without errors were considered to perform better than those that did not. Models with successful covariance steps resulting in parameter standard errors (SEs) were also considered to perform better than those that did not provide any parameter SEs. Individual state sequences were obtained from the MHMM and compared to the clinical classification of ADA in a random subset of individuals (n = 12) from study CDP870-004, since this study had the most consistent sampling within the individuals. Simulated (n = 100 simulations) distributions of the observed variables (\(__}\) and \(__}\)) from the final model were compared to the observed distributions.

In individuals that tested positive for ADA via the assay, the time to positive measurement was calculated and compared with the model predicted first time being in SADA. The best performing model, with regards to individual state prediction and estimation properties, was transferred into two univariate models and re-estimated. The performance of the resulting univariate models was compared with the bivariate model.