Briefly, phase Ⅰa (CTR20170306) was a multi-center, randomized, double-blinded, single ascending dose study of WBP216 in RA patients. The study protocol was designed in full compliance with Good Clinical Practice and the Declaration of Helsinki and approved by the independent ethics committee of Peking Union Medical College Hospital and Beijing hospital (Beijing, China). Subjects eligible to take part in this study were RA patients diagnosed over 6 months, treated with basic RA medications (e.g., methotrexate and leflunomide) stably for at least 28 days, aged from 18 to 70 years and with a body mass index (BMI) of 19–30 kg/m². All subjects signed the Informed Consent Form before their participation.

A total of 36 subjects were enrolled into 5 dose cohorts (10, 30, 75, 150 or 300 mg) respectively and randomized to receive WBP216 or placebo subcutaneously, with 3 subjects receiving active WBP216 and 1 subject receiving placebo in the lowest dose cohort (i.e., 10 mg) and with 6 subjects receiving WBP216 and 2 subjects receiving placebo in other each dose group. Serum samples were collected at pre-dose and 0.083, 0.5, 1, 3, 7, 14, 21, 28, 42, 56, 84, 112, 168 days post-dose to obtain total WBP216 concentration and CRP data. Total WBP216 was assayed by validated methods whose linear calibration ranges were 39.1–10,000 ng/ml for WBP216. The CRP and ESR measurements were conducted by the clinical laboratory center of the hospitals. The LLOQ for CRP and ESR are 0.01 mg/L and 1 mm/h, respectively. DAS28-ESR and ACR scores were estimated by doctors at baseline and on days 7, 28, 56, 84, 112, 140 and 168, with 112/140 days measurements only for the 75–300 mg groups. If any anti-drug antibody (ADA) measurement was positive in any subject during the trial period, this subject was described as ADA positive.

The YTE mutations in the Fc fragment of WBP216 may result in non-standard extended PK properties, which is one of our study focus. In order to explore the PK characteristics, PK parameters were firstly calculated based on individual plasma concentration-time-profiles using non-compartment analysis (NCA) by validated Phoenix WinNonlin version 8.1 software (Pharsight Corporation, Mountain View, CA, USA). PK parameters through NCA were then analyzed and explored using WinNonlin or PRISM (GraphPad 8.0.1, San Diego, CA, USA). A power model (Eq. 1) proposed by Gough et al. was applied to assess the dose proportionality (Kevin Gough, 1995). log ( Y ) = μ + β × log ( dose ) (1) in which Y denotes PK parameters such as area under curve (AUC) or maximum concentration (C_max). This approach assumes that the underlying relationship between log(Y) and log(dose) is linear. β = 1 indicates total dose proportionality. In our study, a less stringent criterion was used (Hummel et al., 2009), given the small sample size in phase Ⅰ across multiple dose groups. The estimate of β together with CI falling completely within the range of [0.5–2] was quantified as dose proportionality.

The relationship between drug exposure and response was evaluated using nonlinear mixed effects models (Phoenix NLME, version 8.1, Certara). First-order conditional estimation, extended least squares method (FOCE-ELS) was used to estimate pop PK/PD model parameters. The final structural model was determined by the objective function value (OFV) and Akaike information criterion (AIC). A sequential modeling strategy was used for fitting the models to the phase Ⅰa data. The population PK model was first developed and then the individual post hoc parameters from the final PK model were used to predict the individual WBP216 concentrations to drive the drug effect on CRP or DAS28 time profiles using appropriate PD models.

We tried one, two or three-compartment PK model with a first order, saturate or sequential zero-first order absorption compartment. The elimination phase was also analyzed by fitting to first order, saturate elimination or target-mediated-drug-disposition (TMDD) models. Based on the mechanism of action, an indirect-response model was chosen as a starting point for PD model development for CRP or DAS28. A linear model, an E_max model or a sigmoidal E_max model (Hill equation) were applied to characterize the relationship between WBP216 serum concentrations and those disease activity measures.

