Study design

Model development was performed using data from a Phase 2a, randomized, double-blind, placebo-controlled, parallel group, multicenter study in subjects with moderate-to-severe plaque psoriasis (ClinicalTrials.gov Identifier: NCT02969018). A schematic of study design is shown in Fig. 1. The first part of the study, following a screening period (up to 6 weeks), was a 4-week induction period with double-blind daily treatment. At the end of week 4, all subjects switched to their predefined double-blind maintenance treatment regimen for week 5 through week 12.

Approximately 200 subjects were planned to be randomized into the study, to allow for approximately 160 evaluable subjects (20 completers per arm). The randomization ratio was 7:1, active: placebo. During the first 4 weeks of the treatment period, subjects received orally either 30 mg once-daily (QD), or 60 mg QD of brepocitinib, or matching placebo. During the 8-week maintenance portion of the treatment period (weeks 5 through 12), subjects received orally either 10 mg QD, or 30 mg QD, or a 100 mg once weekly (QW) regimen of brepocitinib, or matching placebo. Maintenance dose level and regimen were assigned at the initial time of randomization into the study. All subjects, regardless of assigned regimen (i.e., QD or QW) received blinded QD tablets throughout the study treatment period to maintain the study blind.

The duration of study subject participation was approximately 26 weeks, including screening (up to 6 weeks), 12-week treatment period, and 8-week follow up period. For further details on the study design, the reader is referred to the associated publication [9].

Female and male patients between 18 and 75 years of age with a diagnosis of plaque psoriasis for at least 6 months before the start of the study, PASI score of ≥ 12, physician’s global assessment (PGA) score of 3 or 4, and psoriasis covering ≥ 10% of total body surface area were eligible to participate in the study. Key exclusion criteria were non-plaque psoriasis, other skin conditions that would affect the assessment of psoriasis, drug-induced psoriasis, use of corticosteroids, and psychiatric conditions including suicidal ideation or behavior.

The study was conducted in compliance with the Declaration of Helsinki and Good Clinical Practice Guidelines established by the International Council on Harmonisation. The final protocol, amendments, and informed consent documentation were reviewed and approved by the institutional review boards and independent ethics committees of the investigational centers.

Study assessments

PASI scoring, which quantifies body surface area and lesion severity into a single score, was used to assess brepocitinib efficacy in this study. PASI score is calculated by combining the percentage of body areas (head and neck, upper limbs, trunk, and lower limbs) covered, with the severity of erythema, thickness/induration, and desquamation/scaling, for a score between 0 and 72 [10]. PASI scores were assessed: at screening; on week 0 (baseline); on weeks 1, 2, 4, 6, 8, 10, and 12; and during follow-up (weeks 14 and 16). Only data collected during the 12-week active treatment period of the study were included in the analysis. Any data collected during screening and the follow-up period were not included in the analysis.

Blood samples for PK analysis of brepocitinib were collected pre-dose on week 0 (baseline) and week 1; pre-dose and 30 min post-dose on weeks 2, 6, 8, and 10; and pre-dose, 30 min, 1, 2, and 4 h post-dose on weeks 4 and 12. The samples were analyzed using high-performance liquid chromatography-tandem mass spectrometry. The lower limit of quantification for brepocitinib was 0.2 ng/mL.

For additional details and study assessments not related to the current analysis, the reader is referred to the associated publication [9].

Population PK model and derivation of Cave for exposure-response modeling

A population pharmacokinetic (PK) model has been previously developed for brepocitinib using data from five clinical trials, consisting of three Phase 1 and two Phase 2 studies (including the Phase 2a study in psoriasis patients described in this manuscript) [11]. Briefly, brepocitinib PK were described with a one-compartment model with first-order oral absorption and an absorption lag for the tablet formulation (apparent clearance (CL/F) of 18.7 L/h, apparent volume of distribution (V/F) of 136 L, first-order absorption rate constant (ka) of 3.46 h−1 and a lag time (Alag) of 0.24 h). The effect of body weight on CL/F and V/F was included with an allometric relationship (referenced to 70 kg) and the associated coefficients were fixed to 0.75 and 1, respectively. Random inter-individual variability was accounted in the CL/F and V/F parameters using a full variance-covariance matrix (coefficient of variation (CV) of 78% and 60.5% for CL/F and V/F, respectively, and a correlation coefficient of 0.76).

The empirical Bayes estimates (EBEs) of CL/F for all study participants (η-shrinkage < 0.1%) were extracted from the population PK model output and were used to calculate the average concentration (\(_\)) of each individual across the induction and maintenance periods of the study, using Eq. 1, where \(\tau\) refers to the dosing interval (24 h for QD dosing and 168 h for QW dosing) and \(Dose\) refers to the administered brepocitinib dose either during the induction or the maintenance period of the study.

