Characterization of anti-drug antibody dynamics using a bivariate mixed hidden-markov model by nonlinear-mixed effects approach

Biological therapies may act as immunogenic triggers leading to the formation of anti-drug antibodies (ADAs). Population pharmacokinetic (PK) models can be used to characterize the relationship between ADA and drug disposition but often rely on the ADA bioassay results, which may not be sufficiently sensitive to inform on this characterization. In this work, a methodology that could help to further elucidate the underlying ADA production and impact on the drug disposition was explored. A mixed hidden-Markov model (MHMM) was developed to characterize the underlying (hidden) formation of ADA against the biologic, using certolizumab pegol (CZP), as a test drug. CZP is a PEGylated Fc free TNF-inhibitor used in the treatment of rheumatoid arthritis and other chronic inflammatory diseases. The bivariate MHMM used information from plasma drug concentrations and ADA measurements, from six clinical studies (n = 845), that were correlated through a bivariate Gaussian function to infer about two hidden states; production and no-production of ADA influencing PK. Estimation of inter-individual variability was not supported in this case. Parameters associated with the observed part of the model were reasonably well estimated while parameters associated with the hidden part were less precise. Individual state sequences obtained using a Viterbi algorithm suggested that the model was able to determine the start of ADA production for each individual, being a more assay-independent methodology than traditional population PK. The model serves as a basis for identification of covariates influencing the ADA formation, and thus has the potential to identify aspects that minimize its impact on PK and/or efficacy. Supplementary Information The online version contains supplementary material available at 10.1007/s10928-023-09890-8.


Introduction
The administration of biological drugs can invoke an immune response where anti-drug antibodies (ADAs) are formed [1].Since ADA formed against biological drugs may cause treatment failure, characterization of immunogenic response is important.Formed ADAs may differ with respect to their clinical impact, e.g., not all binding antibodies are neutralizing and some may affect the disposition of biologics via formation of drug-ADA complexes, which may be eliminated differently from the unbound drug [2].ADA influence, e.g.ADA are reliably measured but no characterization to classify them as neutralizing available, in such case the occurrence of ADA may add to explain variability in PK but may not have use for clinical interpretation due to the lack of information on whether the ADA is neutralizing or not [8,9].Even if improvements in assay techniques result in an increase in the sensitivity of the assay to detect ADA, the priority should be to identify ADA with a clinical impact and to determine thresholds for ADA levels that are clinically relevant.Such determinations can often be susceptible to bias, particularly when the frequency of positive ADA is low and common statistical methods cannot be used to determine a clinically relevant threshold [10].
Methods to characterize ADA dynamics and patient response to ADA are of value given the potential consequences on clinical outcome.Population pharmacokinetic (PK) and pharmacodynamic (PD) analyses are frequently used to characterize the disposition of drugs and several have been published for anti-TNFα therapeutics [11][12][13].
To account for immunogenicity in such analyses, ADAs are often included as covariates on the clearance (CL) of the therapeutic.However, due to the aforementioned assay aspects, some ADA measurements may be classified as being false positive or false negative.Further, patients are commonly not considered ADA positive until the ADA measurement reaches a threshold level, afterconsequences for drug disposition may already have occurred.Therefore, alternative methods accounting for the masked (or hidden) nature of ADA dynamics and discerning ADA production that influences the PK of the drug may be of interest.
Mixed hidden-Markov models (MHMMs) can be used to determine the relationship between observed and unobservable stochastic processes on a population level [14].MHMMs allow for the characterization of unobservable processes in patients using one or more observed variables.They have been applied in various areas, including drawing inference about relapse in multiple sclerosis (hidden variable) given observations of brain lesions in patients and determining exposure-response in epilepsy where epileptic activity is the hidden variable [14,15].An MHMM using measurements of PK and ADA to inform the immunogenic response (hidden variable), may enable the characterization of the underlying ADA production in patients [16].
The aim of this work was to explore a novel method to characterize ADA dynamics, through the development of a bivariate MHMM accounting for drug PK and ADA measurements to predict the most likely ADA dynamics in an individual, using certolizumab pegol (CZP) data as an example.CZP is a PEGylated Fc-free antibody fragment, which minimizes its potential for complement-dependent cytotoxicity and antibody dependent cell-mediated cytotoxicity [17].Pegylation combats the relatively short half-life caused by the removal of the Fc region, resulting in a welltolerated molecule that has a half-life of approximately 14 days [18].CZP binds to both soluble and membrane-bound TNFα, blocking the TNFα mediated inflammatory cascade and reducing the clinical symptoms of RA and other chronic inflammatory diseases [19].In the current work, historical ADA binding data and PK was available based on the studies used.New ADA binding and PK assays have been developed for to test the drug in other indications, according to the state of the art and lifecycle management of the product however, for the purpose of the model development, the former historical assays were considered appropriate.

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.