Inter-individual variability (IIV) was assumed to follow a log-normal distribution and was described using exponential model (Eq. 2). P i j = θ i ×   e η i j (2) where P _ij represents the individual value of the parameter for the ith parameter in the jth individual, θ _i depicts the population typical value for the ith parameter, and ηij represents random effect in jth individual sampled from a normal distribution with a mean of zero and variance of ω ².

The residual unexplained variability of WBP216 concentration and PD observation was described by a proportional or additive error model, respectively (Eqs 3, 4). Y obs   =   Y pred × ( 1 + ε 1 ) (3) Y obs = Y pred + ε 2 (4) where Y _obs and Y _pred are the observed and predicted serum concentration in plasma or disease activity measures. And ε ₁ is the proportional and ε ₂ is the additive component of the residual error model, respectively. Both of ε ₁ and ε ₂ are assumed to be normally distributed in the range of (0, σ ²).

After collinearity diagnostics and correlation analysis, possible covariates including weight, age, sex, baseline serum albumin (ALB), alanine amiotransferase (ALT), creatinine clearance (CLcr), total bilirubin (TBIL), free IL-6, ADA (negative or positive) etc. were tested on both PK and PD parameters. Continuous covariates were described using the power function and categorical covariate were modeled by exponential function, see Eqs 5 and 6. Effect i = ( Cov i j / Cov median ) θ cov i (5) Effect i = e Cov i j · θ cov i (6) where Effect _i is the multiplicative factor for covariate i, Cov _ij is the covariate value for individual j, Cov _median is the median covariate value, and θcov _i is the exponent or parameter for covariate i model.

Potential covariates were incorporated into the base model one by one using stepwise forward inclusion. When OFV decreased by 6.63 (at p < 0.01), the covariate was selected for inclusion to develop a full model, followed by the backward elimination. The covariates were subtracted one at a time in a stepwise manner as well once OFV increased above 10.83 (at p < 0.001, df = 1), until all remaining covariates in model were statistically significant.

During the process of models building, the goodness of fit (GOF) of different models was compared on the basis of OFV and AIC. GOF was graphically evaluated by inspecting plots of the individual or population predicted vs. observed values, and plots of the conditionally weighted residuals (CWRES) vs. population predictions or time.

Models were also validated internally using prediction-corrected visual predictive checks (pcVPC) as well (Bergstrand et al., 2011). On the basis of 1,000 times pcVPC simulation, the 90% prediction interval (PI) was compared with the 90% interval of the prediction-corrected observations. Bootstrap analysis was also performed for the final model along with a total of 500 data sets resampling randomly from the original data set (Ette, 1997). We reported the calculated 90% confidence interval (CI) of model parameters from successfully minimized runs.

Simulation was conducted in Phoenix NLME (version 8.1, Certara, USA) based on a Monte-Carlo simulation approach. Up to 100 Phase Ib trials were simulated using the uncertainty distribution in parameter estimates. 27 patients in each simulated trial were simulated using the IIV log-normal distribution in both PK and PD parameters. The distribution of covariates still leveraged those of phase Ⅰa data set. Serum CRP levels and DAS28 were simulated for a range of maintenance doses (30–300 mg) under three different dosing frequencies: every 4 weeks (Q4W), every 8 weeks (Q8W), and every 12 weeks (Q12W). The duration of drug effect was simulated up to week 72 with weekly virtual PD sampling. It was assumed that the PK/PD relationship based on the Phase Ⅰ study lasting 24weeks could be extrapolated to longer term studies.

Tocilizumab has proven its successful clinical efficacy, so it was used as the reference for comparison of clinical endpoints. Since mean CRP is decreased by around 90% and mean DAS28 could be reduced by over 56.5% using the recommended dosage of tocilizumab (ACTEMRA^® HIGHLIGHTS OF PRESCRIBING INFORMATION; U.S. Food and Drug Administration), ∆CRP≥90% and ∆DAS28 ≥ 56.5% were set as our target efficacy to optimize WBP216 phase Ⅰb maintenance dose levels and dosing frequencies.