The derived \(_\) for each study participant was the exposure-relevant metric that was used in the exposure-response analysis (see Population exposure-response model section, Eq. 12). Brepocitinib PK is assumed to be at steady state during all efficacy assessments performed over the active treatment period of the study (since brepocitinib has a short terminal half-life, ranging from 3.8 to 7.5 h, and efficacy assessments were performed at least 1 week since first dose either during the induction or maintenance period).

Population exposure-response model

An approach that employs a zero-inflated beta distribution [12, 13] was used to allow modeling of PASI scores without disregarding their bounded nature (range from 0 to 72). To enable this, observed PASI scores were first transformed to the 0 to 1 interval using Eq. 2.

Note that the analyzed dataset although it contains observations on the lower bound (PASI = 0), it does not contain any observations on the upper bound (PASI = 72), see Fig. 2, thus \(y\in \left[0,\right.\left.1\right)\). In the case that the dataset had contained data on both boundaries, the use of a zero- and one-inflated beta distribution [12, 13] would have been appropriate.

The transformed PASI scores (\(y\)) were then assumed to follow a zero-inflated beta distribution with a corresponding probability density function that is given by the mixture in Eq. 3,

$$p\left(y;\, _,\,\alpha ,\,\beta \right)=\left\_ ,i\!f \,\,y=0\\ \left(1-_\right)\cdot f\left(y;\alpha ,\beta \right) , i\!f \,\,0<y<1\end\right.$$

(3)

where \(_\) is the probability of a 0 observation (i.e., PASI = 0) and \(f\left(y;\alpha ,\beta \right)\) is the density function of the beta distribution defined in Eq. 4,

$$f\left(y;\alpha ,\beta \right)=\frac\cdot ^\cdot ^ $$

(4)

where \(\Gamma\) (∎) is the gamma function and \(\alpha\) and \(\beta\) are shape parameters of the beta distribution (with \(\alpha\) > 0 and \(\beta\) > 0). The Nemes approximation to the gamma function was used [14, 15, 16], as defined in Eq. 5, where \(\) represents \(\alpha\), \(\beta\), or \(\alpha +\beta\) in the density function above (note that the GAMLN function can alternatively be used in NONMEM® Version 7.3 [17] onwards).

$$\varGamma \left(\right) \sim }\right)}^\cdot \sqrt}\cdot ^}\right)}^X}$$

(5)

The shape parameters \(\alpha\) and \(\beta\) were parameterized (see Eqs. 6 and 7) with respect to the expected value of the beta distribution \(\mu\) and the precision parameter \(\phi\).

$$\alpha =\mu \cdot \phi $$

(6)

$$\beta =\left(1-\mu \right)\cdot \phi$$

(7)

Under this parameterization, the variance (\(^\)) of the beta distribution is defined in Eq. 8. Thus, the parameter \(\phi\) plays the role of a precision parameter in the sense that for a given value of \(\mu\), variance decreases as the value of \(\phi\) increases.

The precision parameter \(\phi\) was estimated as a fixed effect parameter, while \(\mu\) was parameterized (see Eq. 9) with respect to baseline PASI score (\(BSL\)), the placebo effect (\(_\left(t\right)\)) and the drug effect (\(_\left(t\right)\)).

$$\mu =\frac\cdot \left(BSL-_\left(t\right)-_\left(t\right)\right) $$

(9)

The placebo effect (\(_\left(t\right)\)) was modeled as illustrated in Eq. 10 with an empirical relationship, where \(_\) represents the maximum placebo effect (expressed as fraction of baseline), and \(_\) represents the rate of onset of the placebo effect.

$$_\left(t\right)=BSL\cdot _\cdot \left(1-^_\cdot t}\right) $$

(10)

The drug effect (\(_\left(t\right)\)) was modeled as illustrated in Eq. 11 through a latent variable \(R\) on which the drug elicits its effect via a Type I indirect response model (Eq. 12), similarly to [18].

$$_\left(t\right)=BSL\cdot \left(1-_\right)\cdot \left(1-R\left(t\right)\right) $$

(11)

$$\frac=_\cdot \left(1-\frac_}_+_}\right)-_\cdot R\left(t\right)$$

(12)

It was further assumed that at baseline the latent variable \(R\) takes the value of 1 (i.e., \(R\left(0\right)=1\)) and consequently \(_=_\) (i.e., equal rates of onset and offset of the drug effect).

The model parameterization described in Eq. 9, 10, 11, 12 assures that \(\mu\), as the expected value of the beta distribution, is constrained within the \(\left(\text\right)\) interval. The value of \(\mu\) can be interpreted as the individual PASI score prediction after transformation to the \(\left(\text\right)\) domain (Eq. 9).

The probability of a 0 observation (\(_\)), see Eq. 3, was modeled through Eq. 13 as a function of \(\mu\) and two additional parameters \(_\) and \(_\) (with \(_\) > 0). The rationale for this parameterization is to enforce through a flexible function that the smaller the value of \(\mu\) for a given subject, the higher the probability of an actual PASI=0 observation.