Observed variables
Two observed variables were considered; individually weighted PK residuals (YPK RES ) and ADA measurements (YADA MES ), 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 Y P K RES , 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 Y P K RES with a distribution deviating from the expected.Using the Y P K RES instead of the actual PK observations of CZP as an observed variable was motivated by two factors: Y P K RES have an expected distribution, which can be incorporated in the developed model, and, Y P K RES 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 (YPK RES and Y ADA MES ), 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 (S ADA ) and no production of ADA (S NOADA ) 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 S ADA given that the previous state In studies CDP870-004, CDP870-011 and CDP870-014, no loading dose was administered.Study PHA001 was a single dose study.In studies RA006 and C87041, loading doses were administered at weeks 0, 2, and 4 and followed by the maintenance doses First, the bivariate MHMM for PK RES and ADA MES , 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 S NOADA (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: (5) where P(PK RES , ADA MES ) is the bivariate Gaussian PDF (in the bivariate model) and π ADA−NOADA and π ADA−ADA are the transition probabilities.The individual contribution of the states to the likelihood, Φ NOADA and Φ ADA , are updated at each record and are defined as the entire likelihood divided by the likelihood of S NOADA and S ADA , 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 was S NOADA ( π NOADA−ADA) and probability of transitioning to S NOADA given that the previous state was S ADA (π 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 S NOADA .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): The model was extended to include IIV on several parameters.The individual values for the parameters were modeled assuming a normal distribution, exemplified for PK RES : where i denotes an individual,P K RESs i is individual estimate of PK RES in a certain state (s), P K RESs is the population estimate of the mode of the PK RES, and η P K RES is a random effect assumed N(0, ω 2 P K RES ) 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 : where π NOADA−ADA i is the individual transition probability of going from S NOADA to S ADA and η π NOADA−ADA is a random effect assumed N(0, ω 2 π NOADA−ADA ).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.
However, in study PHA001 no ADA measurements were available as it was a single dose administration study.
Bivariate models with different degree of complexity were developed: 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 S ADA .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.

Results
A total of 6898 CZP plasma concentration measurements and 6557 ADA measurements were available for the analysis.For most records, both ADA and PK were available.3. Simulations of the dependent variables from model 6 are presented in Fig. 3 and in Supplemental 4 and show that Y P K RES were well captured while Y ADA MES were not as well captured.
All models performed well with regards to the observed variable parameter estimates, which were consistent with the expectations (Supplemental 2).However, models 1 to 3 resulted in estimates of π NOADA−ADA at the lower boundary (i.e., ~ 0).Models 4 to 6 had more realistic estimates with π NOADA−ADA of 0.04, 0.02 and 0.03 for these models.In  observations and because the range of ADA measurements was large.
The mean time to first clinically positive ADA measurement (ADA measurement > 2.4 U/mL) was 81.0 days.The developed models 1 to 6 predicted the mean first time points of a transition from S NOADA to S ADA were 85, 78, 83, 64, 84 and 63 days, respectively.
Reducing model 6 to two univariate models showed that the parameter estimates in the two univariate models were, in general, similar to those obtained in the bivariate model, apart from the hidden variable parameters for the univariate ADA MES model (Table 4).With respect to individual state sequence predictions Model 6 and the univariate PK RES model performed similarly.However, the univariate ADA MES model, did not predict a transition to S ADA for six individuals although their PK indicated presence of ADA, and predicted an individual which may have a false-positive ADA record (occurring at ~ 28 h) as an ADA producer (Fig. 2).

Discussion
A common manifestation of immunogenicity is the formation of ADAs that bind to the therapeutic agent and may cause changes in the disposition of the drug and/or prevent binding to the target reducing the drug effect.ADA against biological therapeutics are routinely characterized and reported as they can be associated with treatment failure [24][25][26].The characterization and interpretation of ADA effects on disposition and efficacy of biologics can be limited by the characteristics of the bioassay used to quantify ADAs.For instance, in the presence of excess drug, the commonly used ELISA sandwich assay, may not be drug tolerant resulting in inadequate detection of ADA in the presence of high concentrations of drug [32].The presence of ADA needs to be put in the context of its effects estimated parameters (apart from ρ S ADA , RSE 110%).Models 4 and 6 performed the best with regards to the individual state sequence predictions, for model 6 predictions are shown in Fig. 2. The RSE of π ADA−NOADA was high (111.5%) in model 6, which is plausible given the low parameter estimate (0.003), indicates that transient ADA formation in the current dataset is unlikely and once an individual is in S ADA they are unlikely to transition to S NOADA .Fixing the modes of ADA in both S NOADA and S ADA to values associated with the clinical classification of ADA resulted in improved individual state sequence predictions as indicated by the difference between the predictions from model 5 (Figure S2-4), where the parameters were estimated, and model 6 (Fig. 2), where they were fixed.Simulations from model 6 (Figure S4-1 and S4-2) indicated that Y P K RES were well captured while Y ADA MES were not as well captured since the simulated percentiles did not correspond with the same observed percentiles.This is likely due to the sparse nature of ADA  ADA observations, there is no appropriate distribution as a vast majority of them were below the LOQ and thus set to the LOQ value of 0.6 U/mL.Simulations of the observed variables from the model demonstrated appropriateness of the model for the PK residuals, but simulations of ADA measurements did not perform as well, and underprediction of the number of values associated with high ADA measurements was seen.Since the two residual terms were shared for ADA MES, the very low variance of ADA in S NOADA may have contributed to the relatively low estimate of σ 2 ADA , likely responsible for this behaviour.To correlate ADA MES with PK RES, a bivariate normal distribution was assumed.Further, it was assumed that all individuals start in S NOADA , since the individuals included in the trial data were drug naïve.
Fixed effects parameter estimates for the observed variables in the tested bivariate models were as expected in most cases.In several of the models, the modes of ADA MES were fixed to clinically relevant values.A measurement that is LOQ is most likely associated with no ADA production, while a value greater than 2.4 U/mL is likely to be associated with a true ADA measurement.In our work, estimation of the modes of ADA MES in the respective states were associated with poorer individual state sequence predictions than when fixing those parameter values.Consequently, by fixing the parameters in the developed models, a weighting of the observations is made, as observations > 2.4 U/ mL are more likely to have come from S ADA than S NOADA .For instance, when µ ADA in S ADA was fixed to 1000 U/mL in a model considering only ADA measurements (univariate model), while estimating the rest of the model parameters, the model did not result in any transitions to S ADA and thus any transitions in the bivariate model would be driven solely by PK RES (results not shown).Transition probabilities in the model were generallylow, expected given the few observations associated with ADA formation.Furthermore, on clinical efficacy, safety and disposition of the drug [10].A model, such as the one presented in this work, allows us to gain additional insight into the assessment and characterization of ADA by using all sources of information.Similarly, the developed methodology was applied to ADA against CZP that does not distinguish between antibodies isotypes.In the case where specific isotypes against PEG are measured, the model could be extended to also account for different types of anti-drug antibodies.A caveat in the current analysis is that any discrepancies in PK resulting in changes in the PK residuals are assumed to be attributed to ADA formation; if the model was miss-specified, remaining unexplained PK influences would still be attributed to ADA formation.The parameters of the model could be classified as "observed variable" and "hidden variable" parameters.Two modes were estimated for each of the observed variables, one for each state together, with a single residual term shared across states.A single residual term (per variable) was used in order to reduce model flexibility and assumes that the distribution of the observations were not different depending on which state the observations arose from.If two residual terms would have been estimated for the ADA MES , the residual term related to S NOADA would have been estimated to ~ 0, as most observations would be equal to the LOQ and the other would be some large value, as measurements of ADA are highly variable.Thus, numerical issues arise as the likelihood would approach infinity due to the division of an increasingly small number in Eq. 3.This may be a consequence of the assumption of observations being equal to LLOQ, instead of recognizing that they are not measurable with the analysis instrument used.In this work, we assumed that the observed variables were normally distributed.The Gaussian distribution was deemed to be appropriate, since for well performing PK models the distribution of the residuals (which we used as the observed variable in the MHMM) are assumed to be normal.For the

Conclusion
The immunogenic response to proteins is heterogeneous and depends on different factors inherent to the drug, patient, disease and others such as dose regimen and route of administration.Characterization of ADA and its impact on PK, efficacy and safety is key in drug development of new entities.PKPD models often focus on characterising the impact of ADA on PK and hence clinical consequences.However, they do not focus on what triggers the ADA response.
We developed a model able to characterise the hidden ADA production and state transitions using all the sources of information available (both PK and ADA).To our knowledge, this is the first time that this methodology has been applied to characterize ADA, and this work has served as a basis for other posterior examples that have recently been published [33,34].The results presented in this work should serve as an example for the use of a MHMM in the analysis of immunogenicity that makes use of all available PK and ADA data, and is able to weigh in the information being more assay independent than traditional population PK analysis exploration of ADA.Although not explored here, future uses of the model may include the identification of covariates influencing the transition probabilities in the model that would help in identifying trial design and treatment related aspects that would minimize the ADA formation and thus minimize its impact on PK and/or efficacy.The methodology can be extended to other molecular entities and be used to inform on covariates that trigger immunogenicity response.Further it can be linked back to the clinical implication on safety and efficacy, thereby inform the estimate of π ADA−NOADA was in general much lower than π NOADA−ADA suggesting a lack of transient ADA formation going from S NOADA to S ADA and back again for an individual.
The estimation of IIV was associated with both issues related to parameter estimation and worse individual state sequence predictions.Although inclusion of IIV on the transition probabilities makes sense, this was poorly estimated given the data, with only very few transitions in the total data, and very few individuals having transient ADA formation.This rendered the estimation of additional parameters (such as IIV) not possible.IIV on the transition probabilities would allow covariates influencing the ADA status to be identified.This is of interest, to reduce the probability of transitioning to S ADA , which is in turn associated with a more beneficial treatment outcome.For instance, MTX administration can reduce the incidence of ADA formation, and this may be tested in the MHMM.However, this would require a larger number of transitions within an individual to be viable.The current model was estimated to a large data set.Further research is required to determine the data requirements for state-sequence predictions in smaller data sets generated in clinical practice.
When ADA are included in population PK analyses, they are often included as a dichotomous or continuous covariate on drug CL [28][29][30].These models often do not consider a time-course of ADA positivity since a subject who is ADA positive at a certain timepoint is assumed to remain positive throughout the analysis.It may be possible to use the information obtained from the present model to infer the probabilities of being in an ADA producing state and using those resulting probabilities as a covariate in the PK model.The time-course of ADA can be more readily incorporated since ADA are not the only component driving the PK response.
In the best performing model (model 6, Table 3; Fig. 2), the individual state sequence predictions obtained for a subset of individuals was able to identify potential false positive and negatives, as well as to assign ADA positive status influencing PK before the ADA was detected, and therefore before the onset of ADA effects based on the combined information coming from the PK RES and ADA MES .This highlights that the present model, by using all the information available, may fill the caveats related to the ADA bioassay that make it difficult to obtain a true characterization of immunogenicity in some instances.The modelling results may suggest that ADA are formed earlier than the assay can detect them.In individuals that tested positive for ADA, the mean time to measuring ADA was 81 days while the mean model predicted time to transitioning to the ADA producing state was 63 days.These results depend on the number of observations available, where more observations would result in more individual state sequence predictions.This methodology may be able to identify false-negative 1 3 general, for models including IIV, estimation of imprecision in parameters was unsuccessful.The estimation of ADA MES (models 1, 3 and 5) was associated with poor individual state sequence predictions as these models did not predict a transition to S ADA in individuals 5, 7, 9 and 12 despite the indication of ADA formation based on the PK data (Supplemental 3, FiguresS3-1, S3-3, S3-5).Models 2 and 5 also performed poorly for these individuals, as expected given the lower boundary estimate of π NOADA−ADA (model 2) while less so for model 5 which resulted in relatively well presented in Table

Fig. 2
Fig. 2 Individual state predictions in twelve individuals resulting from model 6 and the same model only considering univariate observed variables (either PK residuals [middle panels] or ADA measurements [lower panel]).CZP concentration measurements (green), individual CZP PK model predictions (green dashed line) and ADA measurements (red) are presented for 12 selected individuals.Each individual

Fig. 1
Fig. 1 The general model structure of the bivariate hidden-Markov model given pharmacokinetic (PK) observations (YPK RES ) and antidrug antibody measurements (YADA MES ).The two hidden states, S NOADA (grey) and S ADA (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

Fig. 3
Fig. 3 Violin plots of the observed variables (shaded areas and circles [indicate medians]) and simulations from the final model (model 6) (dashed lines and squares [indicate medians]) versus time intervals.Simulated ADA measurements below 0.6 were set to 0.6 as LOQ.

Table 1
Data used in the analysis -drug antibody; LD, loading dose; PK, pharmacokinetics; Q2W, every second week; Q4W, every fourth week a Have been included in a PKPD meta-analysis previously b Number of patients with CZP concentrations.2 subjects in CDP870-004 were excluded due to lack of information whereYPK RES and Y ADA MES are observed variables of interest, s is the current state (hidden), which can be either S NOADA or S ADA , P K RESsi and ADA MES si are state specific individual model predictions of the variables, σ 2 P K RESs and σ 2 ADA MESs 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.

•
Model (1) Modes of PK RES and ADA MES estimated.IIV estimated on modes and transition probabilities.• Model (2) Modes of PK RES estimated.ADA MES modes fixed to clinically relevant values of 0.6 and 2.4 U/mL in S NOADA and S ADA , respectively.IIV estimated on modes and transition probabilities.• Model (3) Modes of both PK RES and ADA MES estimated.IIV estimated on transition probabilities.• Model (4) Modes of PK RES estimated.ADA MES modes fixed to 0.6 and 2.4 U/mL in S NOADA and S ADA , respectively.IIV estimated on transition probabilities.• Model (5) Modes of PK RES and ADA MES estimated.No IIV estimated.• Model (6) Modes of PK RES estimated.ADA MES modes fixed to 0.6 and 2.4 U/mL in S NOADA and S ADA , respectively.No IIV estimated.
Model 6was the best performing model as indicated by individual state sequence predictions (Fig.2) and parameter estimation.Parameter estimates of the final model are 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 (YPK RES and Y ADA MES ) from the final model were compared to the observed distributions.

Table 2
Parameters in the developed bivariate hidden-Markov model

Table 3
Parameter estimates of the final bivariate model.FIX indicated parameters which were not estimated

Table 4
Parameter estimates of the bivariate model and the two univariate models handling only PK residuals or ADA measurements.Relative standard errors (RSE%) are reported when they were obtained using covariance step in NONMEM with default settings ADA measurements, and thus better relate ADA to potential clinical impact as well as limit the otherwise noise risked when classical population PK methods without characterisation of the ADA is performed; whereas identifying false-positive measurements may not be as relevant if PK is unaffected.The model performed well to correct for some false-negative ADA measurements and was able to suggest false-positives as well.In the current work, ADA measurements alone could not result in the observed state-sequence predictions.The model requires a second observed variable containing information about the impact of ADA on PK, which could then provide further information on potential clinical relevance (for example, if efficacy was impacted).Without the PK residual as a second observed variable in the current model, drawing any conclusions regarding clinical relevance would be difficult.Further work may be needed to elucidate what types of observations are most informative for clinical relevance purposes in the context of ADA assay results interpretation.