WBP216 would take a long time to achieve steady state exposure and efficacy because of its slow elimination rate constants. Because RA patients would require a rapid pain relief, a loading regimen of WBP216 would be necessary. Hence, we simulated two categories of loading regimen to achieve steady state exposure: Firstly, WBP216 was given in a more intensive frequency at an initial three-administrations, including 0–4–8th week, 0–2–4th week or 0–2–6th week; The second simulation used a loading dose, that doubled the confirmed maintenance dose level.

Up to 36 RA patients took active medicine and placebo in a ratio of 3:1, respectively, wherein 27 patients received WBP216. One of the subjects in the 75 mg group showed a huge fluctuation of CRP level after 21 days from administration, very different from other subjects. The CWRES of there CRP samples were also greater than six during CRP model exploration. This subject was considered as an outlier and excluded from this PK/PD model analysis. Descriptive statistics (baseline values) of potential covariates of 26 patients tested in the PK/PD analysis were summarized in Table 1. A total of 391 PK samples were obtained during the phase Ⅰa and PD data consisted of 384 CRP samples and 241 DAS28 samples. Since fewer than 10% of PK samples were below the lower limit of quantification (LLOQ), they were handled by M1 method (Keizer et al., 2015). PK concentrations whose corresponding ADA was positive were all above the LLOQ, so they were not discarded.

We performed an NCA analysis for different dose groups before population PK model development to understand the PK properties fully because of the unique YTE mutations in WBP216. Figure 1 showed that mean apparent clearance (CL/F) tended to increase with increasing dose, contrary to the clearance change pattern of TMDD, which generally has a high clearance in lower dose groups (Mager, 2006). Large variability of CL/F among individuals was observed in both 75 and 300 mg dose groups. High individual variability of apparent volume (V/F) was also observed in higher dose levels (75–300 mg). Like CL/F, V/F presented an increasing trend over dose levels. The phenomenon of CL/F and V/F changing with dose levels was speculated to result from the changes of either the actual increased CL and V or decreased bioavailability (F). To clarify the real reason, the half-life of WBP216 was analyzed. Half-life (t _1/2 = 0.693*V/CL) is considered not to be impacted by F. In this analysis of half-life, the value of t _1/2 distributed evenly in five dose groups, and the mean t _1/2 remained almost the same in different groups, except for the 150 mg group, likely caused by lower CL/F and higher V/F value in this group. Those changing trends of CL/F, V/F and t _1/2 suggested that the decrease of F with doses was probably the cause of the increase of CL/F and V/F. Mean t _1/2 is 40–60 days, indicating the potential for a long dosing interval in therapeutic use. The power model fitted ln (AUC) over ln (dose) data well (Figure 1D). Parameter β was calculated as 0.652 with 90% confidence interval (CI), [0.43, 0.87], part of which fell outside the lower limit of prespecified range [0.5, 2] (Hummel et al., 2009). This dose-dependent study indicated little lower than proportional increases in exposure (AUC).

According to AIC value and goodness of fit (GOF) plots, the selected final proposed PK/PD model structure is shown in Figure 2. The WBP216 serum concentrations were best described by a two-compartment PK model with sequential zero-first order absorption and first order elimination (see Eqs 7–9). An indirect-response model with an drug E_max inhibition of the CRP or DAS28 zero order rate production constant (K_in) best described the disease measures (see Eqs 10, 11). Residual variability was characterized by a proportional error model for serum concentration, CRP and DAS28. All PK/PD model parameters are summarized in Table 2. d A depot d t = Dose T d − ka × A depot ( t = 0 ,   A depot = 0 ) (7) d A 1 d t = k a × A depot − C L / F × C 1 − Q / F × ( C 1 − C 2 )   ( t = 0 , A 1 = 0 ) (8) d A 2 d t = Q / F × ( C 1 − C 2 )   ( t = 0 ,   A 2 = 0 ) (9) d CRP d t = K in,CRP × ( 1 − E max , C R P × C 1 E C 50 , C R P + C 1 ) − K out,CRP × CRP (10) dDAS 28 d t = K in , DAS 28 × ( 1 − E max , D A S 28 × C 1 E C 50 , D A S 28 + C 1 ) − K out,DAS 28 × DAS 28 (11) Wherein A _depot, A1, A2 representing WBP216 amounts in absorption depot, systemic central and peripheral compartments respectively, were equal to zero when t = 0. C1 and C2 were the concentration in central and peripheral compartment, respectively and equal to zero before dosing. Response was CRP or DAS28, and was equal to its baseline value before drug administration.