$$_=\frac^_-_\cdot \mu \right)}}^_-_\cdot \mu \right)}} $$

(13)

Incorporation of random effects on the parameters of the exposure-response model was assessed using different parameterizations and covariance structures in order to take into account inter-individual variability in the observed response. Since the \(_\) and \(BSL\) parameters need to be constrained at the individual level within a bounded region, inter-individual variability on these parameters was assessed using a generalization of the logit-normal distribution (see Eq. 14) [19],

$$_=_+\left(_-_\right)\cdot \frac^_-_}_-_}\right)+_}}^_-_}_-_}\right)+_}} $$

(14)

where \(_\) is the individual value of parameter \(P\), \(_\) is the typical population value of parameter \(P\), \(_\) is the random effect term with respect to the inter-individual variability in parameter \(P\) assuming to follow a normal distribution with mean of 0 and variance \(^\), and \(_\) and \(_\) are the lower and upper bounds respectively of parameter \(P\). \(_\) and \(_\) were set to 0 and 1 respectively for \(_\) (since it represents a fraction) and to 12 and 72 for \(BSL\) (since baseline PASI score cannot be less than 12 per study’s inclusion criteria).

No covariates were tested in the current exposure-response analysis.

Assessment of model performance

In addition to evaluating the relative standard errors derived from the NONMEM® covariance step, the parameter uncertainty of the final model was also assessed using sampling importance resampling (SIR) [20, 21]. SIR was performed by subjecting the final model covariance step output to five iterations of sampling, resampling, and multivariate Box-Cox transformation. Samples for the sampling and resampling steps used 1000, 1000, 1500, 2000, 2000 and 200, 250, 500, 1000, 1000 samples, respectively. The median and 95% confidence intervals (CIs) for each parameter were calculated from the SIR resamples of the final iteration.

The model’s performance to adequately describe the observed PASI scores was assessed via a visual predictive check (VPC) based on 1000 simulations using the design of the index dataset. Additionally, although the model was developed using the raw PASI scores, it is of particular importance to be also able to describe derived responder metrics, as the latter are important endpoints for future studies. Therefore, the model performance was also evaluated with respect to its capacity to predict the following responder metrics: proportion of participants achieving at least 50%, 75%, 90%, and 100% improvement from baseline PASI (PASI50, PASI75, PASI90, and PASI100, respectively). For this evaluation these responder metrics were derived from the raw simulated PASI scores in each of the 1000 simulated datasets, and associated 95% CIs were generated. The latter were then compared to the responder metrics that were observed in the current study.

Clinical trial simulations for dose selection in future studies

The developed population exposure-response model was used to perform clinical trial simulations to aid dose selection for subsequent trials. Variability in brepocitinib exposure was accounted by sampling \(CL/F\) from the associated population distribution determined in the population PK model after incorporating the effect of weight (weight of each simulated individual was randomly sampled with replacement from the current study dataset). The individual \(CL/F\) values were subsequently used to calculate \(_\) for each simulated individual using Eq. 1. The output \(_\) values were then passed to the developed population exposure-response model to simulate PASI score response in clinical trials.

A set of 2000 clinical trials were simulated assuming that brepocitinib or placebo were administered QD for 16 weeks. Each participant was assumed to be assigned the same treatment throughout a 16-week period (no induction/maintenance periods). Brepocitinib dose levels between 5 and 60 mg were evaluated in 5 mg intervals. Each treatment arm (either brepocitinib or placebo) was assumed to have a sample size of 56 participants. Efficacy assessments (PASI scores) were assumed to be performed at baseline and every 2 weeks until week 16.

The PASI75, PASI90, and PASI100 responder metrics after 16 weeks of treatment are expected to be the endpoints of interest in subsequent studies. Thus, the output from the clinical trial simulations (raw PASI scores) was summarized to derive predictions for these responder metrics on week 16 together with associated 90% CIs across all the evaluated dose levels.

Using the exact framework described above, additional exploratory clinical trial simulations were performed using selected dose levels in order to compare the projected efficacy of flat dosing (same dose throughout the 16-week period) with that of an induction/maintenance dosing paradigm (4 weeks of induction dose followed by 12 weeks of maintenance dose).

Modeling software

The population exposure-response model was developed using non-linear mixed effects methods and NONMEM® Version 7.3 [17]. Population parameter estimation was performed with the Laplacian estimation method. The ADVAN13 subroutine with TOL = 7 was used for solving differential equations. NONMEM® control stream is provided in Online Resource 1.

SIR was conducted using Perl speaks NONMEM® (PsN) (version 5.2.6) [22, 21]. Data visualization, exploratory analyses, model diagnostics, post-processing of NONMEM® output and all clinical trial simulations were generated using the R statistical and programming language [23] (version 3.6.1).