The final PK model fitted the observed concentration data better by assuming that after WBP216 was subcutaneously administered, drug was released from the injection site at a zero order rate to a depot compartment, and was then absorbed to a central compartment. Td is the period of zero order release and Ka denotes the first order absorption rate. The reason of selecting such a more complex absorption model will be discussed later. Based on the mechanism of action, WBP216 neutralizes IL-6 and inhibits CRP production, and further slow down disease progression, so it is plausible to set the drug inhibition on response production rate, K_in (Sharma and Jusko, 1998). As shown in Table 2, the relative standard errors (RSE%) for almost all fixed-effect parameters were ≤30.9%. The uncertainties for random-effect parameters were <27%. Overall, the precision of parameters estimates was acceptable. Parameters showed various inter-individual variability, ranging from 17 to 128.4%. E _max for CRP was set to 1 since its value was always very close to 1 in all tested runs. E_max for DAS28 was also fixed to 1, thus leading to a straightforward convergence of the model. All applicable shrinkage was below 18.1%, which was smaller than reported cut-off value 30% and could assure accurate IIV estimates and avoidance of misleading diagnostic plots (Savic and Karlsson, 2009).

The correlation diagnosis chart of various covariates was presented in Supplementary Figure S1. The potential impacts of demographics and laboratory data baseline on PK/PD parameters of WBP216 were tested using a stepwise covariate modeling procedure. Those statistical significant covariate effects were identified and retained in the final model: ALT on CL, baseline free IL-6 on EC_50,CRP and K_in,CRP (see Eqs 12–14 and Supplementary Table S2). Continuous covariates were described using the power function, centered by the median value. Apparent clearance decreased with increasing ALT with exponent −0.833, which explained around 18% of CL/F inter-individual variability (IIV). Higher baseline free IL-6 levels could result in increased K_in,CRP (exponent 0.695), while lead to its EC_50,CRP reduction (exponent −0.772). Baseline free IL-6 accounted for 59.2% IIV of K_in,CRP and 55.6% IIV of EC_50,CRP. No covariates were found statistically significant in parameters for DAS28. Covariate associated parameters are also presented in Table 2. CL/ F = θ CL/ F × ( ALT 11 ) − 0.833 × e η CL/F (12) K in,CRP = θ Kin,CRP × ( BaseFreeIL-6 24.9 ) 0.695 × e η K in,CRP (13) EC 50 ,CRP = θ EC 50 ,CRP × ( BaseFreeIL-6 24.9 ) − 0.772 × e η EC50,CRP (14)

Goodness-of-fit (GOF) plots for the final PK model in serum are shown in Figure 3. Plots of the population- and individual- predicted concentration vs. observed concentrations demonstrate no major bias. The conditionally weighted residuals (CWRES) were symmetrically distributed about zero axis and most points laid within the acceptable range (−2 to 2), suggesting that little to no bias accompanied with concentration or time. The GOF plots for CRP and DAS28 model were presented in Supplementary Figures S2 and S3, respectively, which also performed well in visual diagnostic.

The predictive performance was evaluated internally by pcVPC. Plots of pcVPC were presented in Figure 4. We can see that the 5th, 50th and 95th percentiles of prediction-corrected observations and predicted data were fairly consistent, especially a better match for DAS28. The 95th percentile of predicted PK data through VPC is a slight under-prediction and the CRP model over-predicted drug inhibition at 5th percentile slightly, which will be discussed in the part of discussion. Despite these small deviations, the 90% prediction interval of simulated data covered most of the observations. The bootstrap results are also shown in Table 2. The typical values of parameters and IIV estimates in the final model were pretty close to the median values from bootstrap validation, and fell within 90%CI of bootstrap parameters completely, which indicated high stability and precision of the final model.

Overall, the good performance of GOF, pcVPC plots and bootstrap estimations reconfirmed that the final PK/PD model was adequately developed and the predictive performance was sufficient to capture PK/PD observations.

Based on the observations, we chose dose levels of 30, 75 and 150 mg as potential maintenance doses so as to rapidly reach the target effect (mean ∆CRP≥90%; ∆DAS28 ≥ 56.5%, i.e. ∆DAS28 ≥ 3 when average baseline DAS28 is 5.3) at steady state, equivalent to tocilizumab (ACTEMRA^® HIGHLIGHTS OF PRESCRIBING INFORMATION; U.S. Food and Drug Administration). Then, simulations of those different maintenance dose levels and three varying dose frequencies Q4W, Q8W and Q12W were performed, respectively. Figure 5 showed the Monte-Carlo simulations of CRP and DAS28 changing from baseline over time under different dose regimens using the final PK/PD model. Consistent with observations, the simulation results showed that CRP would decrease very rapidly to the nadir within the first week, while DAS28 changed slowly and reached a plateau after 24 weeks. Multiple dose regimens could decrease the mean CRP levels more than 90% from baseline and reduce DAS28 by more than 3 units, except for 30 mg Q8W, 30 mg Q12W, 75 mg Q12W scenario. The optimal regimens, in terms of achieving mean ∆CRP≥90% and ∆DAS28 ≥ 3, seemed to be 75 and 150 mg as maintenance doses administered every 8 weeks. Approximately 81% and 92% of virtual patients achieved DAS28-ESR <2.6 (the cutoff of RA remission defined by EULAR), when dosed at 75 mg Q8W and 150 mg Q8W, respectively. The 30 mg dose level did not achieve target efficacy when dosed every 8 or 12 weeks. Dosing every 4 weeks was a bit over-dosing and inconvenient to use in clinical practice, however, Q12W dosing frequency resulted in larger fluctuation in efficacy, especially in CRP, and led to a lower proportion of patients nearby target reference line no matter what dose level.

EULAR (European league against rheumatism) response criteria state (Broeder et al., 2002; Wells et al., 2009) (see Supplementary Table S1): moderate responders were patients with an improvement of >0.6 and a present DAS28 score of <5.1, or an improvement of >1.2 and a present DAS28 score of >5.1. We predefined that once patients start having a moderate response to WBP216, they would begin to feel pain relief intuitively. Patients’ time needed to reach at least moderate response are summarized in Figure 6. The results showed that over half of patients got pain relief within around 10 days, while 12/26 patients did not reach moderate response until at least 28 days from drug administration. WBP216 caused a gradual increase in clinical efficacy during the phase Ⅰa period, which was consistent with the patients’ reports of slow onset of drug action.

Therefore, in order to allow patients to benefit from WBP216 therapy more quickly, loading dose regimens were designed. On the basis that the optimal dosing frequency was Q8W during maintenance therapy, we simulated four kinds of “loading dose schedules”, Figure 7 showed examples of varying loading schedules followed by a maintenance dose level of 150 mg Q8W: (1) Loading dose at weeks 0, 4 and 8 followed by Q8W still took more than 8 weeks to reach steady state exposure. (2) Loading dose at weeks 0, 2 and 4 followed by Q8W caused a steep WBP216 concentration increase, exceeding the peak concentration achieved by the maximum dose level (300 mg) in phase Ⅰa, which may raise potential safety concerns. (3) Dosing at weeks 0, 2 and 6 followed by Q8W would reach plateau exposure at the second drug administration. (4) A loading dose of 300 mg, doubling the maintenance dose level, resulted in steady state instantly at the first dosing. Hence, the last two loading dose schedules are highly recommended for future clinical trials of WBP216. The simulated DAS28 profiles under the optimized loading dose regimen (3) and (4) were presented in Supplementary Figure S4, showing that the target efficacy was reached three weeks in advance.