algorithm can be useful to explore all the design space, but could result in time
points that have no clinical relevance. The Fedorov-Wynn type algorithm is
one of the earliest and most notable algorithms for finding optimal designs in
a discrete space and can be proved to converge to the optimal design for con-
vex criterion functions. However, it does not have good computation speed in
complex design problems [14].

With the increasing use of modeling techniques when analyzing longitudinal
data, there is a need to explore new optimization algorithms to design these
studies appropriately. Multiplicative algorithms optimize in discrete design
space and their convergence for D-optimality and other criteria has been gen-
erally established [15]. Several advantages were highlighted with multiplicative
algorithms, in particular a better computation speed than the Fedorov-Wynn
algorithm [14]. So far, multiplicative algorithms have been mostly explored in
standard individual regression. That is why we aim to explore these algorithms
further in mixed effect models and to illustrate it with different applications.

In this article, we present the principles of a multiplicative algorithm [15] and
how it can be applied to optimize population designs in clinical trials analyzed
by NLMEMs. The presented optimization method is based on expression of
the FIM. We illustrate these methods through two longitudinal data examples.
The first example is inspired by a longitudinal study of a monoclonal antibody
indicated in the treatment of wet age-related macular degeneration (wet-AMD)
[16]. Considering a mixed effect model previously developed for a continuous
outcome, we aim to optimize both the measurement times and doses in patients.
The second example is a trial with repeated binary outcomes in two treatment
groups, modeled by a mixed logistic model, inspired by a previous study [17, 18].
We aim to optimize the measurement times considering either identical design
for the two treatment arms.

In section 2, we introduce the notations for NLMEMs and FIM. In section
3, we present the multiplicative algorithm that we used in this work. Sections
4 show an application of the algorithm. Finally, we discuss the results and
perspectives of this work in section 5.
## 2. Population Fisher Information Matrix

### 2.1. *Design*

Let $N$ be the number of subjects and $\xi_i$ the elementary design of individual $i$, defined by the number $n_i$ of samples and their allocation in time ($t_{i1}, \dots, t_{in_i}$), and, in some cases, additional design variables (like dose $d_i$). The design space for $\xi_i$ is discrete i.e. $\xi_i$ is included in the set of $\xi_q$ with $q = 1, \dots, Q$. The $Q$ possible elementary designs are predefined according to the constraints that are given in terms of dose and time measurements (number and allocation). Usually, the population design $\Xi$ is defined by the total number of subjects $N$ and the $Q$ elementary designs with weight $w_q$ associated with each design $\xi_q$ as: $\Xi_w = \begin{pmatrix} \xi_1 & \dots & \xi_Q \\ w_1 & \dots & w_Q \end{pmatrix}$, so that $\Xi = (N, \Xi_w)$, with the weight set $w = (w_1, \dots, w_Q)$. Each weight represents the proportion of the total number of subjects $N$ to be taken for each elementary design $\xi_q$ (the number of subjects allocated to one elementary design is $N_q = N \times w_q$). Each $w_q$ satisfies the equation $\sum_{q'=1}^Q w_{q'} = 1$, $0 \le w_q \le 1$, $\forall q = 1, \dots, Q$. The algorithm tries to find the optimal weight set $w$, i.e. the optimal proportions of subjects allocated to each possible elementary design $\xi_q$.

In practice, after optimization of the population design (i.e. after optimization of the $w_q$ by the multiplicative algorithm), designs with a very low weight (i.e. $N \times w_q = 0$ after rounding to the closest integer) were discarded. Then, weights corresponding to the $Q'$ remaining designs are summed up and each remaining weight is divided by the sum to obtain new weights $w_{q'}$ whose sum is equal to 1. The number of subjects per group $N_{q'}$ is rounded to the closest integer of $N \times w_{q'}$, as done when using the Fedorov-Wynn algorithm for nonlinear mixed models [13], with $N_1, \dots, N_{Q'}$ satisfying the equation $\sum_{q'=1}^{Q'} N_{q'} = N$.
## 2.2. Nonlinear mixed effect models for continuous data

The NLMEM to describe observations $y_i$ (of length $n_i$) of the $i^{th}$ individual, with its design $\xi_i$ included in the set of $(\xi_q)_{q=1,...,Q}$ has the following form:

$$y_i = f(\theta_i, \xi_i) + \varepsilon_i \quad (1)$$

where the function $f$ provides a $n_i$-vector of model predicted values for individual $i$ and its elementary design $\xi_i$. An example of $f$ is given in section 4.2 (with $f$ a disease progression model for visual acuity). $\theta_i$ is the vector of the parameters for individual $i$, expressed as a function $g(\mu, b_i, z_i, \beta)$ of: the vector of fixed effects $\mu$, the vector of random effects $b_i$ and the vector of covariates $z_i$ for individual $i$, as well as the vector of covariate effects $\beta$. Random effects $b_i$ follow a normal distribution with mean zero and variance-covariance matrix $\Omega$ which accounts for the inter-individual variability, i.e. $b_i \sim N(0, \Omega)$. $\Omega_u$ is the vector containing all unique elements of $\Omega$. Assuming a simple model without covariate effects, $g$ is usually additive, so that $\theta_i = \mu + b_i$, or exponential so that $\theta_i = \mu \times e^{b_i}$. $\varepsilon_i = (\varepsilon_{i1}, ..., \varepsilon_{in_i})^T$ is the vector of random errors for individual $i$ and $n_i$ is the number of its observations. $\varepsilon_i \sim N(0, \Sigma_i)$, with $\Sigma_i$ a diagonal matrix of dimension $n_i$ (assuming independent $\varepsilon_i$ for each subject) defined as:

$$\Sigma_i = \mathrm{diag}(\sigma_{\mathrm{inter}} + \sigma_{\mathrm{slope}} f(\theta_i, \xi_i))^2 \quad (2)$$

where $\sigma_{\mathrm{inter}}$ and $\sigma_{\mathrm{slope}}$ are standard deviations of additive (i.e. the constant part) and proportional terms of the error model, respectively.

The vector $\Psi^T$ of length $P$ contains all population parameters to be estimated ($\Psi^T = (\mu^T, \beta^T, \Omega_u^T, \sigma_{\mathrm{inter}}, \sigma_{\mathrm{slope}})$).

The notations of NLMEMs for discrete data are available in Appendix A.

## 2.3. Population Fisher information matrix

Assuming independence across $N$ individuals, the population FIM, $M_{PF}(\Psi, \Xi_w)$, for a population design $\Xi_w$, is defined by the weighted sum of all the elementary matrices $M_{PF}(\Psi, \xi_q)$,

$$M_{PF}(\Psi, \Xi_w) = \sum_{q=1}^{Q} w_q M_{PF}(\Psi, \xi_q) \quad (3)$$
For individual $i$, the elementary FIM, $M_{PF}(\Psi, \xi_i)$, for the elementary design $\xi_i$, is defined as

$$M_{PF}(\Psi, \xi_i) = E_{y_i, z_i} \left( \frac{\partial \log L(\Psi, y_i, z_i)}{\partial \Psi} \frac{\partial \log L(\Psi, y_i, z_i)^T}{\partial \Psi} \right) \quad (4)$$

where $\log L(\Psi, y_i, z_i)$ is the log-likelihood of the vector of observations of individual $i$ for parameter $\Psi$, given by

$$\log L(\Psi, y_i, z_i) = \log \left( \int_{b_i} p(y_i|b_i, z_i, \Psi) p(b_i|z_i, \Psi) db_i \right) \quad (5)$$

with $p(y_i|b_i, z_i, \Psi)$ the probability density function (p.d.f.) of $y_i$ given parameters $\Psi$, random effects $b_i$ and covariates $z_i$, and with $p(b_i|z_i, \Psi)$ the p.d.f. of $b_i$ given $\Psi$ and $z_i$.

As there is no analytic expression for the likelihood in NLMEMs, approximations such as linearization can be used to approximate the Fisher information matrix [19, 20]. For continuous outcome models, when linearizing the model around a guess value of the fixed effects, the variance of the observations is then independent of the fixed effects, which leads to a block-diagonal expression of the FIM [20]. It has been shown that this expression of the FIM performed better than the full matrix expression given by other linearization approaches, providing results closer to those obtained by clinical trial simulation [21]. This block diagonal expression, computed by default in PFIM program v4.0, is the one used to evaluate the FIM with continuous outcomes in this work (section 4). The method used in this work to evaluate the FIM for discrete outcome models is presented in Appendix A.

The D-optimality criterion $\phi_D$ is a widely used optimality criterion in the field of optimal design, which consists in maximizing the determinant of the FIM $M_{PF}(\Psi, \Xi_w)$ (see equation (3)), normalized by the number of model parameters to be estimated $P$:

$$\phi_D(\Xi_w) = |M_{PF}(\Psi, \Xi_w)|^{\frac{1}{p}} \quad (6)$$

The relative D-efficiency $E_D(\Xi_w)$ of a design $\Xi_w$, with respect to the reference design $\Xi_w^*$, is computed as:

$$E_D(\Xi_w) = \frac{\phi_D(\Xi_w)}{\phi_D(\Xi_w^*)} \quad (7)$$
### 3. Multiplicative algorithm for design optimization

The monotonic convergence of this class of multiplicative algorithm has generally been established under certain conditions by several optimality criteria, including the D-criterion [15]. In this paper, the proposed algorithm is an extension [15, 22] to NLMEMs, based on the D-optimality criterion. Moreover, this class of algorithm has been shown to converge for several applications, including design optimization with pre-specified possible design points [22], as it is the case in this manuscript.

#### 3.1. Equivalence theorem

Optimizing a population design $\Xi_w$ means optimizing the weight set $w$, i.e the weight $w_q$ of each elementary design $\xi_q$, knowing that all the possible elementary designs are defined beforehand. The equivalence theorem can be used to assess whether any given design is optimal for a given criterion [23]. The general equivalence theorem [24] states that the following conditions for a D-optimal design are equivalent, assuming the concavity of the D-criterion $\phi_D(\Xi_w)$ as formulated by Yu [15]:

1. $\Xi_w^*$ is the D-optimal design

2. $\Xi_w^*$ maximizes $\phi_D(\Xi_w) = |M_{PF}(\Psi, \Xi_w)|^{1/q}$

3. $\Xi_w^*$ minimizes $\max_{1 \le q \le Q} (d(w_q))$, $\forall q = 1, ..., Q$, where $d(w_q)$ is the maximal rate of increase of $\phi_D$ at $w_q$:

$$d(w_q) = \frac{\partial \phi_D(\Xi_w)}{\partial w_q} = tr(\frac{1}{P} |M_{PF}(\Psi, \Xi_w)|^{1/q} M_{PF}(\Psi, \Xi_w)^{-1} M_{PF}(\Psi, \xi_q)) \quad (8)$$

#### 3.2. Multiplicative algorithm and convergence criterion

The algorithm is defined as follows. The vector of weights $w = (w_1, ..., w_Q)$ is initialized as $w^{(0)}$. At each iteration $k$, each element of this vector is multiplied by a factor, depending on the results of the previous iteration:

$$w_q^{(k+1)} = m_q^{(k)} w_q^{(k)} \quad (9)$$
Iterations process until convergence to the optimal design or until the maximum number of iteration K is reached.

The multiplier factor at each iteration k is defined as:

$$m_q^{(k)} = \frac{d(w_q^{(k)})^\lambda}{\sum_{j=1}^{Q} w_j^{(k)} d(w_j^{(k)})^\lambda} \quad (10)$$

where $\lambda \in (0, 1]$, denotes a calibration parameter (in case of D-optimality, the convergence can be shown if $\lambda = 1$ [25]). This parameter is set at 0.99 in our examples.

The general equivalence theorem ensures the monotonicity of the multiplicative algorithm to a D-optimal design. Further proofs can be found in the article of Yu, especially theorem for monotonicity and for convergence [15]. Notice that at this step, for the algorithm itself the space of interest is the weight set $w$. Each element $w_q$ of this weight set corresponds to each possible elementary design $\xi_q$ given before performing the algorithm. As D-criterion $\phi_D$ is increasing and concave, Yu proposed to use the following convergence criterion: $w^{(k)}$ at the $k^{th}$ iteration is considered optimal if each weight $w_q$ satisfies the inequality:

$$\max_{1 \le q \le Q} d(w_q^{(k)}) \le (1 + \delta) \sum_{j=1}^{Q} w_j^{(k)} d(w_j^{(k)}) \quad (11)$$

where $\delta$ is a small positive constant ($\delta = 10^{-4}$ according to Yu [15]).

The D-optimal design $\Xi_w^*$ is found if at certain interation, the weight set $w$ satisfies the equation (11).

## Heuristic explanation

In order to maximize the determinant of the moment matrix, the algorithm proposes to find each optimal weight $w_q$ (corresponding to each pre-defined elementary design $\xi_q$) by multiplying it by the factor $m_q$ at each iteration, which depends on $d(w_q)$. The value of $d(w_q)$ indicates the amount of gain in information, as measured by $\phi_D$, by a slight increase in $w_q$, the weight on the $q^{th}$ elementary design. Then, the equation (9) can be seen as adjusting the
weight vector so that relatively more weight is placed on elementary designs
whose increased weight may result in a larger gain in $\phi_D$.

### 3.3. *Implementation*

Two versions of multiplicative algorithm implementation are proposed. One is implemented in an R function, the other is implemented in C (see Appendices B and C, respectively) for the loop part of the algorithm which needs to be compiled as a dynamic link library file and then to be called by R. Given that the exact CPU time varies with each time of execution and depends on the computer's hardware, values provided are meant to compare relatively with each other. Arguments of the function are presented in Figure 1 and commented codes are available in Appendices B and C. In order to account for several methods of FIM evaluation, the list of evaluated FIMs is considered as an argument of the function. Each FIM in that list quantifies the information of each elementary design generated from the design space given experimental constraints. The codes for FIM evaluation using the FO method and for the algorithm application corresponding to the illustrative example of section 4 are available in Appendices D and E, respectively. The presented function returns the weight vector at each iteration (in the form of a dynamic graph), the number of iterations required, whether the convergence of the algorithm is achieved, the optimal design and the D-criterion value corresponding to the final weight vector.

# 4. Application to design optimization for a longitudinal study with continuous outcome

## 4.1. *Context and objectives*

This example is inspired by the clinical development of a monoclonal anti-body indicated in the treatment of wet age-related macular degeneration (wet-AMD) [16]. The treatment reduces vessel leakiness and improves visual acu-ity (VA) by neutralizing the vascular endothelial growth factor (VEGF) in the
### Function in R

* multiplicative_algo(FIMs, version, designSpace, nbTot, lambda, delta, iteration)

The arguments for this function represent:

<table><tr><td>FIMs</td><td>the list of elementary FIM corresponding to all possible elementary designs according to given constraints</td></tr><tr><td>version</td><td>the version of implementation for the loop in the algorithm, either R or C</td></tr><tr><td>designSpace</td><td>the list or vector of elementary design space corresponding to the list of FIMs</td></tr><tr><td>nbTot</td><td>total number of subject considered in the study</td></tr><tr><td>lambda</td><td>lambda's value in algorithm</td></tr><tr><td>delta</td><td>the constant's value for stop condition</td></tr><tr><td>iteration</td><td>maximum number of iterations per call of function</td></tr></table>

justification=raggedright,singlelinecheck=false

<caption>Figure 1: Arguments of the multiplicative algorithm function in R.</caption>

retina. In wet-AMD, disease progression and treatment effect are assessed using a VA test based on the early treatment diabetic retinopathy study (ETDRS) chart which contains 14 lines (70 letters). The ETDRS chart is used to measure the number of letters successfully read by the patient at a given visit. The trial duration is set to 24 months, with several VA measurements per patient. The patients were randomized to different dose groups.

Both VA measurement times and doses were optimized in these patients for a hypothetical future study, based on a published model of disease progression and the anti-VEGF treatment effect, which have been used in a previous methodological study [26].

## 4.2. Model and parameters

The following model *f* describes the time course of VA (a continuous scale from 0 to 70 letters read by the patient) in wet-AMD patients, with or without anti-VEGF treatment, where *f* provides a predicted scalar value at time *t*:

$$f(\theta, d, t) = VA_0 + (1 - e^{-kt}) \times \left( \frac{E_{max} d}{(ED_{50} + d)} - \beta VA_0 \right) \quad (12)$$

Parameter vector *θ* is composed of 5 elements: $\theta = (VA_0, k, \beta, E_{max}, ED_{50})$. VA₀ is the visual acuity at baseline and in the absence of treatment. VA exponentially decreases over time *t* in untreated patients, reaching an asymptote VA₀ × (1 − *β*) at a rate *k*. Based on *k*, one can derive the average time to steady
state ($5 \times \log(2)/k$), which equals approximately 24 months and corresponds to the end of trial. The model mimics a mixture of both disease-modifying and symptomatic effects. The dose-response relationship is assumed to be an Emax function, where $E_{max}$ represents the maximum number of letters that an individual can gain and $ED_{50}$ the dose at which 50% of the maximal effect is achieved. For a very high dose and at steady state, the predicted visual acuity is $VA_0 \times (1 - \beta) + E_{max}$.

The $E_{max}$ parameter follows a normal distribution, $VA_0$, $k$, and $\beta$ follow a log-normal distribution. No inter-individual variability of $ED_{50}$ is considered. An additive residual error model is assumed (section 2.2). The parameter values are reported in Table 1.

<table><caption>Table 1: Parameter values of the model for visual acuity (with $\mu$ the fixed effects, $\omega^2$ the variance of the random effects, and $\sigma^2$ the variance of the residual error).</caption><thead><tr><th>Parameter</th><th>distribution</th><th>μ</th><th>ω<sup>2</sup></th><th>σ<sup>2</sup></th></tr></thead><tbody><tr><td>VA<sub>0</sub> (Letters)</td><td>log-normal</td><td>55</td><td>0.07</td><td>-</td></tr><tr><td>k (Day<sup>-1</sup>)</td><td>log-normal</td><td>0.005</td><td>0.5</td><td>-</td></tr><tr><td>β</td><td>log-normal</td><td>0.2</td><td>1.0</td><td>-</td></tr><tr><td>E<sub>max</sub> (Letters)</td><td>normal</td><td>30</td><td>150</td><td>-</td></tr><tr><td>ED<sub>50</sub> (μg)</td><td>log-normal</td><td>150</td><td>-</td><td>-</td></tr><tr><td>σ<sup>2</sup> (Letters<sup>2</sup>)</td><td>-</td><td>-</td><td>-</td><td>28</td></tr></tbody></table>

The predicted VA time course curves from 0 to 24 months with four dose groups (0, 50, 100 and 500 µg) are shown in Figure 2 for illustration of the model shape, using the fixed effects values of Table 1.
<figure><img src="image_15.png" /><figcaption>Figure 2: Median time course of the visual acuity (VA, in *Letters*) predicted by the model (12), for different dose groups: 0, 50, 100 and 500 µg respectively.</figcaption></figure>

### 4.3. *Design optimization*

The multiplicative algorithm was used to optimize the population design $\Xi_w$, based on the D-optimality criterion. The population FIM was evaluated using the first-order approximation of the model [9].

The population design $\Xi$ was composed of $N = 300$ patients included in total and assigned to different elementary designs $\xi_q$, chosen among a finite number $Q$ of possible elementary designs. The measurement times were optimized among 26 possible times: at baseline (day 0, fixed), day 7 and every 28 days from the 28$^{th}$ to the 672$^{nd}$ day. The optimal doses were chosen from among 7 possible doses: 0, 25, 50, 100, 150, 300 and 500 µg. Given these constraints, each elementary design $\xi_q$ was composed of a vector $t_q$ of several measurement times and one dose $d_q$: $\xi_q = (t_q, d_q)$.

First, a reference-rich design was evaluated with all 26 measurement times in 4 dose groups (0, 150, 300, 500 µg) defined in a previous study [26] and with the same number of subjects in each group ($\xi_q = (t_q = \{0, 7, \text{seq}(28, 672, 28)\}, d_q)$). In other words, this population design is composed of 4 elementary designs with the same "weight" allocated to each design (Table 2, Reference 26 times scenario).

This reference design was compared to the one with optimized dose groups
<table><caption>Table 2: Reference and optimal design for the model of visual acuity, with 26 measurement times on top and 4 measurement times on bottom with 300 subjects. With $N_{q'}$, the number of subjects in each group, Q' the number of elementary designs in the final population design, Q the number of possible elementary designs in the design space, $\phi_D$ the D-criterion, $E_D$ the D-efficiency, RSE the predicted relative standard error (for fixed effect parameters $\mu$, variance of random effect $\omega^2$ and standard error of the residual error $\sigma$) and CPU Time the computing time (FIM calculation running time + optimization running time in R and in C, respectively). Design elements in bold were optimized (for columns "Dose", "Measurement times" and "$N_{q'}$").</caption><thead><tr><th>Scenario</th><th>Dose</th><th>Measurement times</th><th>N<sub>q'</sub></th><th>Q'(Q)</th><th>φ<sub>D</sub><br>(E<sub>D</sub>%)</th><th>RSE(%) of (μ, ω<sup>2</sup>), σ</th><th>CPU Time<br>(sec)</th></tr></thead><tbody><tr><td rowspan="5">Reference<br>26 times</td><td>0</td><td>0, 7, seq(28, 672, 28)</td><td>75</td><td rowspan="5">4 (4)</td><td rowspan="5">33.7<br>(100%)</td><td>VA<sub>0</sub></td><td>1.6, 8.4</td><td rowspan="5">FIM 0.06</td></tr><tr><td>150</td><td>0, 7, seq(28, 672, 28)</td><td>75</td><td>k</td><td>7.7, 26.4</td></tr><tr><td>300</td><td>0, 7, seq(28, 672, 28)</td><td>75</td><td>β</td><td>12, 16.4</td></tr><tr><td>500</td><td>0, 7, seq(28, 672, 28)</td><td>75</td><td>E<sub>max</sub></td><td>16.2, 41.1</td></tr><tr><td></td><td></td><td></td><td></td><td>ED<sub>50</sub></td><td>47.5</td></tr><tr><td rowspan="5">Fix 26 times<br>optimize<br>doses</td><td>0</td><td>0, 7, seq(28, 672, 28)</td><td>116</td><td rowspan="5">3 (7)</td><td rowspan="5">38.1<br>(113.3%)</td><td>VA<sub>0</sub></td><td>1.6, 8.5</td><td rowspan="5">FIM 0.18<br>optim R 0.50<br>optim C 0.12</td></tr><tr><td>100</td><td>0, 7, seq(28, 672, 28)</td><td>50</td><td>k</td><td>7.2, 22.97</td></tr><tr><td>500</td><td>0, 7, seq(28, 672, 28)</td><td>134</td><td>β</td><td>9.6, 13.1</td></tr><tr><td></td><td></td><td></td><td>E<sub>max</sub></td><td>12.99, 35.8</td></tr><tr><td></td><td></td><td></td><td>ED<sub>50</sub></td><td>43.8</td></tr><tr><td rowspan="5">Reference<br>4 times</td><td>0</td><td>0, 224, 420, 672</td><td>75</td><td rowspan="5">4 (4)</td><td rowspan="5">15.5<br>(100%)</td><td>VA<sub>0</sub></td><td>1.6, 9.09</td><td rowspan="5">FIM 0.05</td></tr><tr><td>150</td><td>0, 224, 420, 672</td><td>75</td><td>k</td><td>14.8, 98.5</td></tr><tr><td>300</td><td>0, 224, 420, 672</td><td>75</td><td>β</td><td>13.6, 19.8</td></tr><tr><td>500</td><td>0, 224, 420, 672</td><td>75</td><td>E<sub>max</sub></td><td>17.5, 48.1</td></tr><tr><td></td><td></td><td></td><td>ED<sub>50</sub></td><td>50.9</td></tr><tr><td rowspan="5">Fix 4 times<br>optimize<br>doses</td><td>0</td><td>0, 224, 420, 672</td><td>113</td><td rowspan="5">3 (7)</td><td rowspan="5">17.6<br>(114.2%)</td><td>VA<sub>0</sub></td><td>1.6, 9.1</td><td rowspan="5">FIM 0.15<br>optim R 0.25<br>optim C 0.05</td></tr><tr><td>100</td><td>0, 224, 420, 672</td><td>56</td><td>k</td><td>13.6, 84.6</td></tr><tr><td>500</td><td>0, 224, 420, 672</td><td>131</td><td>β</td><td>11.2, 16.1</td></tr><tr><td></td><td></td><td></td><td>E<sub>max</sub></td><td>13.9, 41.9</td></tr><tr><td></td><td></td><td></td><td>ED<sub>50</sub></td><td>45.4</td></tr><tr><td rowspan="5">Optimize<br>both doses<br>and 4 times</td><td>0</td><td>0, 168, 196, 672</td><td>107</td><td rowspan="5">4 (1932)</td><td rowspan="5">19.4<br>(125.5%)</td><td>VA<sub>0</sub></td><td>1.6, 8.9</td><td rowspan="5">FIM 4.96<br>optim R 431.59<br>optim C 30.95</td></tr><tr><td>50</td><td>0, 7, 644, 672</td><td>13</td><td>k</td><td>11.1, 58.2</td></tr><tr><td>100</td><td>0, 7, 644, 672</td><td>48</td><td>β</td><td>11.5, 16.4</td></tr><tr><td>500</td><td>0, 168, 196, 672</td><td>132</td><td>E<sub>max</sub></td><td>13.4, 43.3</td></tr><tr><td></td><td></td><td></td><td>ED<sub>50</sub></td><td>44.4</td></tr><tr><td></td><td></td><td></td><td>σ</td><td>3.51</td><td></td><td></td><td></td></tr></tbody></table>
(among 7 proposed doses) with the 26 same measurement times.

Second, a reference sparse design was evaluated with 4 measurement times quasi-equally spaced between 0 and 672 days and with the same 4 doses groups (0, 150, 300, 500 µg) with the same number of subjects in each group ($ξ_q = (t_q = \{0, 224, 420, 672\}, d_q)$). This design was compared to two optimized designs: one optimizing only the doses among 7 possible doses: 0, 25, 50, 100, 150, 300 and 500 µg and another optimizing both times (with each possible elementary design composed of 2 fixed times: 0 and 672 days and 2 times optimized among 24 proposed times between 7 and 644 days) and doses levels among 7 possible doses: 0, 25, 50, 100, 150, 300 and 500 µg ($ξ_q = (t_q = \{t_1, t_2, t_3, t_4\}, d_q)$).

The multiplicative algorithm parameters were set as follows: $λ = 0.99$, $δ = 0.00001$ and maximum iteration number 10000.

All results were compared in terms of doses, measurement times, number of subjects per group ($N_q$), D-criterion, D-efficiency with respect to the reference design, relative standard errors (%) of all population parameters, number of possible elementary designs in the design space (Q), number of elementary designs in the final population design (Q') and CPU time for running the optimization program. The multiplicative algorithm program can be separated into two parts: first it evaluates FIMs of each elementary design and saves them, then it estimates weights of these elementary designs in regard to their FIMs. CPU time was therefore reported separately for evaluating FIMs and for estimating empirical weights.

## 4.4. Results

Table 2 reports results of all the designs evaluated and optimized. With a rich design (26 measuring times), optimization of the doses led to a population design composed of 3 groups, which was, as expected, better than the reference design (D-criterion of 38.12 vs 33.66 with the same total number of subjects) and therefore led to a slightly better precision of estimation (for example, the relative standard error (RSE) of $ED_{50}$ decreased from 48% with the non-optimized design to 44% with the optimal design).
Reducing the number of measurements from 26 to 4 in the two reference designs led to a loss of half of the information measured by the D-optimality criterion (from 33.66 to 15.45) whereas the number of measurements was divided by 6.5. With fixed measurement times, the information gain was approximately 14% after optimization of doses. When both doses and measurement times were optimized (last scenario), 1932 elementary designs (or groups) were pre-defined. After optimization of the weights corresponding to these groups, four elementary designs with non-null number of subjects were finally included in the population design. With this design the gain in information was up to 25% as compared to the reference 4 times design (with non-optimized weights). With the optimal design in both dose and times, the estimation precision was also improved for all parameters, and RSE values were below 20% except for ED₅₀ (44%). However, the predicted RSE value of ED₅₀ was lower than in the non-optimized rich design with 26 samples, whereas the number of visits was 6.5 times lower (with the same total number of subjects).

With respect to CPU time running for the optimization, the more elementary designs were involved (i.e. the larger *Q*), the longer the function ran. In addition, running time for the evaluation of FIMs was lower than the optimization running time, especially when implemented in R and when *Q* was large. As seen in Table 2, implementation in C ran much faster than R. The full code of this application can be found in Appendix E.

## 5. Discussion

In this work, a multiplicative algorithm based on the FIM is proposed and implemented for design optimization for studies with longitudinal data analyzed with nonlinear mixed effect models. The advantages of this algorithm are its simplicity and that it has proof of convergence for several optimality criteria [15]. It increases information after each iteration, which guarantees the quality of its optimal design result, given clinical constraints.

This algorithm was applied to two examples. The first application aims
to design a new study by optimizing both measurement times and doses in a study of a monoclonal antibody indicated in the treatment of wet-AMD. To our knowledge, this is the first time that the multiplicative algorithm is used into such a model with random effects and on top of that, both time measurements and dose administration are optimized. The second application, presented in Appendix A, is an additional test of the algorithm and makes sure that this extension of the algorithm can be applied to any model. This is also the first application of this type of algorithm to a logistic model with repeated discrete data. Here, two specific models and multivariate normal distribution of random effects were considered. However, the multiplicative algorithm could be applied with any structural model and any distribution of random effects as long as the FIM can be computed [27].

The implemented algorithm uses the D-optimality criterion. However, other criteria can be considered. In the A-criterion, the total variance of the parameter estimates is minimized, and this is equivalent to minimizing the average variance [6]. If specific functionals of the parameter, such as the area under the curve, the minimum effective dose or the maximum of the curve, are of interest, the C-criterion is used ($\phi_C = -c^T M c$, where c is a non-zero constant vector) [15, 28]. The Ds-criterion can also be of interest, in order to estimate a subset of the parameters as precisely as possible [6]. They all have proved ability to enable the algorithm to converge to the optimal result monotonically. A straightforward extension of the present work would be to use the algorithm with these different criteria.

Here, we only considered locally optimal designs, which depend on knowledge of the model and parameters. Another possible extension of the multiplicative algorithm is to combine this algorithm with robust criteria, with respect to the model and/or parameters [29, 18].

In conclusion, the extension of the multiplicative algorithm to mixed effect models proposed here allows efficient optimization of population designs, with both continuous and discrete outcomes. The presented version of the algorithm is currently being included in the next release of PFIM.
## 6. Acknowledgement and declarations

The authors thank Dr. Sebastian Ueckert (Department of Pharmaceutical Biosciences, Uppsala University, Uppsala, Sweden) for providing the R code implementing the QRMC-AGQ (quasi-random Monte Carlo - adaptive Gaussian quadrature) method [30].

The authors have no conflict of interest to declare.

## References

[1] P. Van Der Graaf, CPT: Pharmacometrics and Systems Pharmacology, CPT: Pharmacometrics & Systems Pharmacology 1 (9) (2012) e8.

[2] D. Mould, R. Upton, Basic concepts in population modeling, simulation, and model-based drug development, CPT: Pharmacometrics & Systems Pharmacology 1 (9) (2013) e6.

[3] D. Mould, R. Upton, Basic concepts in population modeling, simulation, and model-based drug development - part 2: introduction to pharmacodynamic modeling methods, CPT: Pharmacometrics & Systems Pharmacology 2 (4) (2014) e38.

[4] R. Upton, D. Mould, Basic concepts in population modeling, simulation, and model-based drug development - part 3: introduction to pharmacodynamic modeling methods, CPT: Pharmacometrics & Systems Pharmacology 3 (1) (2014) e88.

[5] F. Mentré, M. Chenel, E. Comets, J. Grevel, A. Hooker, M. O. Karlsson, M. Lavielle, I. Gueorguieva, Current use and developments needed for optimal design in pharmacometrics: a study performed among DDMoRe's european federation of pharmaceutical industries and associations members, CPT: Pharmacometrics & Systems Pharmacology 2 (2013) e46.
[6] A. C. Atkinson, A. N. Donev, R. D. Tobias, Optimum Experimental Designs, with SAS, Oxford Statistical Science, Oxford University Press 2009, Oxford.

[7] V. V. Fedorov, S. L. Leonov, Optimal Design for Nonlinear Response Models, Chapman and Hall/CRC Press 2014, Boca Raton.

[8] L. Pronzato, A. Pázman, Design of Experiments in Nonlinear Models, Springer 2013, New York.

[9] C. Dumont, G. Lestini, H. Le Nagard, F. Mentré, E. Comets, T. T. Nguyen, PFIM Group, PFIM 4.0, an extended R program for design evaluation and optimization in nonlinear mixed-effect models, Computer Methods and Programs in Biomedicine 156 (2018) 217-229.

[10] J. Nelder, R. Mead, A Simplex method for function minimization, Computer Journal 1 (1965) 308-313.

[11] V. Fedorov, Theory of Optimal Experiments, Academic Press 1972, New York.

[12] H. P. Wynn, Results in the theory and construction of D-optimum designs, Journal of the Royal Statistical Society B 34 (1972) 133-147.

[13] S. Retout, E. Comets, A. Samson, F. Mentré, Design in nonlinear mixed effects models: optimization using the Federov-Wynn algorithm and power of the Wald test for binary covariates, Statistics in Medicine 26 (2007) 5162-5179.

[14] R. Harman, Computing D-optimal approximate designs of experiments on finite spaces: a survey and comparison of algorithms, Presented at Workshop on Design of Experiments: New Challenges, Marseilles, France (2018).

[15] Y. Yu, Monotonic convergence of a general algorithm for computing optimal designs, Annals of Statistics 38 (3) (2010) 1593-1606.
[16] P. J. Rosenfeld, D. M. Brown, J. S. Heier, D. S. Boyer, P. K. Kaiser, C. Y. Chung, R. Y. Kim, Ranibizumab for neovascular age-related macular degeneration, New England Journal of Medicine 355 (14) (2006) 1419-1431.

[17] G. Lestini, S. Ueckert, F. Mentré, Robust design in model based analysis of longitudinal clinical data, Presented at PODE meeting, Uppsala, Sweden (2016).

[18] J. Seurat, T. T. Nguyen, F. Mentré, Robust designs accounting for model uncertainty in longitudinal studies with binary outcomes, Statistical Methods in Medical Research 29 (3) (2020) 934-952.

[19] F. Mentré, A. Mallet, D. Baccar, Optimal design in random effect regression models, Biometrika 84 (1997) 429-442.

[20] T. Mielke, R. Schwabe, Some considerations on the fisher information in nonlinear mixed effects models, mODa 9-Advances in Model-Oriented Design and Analysis (2010) 129-136.

[21] J. Nyberg, C. Bazzoli, K. Ogungbenro, A. Aliev, S. Leonov, S. Duffull, A. Hooker, F. Mentré, Methods and software tools for design evaluation for population pharmacokinetics-pharmacodynamics studies, British Journal of Clinical Pharmacology 79 (2014) 6-17.

[22] T. Holland-Letz, On the combination of c- and D-optimal designs: general approaches and applications in dose-response studies, Biometrics 73 (1) (2017) 206-213.

[23] F. Pukelsheim, Optimal design of experiments, SIAM 2006, Philadelphia.

[24] P. Whittle, Some general points in the theory of optimal experimental design, Journal of the Royal Statistical Society (1973).

[25] A. Pázman, Foundations of optimum experimental design, Springer 1986, Dordrecht.
[26] S. Buatois, S. Ueckert, N. Frey, S. Retout, F. Mentré, Comparison of model averaging and model selection in dose finding trials analyzed by nonlinear mixed effect models, American Association of Pharmaceutical Scientists Journal 20 (3) (2018) 56.

[27] C. Proust, H. Jacqmin-Gadda, Estimation of linear mixed models with a mixture of distribution for the random effects, Computer Methods and Programs in Biomedicine 78 (2) (2005) 165-173.

[28] T. Holland-Letz, H. Dette, D. Renard, Efficient algorithms for optimal designs with correlated observations in pharmacokinetics and dose-finding studies, Biometrics 68 (1) (2012) 138-145.

[29] F. Loingeville, T. T. Nguyen, M.-K. Riviere, F. Mentré, Robust designs in longitudinal studies accounting for parameter and model uncertainties—application to count data, Journal of Biopharmaceutical Statistics 30 (1) (2020) 31-45.

[30] S. Ueckert, F. Mentré, Computation of the Fisher information matrix for discrete nonlinear mixed effect models, Presented at 8th International Conference on Computational and Methodological Statistics, London, UK (2015).

[31] T. T. Nguyen, F. Mentré, Evaluation of the Fisher information matrix in nonlinear mixed effect models using adaptive Gaussian quadrature, Computational Statistics & Data Analysis 80 (2014) 57-69.

[32] S. Ueckert, F. Mentré, A new method for evaluation of the Fisher information matrix for discrete mixed effect models using Monte Carlo sampling and adaptive Gaussian quadrature, Computational Statistics & Data Analysis 111 (2016) 203-219.
# Appendices

## Appendix A: Application to design optimization for a longitudinal study with binary outcome

### *Nonlinear mixed effect models for discrete data*

For discrete data, the conditional probability for the vector of observations $y_i$ of individual $i$ is modeled as follows:

$$p(y_i|b_i, z_i, \Psi) = h(y_i, \xi_i, \theta_i) \quad (13)$$

In equation (13), $h$ is the function describing the probability of $y_i$ given individual random effects $b_i$ for a design $\xi_i$ and a vector of subject-specific parameters $\theta_i$. An example of $h$ is given below in this Appendix (with $h$ a logistic model). $\theta_i$ is expressed as a function $g(\mu, b_i, z_i, \beta)$ of: the vector of fixed effects $\mu$, the vector of random effects $b_i$ and the vector of covariates $z_i$ for individual $i$, as well as the vector of covariate effects $\beta$. $b_i$ follows a normal distribution with mean zero and variance-covariance matrix $\Omega$ which accounts for the inter-individual variability, i.e. $b_i \sim N(0, \Omega)$. $\Omega_u$ is the vector containing all unique elements of $\Omega$. The vector $\Psi^T$ contains all population parameters to be estimated ($\Psi^T = (\mu^T, \beta^T, \Omega_u^T)$).

### *Fisher information matrix by quasi random Monte Carlo-adaptive Gaussian quadrature*

When linearization is to be avoided, for example in models with discrete outcomes, alternatives such as adaptive Gaussian quadratures can be used [31, 32]. To evaluate the FIM (equation (4)), one needs to compute two types of integrations: one with respect to $b_i$ and one with respect to $y_i$. Ueckert and Mentré [32] evaluated the former by AGQ and the latter by MC and showed by CTS that this is a suitable approach with a moderate number of random effects. They then extended this approach by using quasi-random Monte Carlo (QRMC) instead of MC, which leads to huge gains in computing time [30]. The principle of quasi-random MC integration is the following: instead of naively
sampling replicates of $b_i$ from $p(b_i|\Psi)$ then sample replicates of $y_i$ from $p(y_i|b_i)$, samples of $b_i$ in $[0, 1]^p$ were generated from a Sobol sequence and samples were transformed using $p^{-1}(b_i|\Psi)$ and $p^{-1}(y_i, b_i)$.

This approach implemented in an R code is the one used to evaluate the FIM with binary outcomes in this work.

### Context and objectives

This example, inspired by a previous work [17, 18], was a one-year clinical trial with repeated binary outcomes (responder vs non-responder) in two balanced treatment groups. The response probability over time was described by a mixed logistic model, including a covariate for the treatment effect. The objective in this application was to determine the optimal location of a limited number of measurement times chosen from 0 to 12 months.

### Model and parameters

The model studied describes binary responses ($y = 0$ or $1$) recorded over time ($t = 0$ to 12 months) in two different treatment groups. The response probability at time $t$ is given as

$$p(y = 1|b, 1_T) = \frac{\exp(\theta_1 + \theta_2(1 + \beta \times 1_T)t)}{1 + \exp(\theta_1 + \theta_2(1 + \beta \times 1_T)t)}, \quad (14)$$

where $\theta_1$ is the baseline and $\theta_2$ the slope parameter. $1_T$ is a treatment group indicator variable (with $1_T = 0$ if control group, $1_T = 1$ if experimental group) and $\beta$ is the treatment effect on $\theta_2$. $\theta_1$ and $\theta_2$ follow a normal distribution with a diagonal variance matrix: $\theta_l = \mu_l + b_l$ where $b_l \sim N(0, \omega_l^2)$ for $l = \{1, 2\}$. Figure 3 shows the predicted response curve from 0 to 12 months for each treatment group using nominal values in Table 3.

### Design optimization

The multiplicative algorithm was used to optimize the population design $\Xi_w$, based on the D-optimality criterion. The population FIM was evaluated
<table><caption>Table 3: Parameter values of the model for repeated binary outcome (with μ the fixed effects, ω² the variance of the random effects, and β the treatment effect).</caption><thead><tr><td>Parameter</td><td>distribution</td><td>μ</td><td>ω<sup>2</sup></td><td>β</td></tr></thead><tbody><tr><td>θ<sub>1</sub></td><td>normal</td><td>-2</td><td>0.49</td><td>-</td></tr><tr><td>θ<sub>2</sub></td><td>normal</td><td>0.09</td><td>0.03</td><td>-</td></tr><tr><td>β</td><td>-</td><td>-</td><td>-</td><td>5</td></tr></tbody></table>

<figure><img src="image_16.png" /><figcaption>Figure 3: Mean of the logit of the response probability over 12 months predicted by the model (14). The corresponding parameter values are indicated in Table 3. Solid lines represent the control group and dashed lines the experimental group.</figcaption></figure>

using the QRMC/AGQ method. The number of QRMC samples was fixed to 1000 and the number of AGQ nodes to 7.

The population design $\Xi$ was composed of $N = 100$ patients (50 in each treatment group) included in total and assigned to different elementary designs $\xi_q$, chosen among a finite number $Q$ of possible elementary designs. Measurement times were optimized among several possible times between the beginning (at month 0, fixed) and the end of the study (month 12, fixed). Given these constraints, each elementary design was composed of a vector $t_q$ of four measurement times ($\xi_q = t_q = \{0, t_2, t_3, 12\}$)).

First, a reference equispaced non-optimized design was evaluated with 4 measurement times ($\xi_q = (t_q = \{0, 4, 8, 12\})$) and $N = 100$ subjects: 50 in the
two treatment groups.

Second, two out of the four measuring times of the design ($\xi_q = (t_q = \{t_1, t_2, t_3, t_4\})$) were optimized, with the same number of subjects (and same measurement times) in the two treatment groups. Assuming that $t_1$ and $t_4$ were fixed to 0 and 12 months (i.e. the beginning and the end of the study respectively), only two times $t_2$ and $t_3$ were optimized among the set of 11 possible times (1,2,3,4,5,6,7,8,9,10,11). No repetition was considered, i.e. four sampling times must be different.

The multiplicative algorithm parameters were set as follows: $\lambda = 0.99$, $\delta = 0.00001$ and maximum iteration number 10000.

The optimal design was thus compared to the reference equispaced design in terms of measurement times, D-efficiency, RSE (notably the treatment effect $\beta$, as the expected Wald test power to detect this covariate effect depends on the predicted SE of $\beta$) and CPU time. The D-efficiency was calculated in regard to the equispaced design.

## Results

Table 4 reports results of the evaluated and optimized designs.

According to the pre-defined clinical constraints, optimal weights of 55 elementary designs were searched by the algorithm. The optimal design was composed of three groups which differed in terms of measurement times. These three groups were respectively composed of 52, 6 and 42 individuals (26, 3 and 21 individuals in each treatment arm, placebo or treated). Compared to the reference design (non-optimized), the information of the optimized design was increased by 18% (as the D-criterion was 34.3 for the equispaced design ($t_2 = 4$ months and $t_3 = 8$ months) and 40.6 for the optimized design), with the same total number of individuals and the same number of measurement times per individual. With the optimal design, predicted RSE was notably improved for the baseline $\theta_1$ and for the treatment effect $\beta$ parameters.

Compared with the other application (section 4), the CPU time increased, mainly because of the evaluation FIM method (QRMC/AGQ, detailed above.
<table><caption>Table 4: Reference and optimal design for the model for repeated binary outcome, with 4 measurement times and 100 subjects (50 in each treatment group). With $N_{q'}$ the number of subjects, Q' the number of elementary designs in the final population design, Q the number of possible elementary designs in the design space, $\phi_D$ the D-criterion, $E_D$ the D-efficiency, RSE the predicted relative standard error (for fixed effect parameters $\mu$, variance of random effect $\omega^2$ and treatment effect $\beta$) and CPU Time the computing time (FIM calculation running time + optimization running time in R and in C respectively). Design elements in bold were optimized (for columns "Measurement times" and "$N_{q'}$").</caption><thead><tr><th>Scenario</th><th>Measurement times</th><th>N<sub>q'</sub></th><th>Q'(Q)</th><th>φ<sub>D</sub><br>(E<sub>D</sub>%)</th><th colspan="2">RSE(%) of (μ, ω<sup>2</sup>), β</th><th>CPU Time<br>(sec)</th></tr></thead><tbody><tr><td>Reference<br>4 times</td><td>0, 4, 8, 12</td><td>100</td><td>1 (1)</td><td>34.3<br>(100%)</td><td>θ<sub>1</sub><br>θ<sub>2</sub><br>β</td><td>16.1, 156.5<br>35.8, 62.6<br>46.9</td><td>FIM 43.4</td></tr><tr><td>Optimize<br>4 times</td><td><b>0, 2, 3, 12</b><br><b>0, 2, 9, 12</b><br><b>0, 3, 7, 12</b></td><td><b>52</b><br><b>6</b><br><b>42</b></td><td>3 (55)</td><td>40.6<br>(118.4%)</td><td>θ<sub>1</sub><br>θ<sub>2</sub><br>β</td><td>14.6, 126<br>36.8, 67.6<br>42</td><td>FIM 4588.85<br>optim R 4.68<br>optim C 0.16</td></tr></tbody></table>

Once the FIM was evaluated for the 55 designs, the multiplicative algorithm performed in 4.68s by program in R, which reduced to 0.16s in C.
Appendix B: Multiplicative function in R

```R
1 multiplicative_algo<-function (FIMs, version, nbParameters,
                                   designSpace, nbSubjects, lambda, delta, iteration) {

  2 numberFIMs <- length(FIMs)
  3 w<- rep (1/numberFIMs, numberFIMs)
  4
  5 if( version == "R")
  6 {
  7   for (k in 1:iteration) {
  8     #Calculate the sum of matrix with weight associated
  9     #idea: Mw<-Mw+w[i]*FIMs[[i]]
60    MFw<-Map("*",FIMs,w) #associate weights
61    Mw<-Reduce("+",MFw) #sum of all matrices
62    dm<-dim(Mw) [1]
63    Dphi <- det(Mw)^(1/dm) * solve(Mw)/dm # calculate
64    derivatives of function phi_D
65    d<-sapply(FIMs, FUN = function (x) sum(diag(Dphi %*% x)))
66    #the vector of multiplier
67    w<- w*d^lambda/sum(w*d^lambda) #develop weight
68    if(max(d)<(1+delta)*sum(w*d)){break} #stop criterion
69   }
70   Dcriterion <- det(Mw*nbSubjects)^(1/dm)
71 }
72 if( version == "C")
73 {
74   dyn.load("algomult.dll")
75   q<- .C("multiplicativeAlgorithm", originMatrices=as double(
76   unlist(FIMs)), weights = as double(w), dim=as integer(dim(FIMs
77   [[1]]) [1]),
78   n=as integer(numberFIMs), dimA = as integer nbParameters
79   ), lambda = as double(lambda), delta = as double delta),
80   it = as integer(iteration))
81   dyn.unload("algomult.dll")
82   w <- q$weights
83   k <- q$it
84   MFw<-Map("*",FIMs,w)
```
Mw<-Reduce("+",MFw)
Dcriterion <- det(Mw*nbSubjects)^(1/q$dim)
}

#by default, we take into count weights which give meaningful
subject numlber

v<-which(w>1/nbSubjects) # v<-which(w>mean(w)) ##if we want to
get weights bigger than the mean

38 cat("\n\n*****************************\n")
39 cat("*** Result of Algorithm **************\n")
40 cat("\n\n IMPORTANT SUBJECT NUMBER ****************\n", w[v]*
nbSubjects)
41 cat("\n\n IMPORTANT DESIGN VALUES ****************\n")
42 print(designSpace[,v])
43
44 cat("\n\n NUMBER OF INTERATIONS ****************\n",k,
"\n\n D-CRITERION ************************\n",Dcriterion)
45
46
47 plot(w,ylim=c(0,1))
48 return(list(w,v))
49
50 }

51
Appendix C: Multiplicative function in C

/* FILE NAME: algomult.cpp

*

* PURPOSE: This file contains function which implement the multiplicative algorithm. (Yaming Yu 2010)

*

* DEVELOPERS: PFIM group, pfim@inserm.fr

* Yuxin Tang, INSERM UMR 1137 BIPID team

*

* INPUT VARIABLES FOR THE ENTRY POINT FUNCTION (multiplicativeAlgorithm) FROM R:

* Name Type Description

* __________ __________ __________

* originMatrices double* a list of fisher information matrices (FIM) computed by R functions

* and read as a vector of values column by column in the following functions

* weights double* a vector of weight corresponding to input FIMs and will evolve in the end

* of each iteration of algorithm

* dim int* the dimension of matrices (all matrices have same dimension)

* n int* total number of input matrices

* dimA int* the dimension of block A of matrix ( all matrices have same dimension of

* block A)

* lambda double* a parameter of algorithm between 0 and 1

* delta double* a precision of stop criterion of algorithm

* iteration int* the algorithm maximum number of iteration

ALGORITHM

* associate each origin FIM input from R with a weight value
```c
27 * sum all weighted matrices as one matrix
28
29 * calculate determinant of the summed matrix, and also separately
30 the determinant of block A and block B
31 of the summed matrix
32
33 * calculate a coefficient as (D-criterion of the summed matrix)
34 divided by (the total dimension of a matrix)
35
36 * inverse the summed matrix
37
38 * multiply each element of the inverse matrix with the
39 coefficient, the matrix obtained is the derivative
40 of phi
41
42 * multiply every origin FIM input from R with the derivative of
43 phi, then compute the trace of each matrix
44 product
45
46 * evolve weights by multiply corresponding (trace^lambda), and
47 divided by its sum
48
49 * algorithm stops when criterion about delta is satisfied
50
51 #include "algomult.h"
52
53 #ifdef _WINDOWS
54
55 #include <Windows.h>
56
57 // windows dll entry point
58
59 BOOL APIENTRY DllMain( HMODULE hModule, DWORD ul_reason_for_call,
```
```c
LPVOID lpReserved)
{
    switch (ul_reason_for_call)
    {
        case DLL_PROCESS_ATTACH:
        case DLL_THREAD_ATTACH:
        case DLL_THREAD DetACH:
        case DLL_PROCESS DetACH:
            break;
    }
    return TRUE;
}

```
*/
void sum(double *matrices, int size, int n, double *result)
{
    int i, j;

    // initialise the result matrix
    // copy the first matrix in the list to the result matrix
    for (j = 0; j < size; j++) {
        *result = *matrices;
        result++;
        matrices++;
    }
    result -= size;

    // add rest matrices in the list to the result matrix
    for (i = 1; i < n; i++) {
        for (j = 0; j < size; j++) {
            *result += *matrices;
            result++;
            matrices++;
        }
        result -= size;
    }
}

/* FUNCTION NAME: MatrixTimesWeight
*
* DESCRIPTION: This function multiplies each element in the matrix
by a weight.
*
* INPUT:
* Name		Type		Description
* _________
* matrices	double*		a vector of n matrices column
by colunm of dimension dim
* weights	double*		n values, one for each matrix
* size		int	
total number of element of a
matrix (size = dim * dim)
* n		int	
total number of matrices
```c
128 * 
129 * RETURN VALUE:
130 * Name		Type		Description
131 * ________
132 * result		double* 	weighted matrices, same number
133 			of matrices as input
134 */
135 void MatrixTimesWeight(double *matrices, double *weights, int size,
136 int n, double *result)
137 {
138 int size_n = size * n;
139
140 for (int i = 0; i < size_n; i++)
141 {
142   *result = (*matrices) * (*weights);
143   matrices++;
144   result++;
145   if ((i + 1) % size == 0)
146   weights++;
147 }
148 }
149
150 /* FUNCTION NAME: decompositionLU
151 *
152 * DESCRIPTION: This function decompose a matrix as a lower
153 triangular matrix and an upper triangular matrix.
154 *
155 * INPUT:
156 * Name		Type		Description
157 * ________
158 * matrix		double* 	 symmetric matrix, can be
159 		divided into two blocks independent non-zero
160 * 						 (block A,B)
161 * dim						 dimension of an entire matrix
162 * dimA					 dimension of block A of the
163 	matrix (dimension of block B = dim - dimA)
```
```c
        } while (i < dimA);

        // decomposition LU for block B, if dimA < dim
        while (i < dimA) {
            i_dim = i * dim;
            for (j = dimA; j < i; j++) {
                jj_dim = j * dim + j;
                if (U[jj_dim] == 0) {
                    i = dim;
                    break;
                }
                i_j_dim = i_dim + j;
                L[i_j_dim] = matrix[i_j_dim];
                for (k = dimA; k < j; k++) {
                    L[i_j_dim] -= L[i_dim + k] * U[k * dim + j];
                }
                L[i_j_dim] /= U[jj_dim];
            }

            for (j = i; j < dimA; j++) {
                i_j_dim = i_dim + j;
                U[i_j_dim] = matrix[i_j_dim];
                for (k = dimA; k < i; k++) {
                    U[i_j_dim] -= L[i_dim + k] * U[k * dim + j];
                }
            }
            i++;
        };

        }

        /* FUNCTION NAME: inverse
         *
         * DESCRIPTION: This function inverses a symmetric matrix
         *
         * INPUT:
```
234 * Name Type Description
235 * ________ ________ ________
236 * m double* a symmetric matrix to be
237 * ________ ________ ________
238 * dim int dimension of an entire matrix
239 * dimA int dimension of block A of the
240 * ________ ________ ________
241 * detA double determinant of block A
242 * detB double determinant of block B
243 *
243 * RETURN VALUE:
244 * Name Type Description
245 * ________ ________ ________
246 * invM double* inversed matrix
247 */
248 void inverse(double *m, int dim, double *invM, int dimA, double
249 detA, double detB)
250 {
251 int i, j, ii, jj, il, j1, dim1, il_dim, ii_dim, k;
252 double det, *c, *d, *l, *u;
253 // prepare matrix L and U
254 l = (double *)calloc(dim * dim, sizeof(double));
255 u = (double *)calloc(dim * dim, sizeof(double));
256 // initialise matrix L with all diagonal
257 for (i = 0; i < dim; i++)
258 l[i*dim + i] = 1;
259 //calculate inverse for block A
260 dim1 = dimA - 1;
261 c = (double *)calloc((dim1) * (dim1), sizeof(double));
262
263 for (j = 0; j < dimA; j++) {
264 for (i = 0; i < dimA; i++) {
265 il = 0;
266 for (ii = 0; ii < dimA; ii++) {
267 if (ii == i)
268
```c
269         continue;
270     i1_dim = i1 * dim1;
271     ii_dim = ii * dim;
272     j1 = 0;
273     for (jj = 0; jj < dimA; jj++) {
274         if (jj == j)
275             continue;
276         *(c + i1_dim + j1) = *(m + ii_dim + jj);
277         j1++;
278     }
279     i1++;
280 }
281
282     decompositionLU(c, l, u, dim1, dim1);
283     det = 1;
284     for (k = 0; k < dim1; k++)
285         det *= u[k * dim1 + k];
286     invM[i + j * dim] = pow(-1.0, i + j + 2.0) * (det) / (detA);
287 }
288 }
289
290 //calculate inverse for block B
291     dim1 = dim - dimA - 1;
292     d = (double *)calcol((dim1) * (dim1), sizeof(double));
293
294     for (j = dimA; j < dim; j++) {
295         for (i = dimA; i < dim; i++) {
296             i1 = 0;
297             for (ii = dimA; ii < dim; ii++) {
298                 if (ii == i)
299                     continue;
300             i1_dim = i1 * dim1;
301             ii_dim = ii * dim;
302             j1 = 0;
303             for (jj = dimA; jj < dim; jj++) {
304                 if (jj == j)
305                     continue;
306             }
```
```c
307         *(d + i1_dim + j1) = *(m + ii_dim + jj);
308         j1++;
309         }
310         i1++;
311     }
312
313     decompositionLU(d, l, u, dim1, dim1);
314     det = 1;
315     for (k = 0; k < dim1; k++)
316         det *= u[k * dim1 + k];
317     invM[i + j * dim] = pow(-1.0, i + j + 2.0) * (det) / (detB);
318
319     }
320 }
321     free(c);
322     free(d);
323     free(l);
324     free(u);
325 }
326
327
328     /* FUNCTION NAME: computeTrace
329     *
330     * DESCRIPTION: This function multiplies each matrix in the list by
331                  the same matrix, then compute each
332                  product's trace
333     *
334     * INPUT:
335     * Name                Type                  Description
336     * ____________________  ____________________  ____________________
337     * matrices           double*              a list of matrices to produce
338     * with theMatrix
339     * theMatrix           double*              the matrix to be multiplied
340     * with each matrices in the list
341     * dim                 int                   dimension of an entire matrix
342     * n                   int                   number of matrices in the list
343     *
344     * RETURN VALUE:
```
<table><thead><tr><th></th><th>Name</th><th>Type</th><th>Description</th></tr></thead><tbody><tr><td>342</td><td>*</td><td></td><td></td></tr><tr><td>343</td><td>*</td><td></td><td></td></tr><tr><td>344</td><td>trace<br>product</td><td>double*</td><td>the trace vector of each matrix</td></tr><tr><td>345</td><td>*/</td><td></td><td></td></tr><tr><td>346</td><td>void computeTrace(double *matrices, double *theMatrix, int dim, int n, double *trace)</td><td></td><td></td></tr><tr><td>347</td><td>{</td><td></td><td></td></tr><tr><td>348</td><td>int im, il, in, i;</td><td></td><td></td></tr><tr><td>349</td><td>double *productMatrices = NULL;</td><td></td><td></td></tr><tr><td>350</td><td>int size = dim * dim, i_size;</td><td></td><td></td></tr><tr><td>351</td><td>// intermediate matrices</td><td></td><td></td></tr><tr><td>352</td><td>productMatrices = (double *)calloc(n * size, sizeof(double));</td><td></td><td></td></tr><tr><td>353</td><td></td><td></td><td></td></tr><tr><td>354</td><td>for (i = 0; i &lt; n; i++) {</td><td></td><td></td></tr><tr><td>355</td><td>*trace = 0;</td><td></td><td></td></tr><tr><td>356</td><td>for (in = 0; in &lt; dim; in++)</td><td></td><td></td></tr><tr><td>357</td><td>{</td><td></td><td></td></tr><tr><td>358</td><td>for (im = 0; im &lt; dim; im++)</td><td></td><td></td></tr><tr><td>359</td><td>{</td><td></td><td></td></tr><tr><td>360</td><td>*productMatrices = 0;</td><td></td><td></td></tr><tr><td>361</td><td>for (il = 0; il &lt; dim; il++)</td><td></td><td></td></tr><tr><td>362</td><td>{</td><td></td><td></td></tr><tr><td>363</td><td>*productMatrices += (*matrices) * (*theMatrix);</td><td></td><td></td></tr><tr><td>364</td><td>matrices++;</td><td></td><td></td></tr><tr><td>365</td><td>theMatrix += dim;</td><td></td><td></td></tr><tr><td>366</td><td>}</td><td></td><td></td></tr><tr><td>367</td><td>matrices -= dim;</td><td></td><td></td></tr><tr><td>368</td><td>theMatrix -= size;</td><td></td><td></td></tr><tr><td>369</td><td>if (im == in)</td><td></td><td></td></tr><tr><td>370</td><td>*trace += *productMatrices;</td><td></td><td></td></tr><tr><td>371</td><td>productMatrices++;</td><td></td><td></td></tr><tr><td>372</td><td>theMatrix++;</td><td></td><td></td></tr><tr><td>373</td><td>}</td><td></td><td></td></tr><tr><td>374</td><td>matrices += dim;</td><td></td><td></td></tr><tr><td>375</td><td>theMatrix -= dim;</td><td></td><td></td></tr><tr><td>376</td><td>}</td><td></td><td></td></tr><tr><td>377</td><td>trace++;</td><td></td><td></td></tr></tbody></table>
```c
}
productMatrices -= n * size;
free(productMatrices);
}

/* FUNCTION NAME: multiplicativeAlgorithm

 * DESCRIPTION: This is the main function of multiplicative
 algorithm, which is exported and used in PFIM

 * implemented in R-S4

 * INPUT:
 * Name		Type		Description
 * ____________

 * originMatrices	double* 		a list of matrices computed by
 R functions

 * weights	double* 		a vector of weight
 corresponding to each matrix in
 the originMatrices list

 * dim		int		dimension of an entire matrix (
 all same for matrices in the list)

 * n		int* 		number of matrices in the list
 * dimA		int* 		dimension of block A
 * lambda	double* 	convergence power of the
 algorithm (0, 1]

 * delta	double* 	algorithm stop criterion
 parameter

 * iteration	int* 	algorithm maximum iteration

 *

 * RETURN VALUE:
 * Name		Type		Description
 * ____________

 * weights	double* 	the vector of weights evolved
 by algorithm

 * iteration	int* 	real iteration run number

 */

void multiplicativeAlgorithm(double *originMatrices, double *
weights, int *dim,
```
```c
408 int *n, int *dimA, double *lambda, double *delta, int *iteration
409 }
410 double det, coef, sumw, sumwd, maxTrace, *weightedMatrices, *
411 sumWeightedMatrices,
412 *inverseSumWeightMat, *Dphi, *tr, *l, *u, detA, detB;
413 int idim = *dim, in = *n, idimA = *dimA, isize = idim * idim, i,
414 it;
415
416 weightedMatrices = (double *)calloc(in * isize, sizeof(double));
417 if (!weightedMatrices) {
418   it = -2;
419   return;
420 }
421 sumWeightedMatrices = (double *)calloc(isize, sizeof(double));
422 inverseSumWeightMat = (double *)calloc(isize, sizeof(double));
423 Dphi = (double *)calloc(isize, sizeof(double));
424 tr = (double *)calloc(in, sizeof(double));
425 l = (double *)calloc(isize, sizeof(double));
426 u = (double *)calloc(isize, sizeof(double));
427
428 // initialise diagonal elements of the low matrix for
429 decomposition LU as 1
430 for (i = 0; i < idim; i++)
431   l[i*idim + i] = 1;
432
433 for (it = 0; it < (*iteration); it++)
434 {
435   // associate each FIM calculated in R with a weight value
436   MatrixTimesWeight幸originMatrices, weights, isize, in,
437   weightedMatrices);
438   // sum all weighted matrices
439   sum(weightedMatrices, isize, in, sumWeightedMatrices);
440
441 //calculate determinant of Mw, and separately of block A and B
442 of Mw
443 detA = 1;
444 detB = 1;
445
```
decompositionLU (sumWeightedMatrices, l, u, idim, idimA);
for (i = 0; i < idimA; i++)
  detA *= u[i * idim + i];
for (i = idimA; i < idim; i++)
  detB *= u[i * idim + i];
det = detA * detB;
if (det < 0 || isnan det) {
  it = -1;
  break;
}

// calculate a coefficient as (D-criterion)/(dimension of matrix)
coef = pow det, (1.0 / idim) / idim;
// inverse the sum matrix of weighted matrices
inverse (sumWeightedMatrices, idim, inverseSumWeightMat, idimA,
  detA, detB);
// multiply each element of the inverse matrix with the coefficient,
// the result obtained is the derivative of phi
MatrixTimesWeight (inverseSumWeightMat, &coef, isize, 1, Dphi);
// multiply every FIM calculated in R with the derivative of phi,
// then obtain the trace of each product
computeTrace (originMatrices, Dphi, idim, in, tr);
// evolve weights
sumw = 0.0; // sum of weights
sumwd = 0.0; // sum of lambda exponential weights
for (i = 0; i < in; i++)
{
  sumwd += weights[i] * tr[i];
  weights[i] *= pow (tr[i], *lambda);
  sumw += weights[i];
}
maxTrace = tr[0];
for (i = 0; i < in; i++)
{
  weights[i] /= sumw; // make sure weights are proper
fractions
    if (maxTrace < tr[i])
        maxTrace = tr[i];
}

if (maxTrace < (1 + *delta)*sumwd) // stop criterion of the
    algorithm
    {
        printf("loop finished with iteration = %d\n", it);
        break;
    }
}

*iteration = it;

free(weightedMatrices);
free(sumWeightedMatrices);
free(inverseSumWeightMat);
free(Dphi);
free(tr);
free(1);
free(u);
}
Appendix D: Fisher Information Matrix evaluation function

funFIMpop <- function (equation, parameters, beta, o, sigma, t_group, Trand, d, PropSubjects, nbTot) {
  #List of names of the fixed effects parameters and sampling times
  ParametersInEquation <- c (parameters, "t")
  PSI <- c (parameters, "sig.inter", "sig.slope")
  lpsi <- length (PSI)
  lengthParameters <- length (parameters)
  # (Diagonal Matrix of) variance for inter-subject random effects:
  omega <- diag (o)
  #Parameters of standard deviation model of residual error (sig.
inter+sig.slope*f)^2:
  sig.inter <- sigma [1]
  sig.slope <- sigma [2]
  #get derivatives of sigma
  Vmodel <- parse (text = paste0 ("( sig.inter + sig.slope * fixed ) ^2"))
  dv <- deriv (Vmodel [[1]], PSI)
  #gather all groups of protocol
  M_f <- list ()
  M_F <- matrix (rep (0), nrow=lpsi+lengthParameters, ncol=lpsi+
lengthParameters)
  for (q in 1: length (t_group)) {
    #design of sampling times
t <- t_group [[q]]
    #dose value
dose <- d[q]
    #calculate matrix E for n individuals
equatf <- parse (text = equation, n=-1)
    f <- function (ParametersInEquation) { eval (equatf [[1]]) }
```R
#Fixed effects parameters values
for (i in 1:lengthParameters){
  assign (parameters[i], beta[i])
}
```

```R
#calculate the observations with personal parameters
parameterValues <- c beta, t
fixed <- f (parameterValues)
```

```R
#calculate variance Var
var <- diag ((sig.inter + sig.slope * fixed) ^ 2)
```

```R
#get derivatives for fixed parameters
df <- deriv (equatf [[1]], PSI)
mdf <- attributes (eval(df)) $gradient
```

```R
#delete the last two columns (correspond to sig.inter and sig.slope)
mdfi <- mdf[,-c(lpsi-1, lpsi)]
```

```R
#complete derivative for random effect model
Random effect model (1) = additive (2) = exponential
```

```R
beta0 <- beta
```

```R
beta0[which (Trand==1)] <- 1
```

```R
mdfie <- mdfi %*% diag (beta0)
```

```R
#calculate variance Vi
Vi <- mdfie %*% omega %*% t (mdfie) + var
```

```R
#inverse of matrix Vi
SVi <- solve (Vi)
```

```R
#compute derivatives of sigma
mdv <- attributes (eval (dv)) $gradient
```

```R
#calculate matrix part A
MA <- t (mdfi) %*% SVi %*% mdfi
```
```R
#complete the rest of the matrix with 0
MA <- cbind(MA, matrix(rep(0), ncol=lpsi, nrow(lengthParameters))
)
MA <- rbind(MA, matrix(rep(0), ncol=lpsi+lengthParameters, nrow=lpsi))
```

```R
#calculate matrix part B
initialize the matrix with 0
```

```R
MB <- matrix(rep(0), nrow=lpsi+lengthParameters, ncol=lpsi+
lengthParameters)
#prepare a list of matrix of derivatives of sigma to simplify
usage
```

```R
if (length(t)==1){
  m <- lapply(c((lpsi-1), lpsi), function(i,mdv) mdv[i],mdv=mdv)
} else {
  m <- lapply(c((lpsi-1), lpsi), function(i,mdv) diag(mdv[,i]),
mdv=mdv)
}
```

```R
#calculate first three rows of part B
for(i in 1:lengthParameters){
  mdfiei <- (mdfie[,i] %*% t(mdfie[,i]))
```

```R
MB[ lengthParameters+i, seq(lengthParameters+1,
  lengthParameters+lengthParameters)] <- sapply(
  lapply(seq(1, lengthParameters), function(i,mdfie) mdfie[,i
  ],mdfie=mdfie),
  function(x) 1/2 * sum(diag((mdfiei) %*% SVi %*% (x %*% t(x)
  ) %*% SVi)))
```

```R
MB[ lengthParameters+i, c(lengthParameters+lpsi-1,
  lengthParameters+lpsi)] <- sapply(
  m, function(x) 1/2 * sum(diag((mdfiei) %*% SVi %*% x %*% SVi
  )))
}
#calculate the last two rows of partB
for(i in (lpsi-1):lpsi){
```
MB [lengthParameters +i, seq (lengthParameters +1,
lengthParameters +lengthParameters )] <- sapply(
  lapply (seq (1, lengthParameters ), function (i, mdfie) mdfie [,
  i
  ], mdfie =mdfie ),
  function (x) 1/2 * sum (diag ((m [[i − lengthParameters ]])) %*%
  SVi %*% (x %*% t (x)) %*% SVi ))

MB [lengthParameters +i, c (lengthParameters +lpsi − 1,
lengthParameters +lpsi )] <- sapply(
  m, function (x) 1/2 * sum (diag (((m [[i − lengthParameters ]]))
  %*% SVi %*% x %*% SVi )))
}

M.f [[q]] <- (MA+MB)*PropSubjects [q]
M.F <-M.F+M.f [[q]]
}

M.F <- M.F *nbTot

#set names for vectors
fname <-c (sapply (1: lengthParameters , function (x) paste0("u_",
  parameters [x] ))), sapply (
  1: lengthParameters , function (x) paste0("w2_", parameters [x] )),
  "sig .inter", "sig .slope")

rownames (M.F) <- fname
colnames (M.F) <- fname

if (0 %in% c (o, sigma)) {
  M.F<-M.F[-c (lengthParameters +which (c (o, sigma)==0)), -c (
  lengthParameters +which (c (o, sigma)==0))]
  PSI<- PSI[-c (lengthParameters +which (c (o, sigma)==0))]
}

if (length (t)==1) {
  return (list (M.F, det (M.F)))
} else {
  tryCatch (
    deterFim <- det (M.F)
SE <- sqrt(diag(solve(MF)))
RSE <- 100 * SE / c beta, o, sigma) [which(c beta, o, sigma) != 0)]
CritereDopt <- deterFim ^ (1 / length(c beta, o, sigma) [which(c beta, o, sigma) != 0)]))

\# write the output in console
\# cat("********************* FISHER INFORMATION MATRIX
***********************\n")
\# print(MF)
\# cat("\n\n********************** DETERMINANT OF THE
MATRI X ************************\n", deterFim,
\# "\n\n********************** CRITERION
*******************************\n", CritereDopt,
\# "\n\n********************** STANDARD ERROR
*******************************\n", SE,
\# "\n\n********************** RELATIVE STANDARD ERROR
***************************\n", RSE, "\n\n")
return(list(MF, deterFim, CritereDopt, SE, RSE))
}
error = function(e) {
return(list(MF, det(MF)))
}

\# draw the graphical illustration of model
\# t<-seq(min(unlist(t_group)),max(unlist(t_group)),(max(unlist(t_group))-min(unlist(t_group)))/100)
\# x<-c beta, t)
\# fvalue<-f(x)
\# lines(t, fvalue) #then plot the line of function
}

}
Appendix E: Application of multiplicative algorithm

1 model <- "exp(VA)+(1-exp(-exp(k)*t))* (emax*dose/(exp(ED50)+dose)-exp(be)*exp(VA))"

2 paramName <- c("VA", "k", "be", "emax", "ED50")

3 mu <- c(log(55), log(0.005), log(0.2), 30, log(150))

4 omega <- c(0.07, 0.5, 1, 150, 0)

5 sigma <- c(sqrt(28), 0)

6 allTimes <- c(0, 7, seq(28, 672, 28))

7 equiSpacedTimes <- c(0, 224, 420, 672)

8 fixed_times <- list(c(1,0),c(4,672))

9 Trand <- c(1,1,1,1,1)

10 dose <- c(0, 25, 50, 100, 150, 300, 500)

11 lambda <- 0.99

12 delta <- 0.00001

13 iteration <- 10000

14 nbSubjects <- 300

15

16

17 # Reference 26 times

18 start <- Sys.time()

19 ResultEval <- funFIMpop(model, paramName, mu, omega, sigma, list( allTimes, allTimes, allTimes, allTimes), Trand, c(0,150,300,500), c(75,75,75,75)/ nbSubjects, nbSubjects)

21 end <- Sys.time()

22 end-start

23

24 # Fix 26 times and optimize doses

25 # Create design space and elementary FIM corresponding

26 start <- Sys.time()

27 nbTimes <- 26

28 Designs_ini <- combn(allTimes, nbTimes)

29 FIMs<-list()

30 Designs<-matrix(nrow=nbTimes+1)

31 n<-1

32 for(i in 1:dim(Designs_ini)[2]) {
```R
for (j in 1:length(dose)) {
  # the number of design groups need to change in accordance
  with the number of doses
  FIMs[[n]] <- funFIMpop(model, paramName, mu, omega, sigma, list(
    Designs_init[, i]), Trand, dose[j], 1, 1)[[1]]
  Designs <- cbind(Designs, c(dose[j], Designs_init[, i]))
  n <- n + 1
}
Designs <- Designs[, -1]
end <- Sys.time()
end-start
```

```R
# Optimize using multiplicative algorithm implemented in R or C
start <- Sys.time()
weights <- multiplicative_algo(FIMs, "R", length(paramName),
    Designs, nbSubjects, lambda, delta, iteration)
end <- Sys.time()
end-start
```

```R
## equi-spaced times
start <- Sys.time()
resultEval <- funFIMpop(model, paramName, mu, omega, sigma, list(
  equiSpacedTimes, equiSpacedTimes,
  equiSpacedTimes, equiSpacedTimes),
  Trand, c(0,150,300,500), c(75,75,75,75)/
  nbSubjects, nbSubjects)
end <- Sys.time()
end-start
```

```R
# Reference 4 equi-spaced times
start <- Sys.time()
resultEval <- funFIMpop(model, paramName, mu, omega, sigma, list(
  equiSpacedTimes, equiSpacedTimes,
  equiSpacedTimes, equiSpacedTimes),
  Trand, c(0,150,300,500), c(75,75,75,75)/
  nbSubjects, nbSubjects)
end <- Sys.time()
end-start
```
```R
97 FIMs<-list()
98 Designs<-matrix(nrow=nbTimes+1)
99 o<-1
100 for(l in 1:dim(Designs_ini)[2]){
101   for(m in 1:length(dose)){
102     #the number of design groups need to change in accordance
103     with the number of doses
104     FIMs[[o]] <-funFIMpop(model,paramName,mu,omega,sigma list(
105                       Designs_ini[,1]),Trand,dose[m],1,1)[[1]]
106     Designs<-cbind(Designs,c(dose[m], Designs_ini[,1]))
107     o<-o+1
108   }
109 } Designs <- Designs[, -1]
110 end <- Sys.time()
111 end-start
112 # Optimize using multiplicative algorithm implemented in R or C
113 start <- Sys.time()
114 weights <- multiplicative_algo(FIMs, "R", length(paramName),
115                             Designs, nbSubjects, lambda, delta, iteration)
115 end <- Sys.time()
116 end-start
```
### 2.2.2. Contribution à la nouvelle version de PFIM

Le programme PFIM (www.pfim.biostat.fr) est conçu pour l'évaluation et l'optimisation de protocoles d'études longitudinales basées sur la FIM (individuelle, de population ou bayésienne). Sa version actuellement disponible en ligne est la version 4.0. (Dumont et al. 2018). Cette version, codée en R par plusieurs membres du groupe PFIM successifs, inclut beaucoup de fonctionnalités mais sa conception peut entraîner des difficultés à ajouter de nouvelles méthodes, voire à utiliser PFIM pour des modèles ou des protocoles complexes.

Une nouvelle version de PFIM, PFIM 5.0., ré-écrite entièrement, est donc en cours de développement. L'objectif de ce travail est d'obtenir un programme modulable, de manière à faciliter l'évolutivité et l'utilisation de PFIM. Pour cela, PFIM est écrit en R S4, langage orienté-objet. La programmation orientée-objet (POO) offre certains avantages par rapport à la programmation traditionnelle (procédurale). En concevant le programme en ensemble d'objets qui communiquent entre eux, la POO permet notamment de gagner en compréhensibilité, d'éviter les duplications de code et de faciliter la réutilisation.

PFIM 5.0. est donc organisé en différentes classes et différents objets. Comme PFIM 4.0., il permet l'évaluation ou l'optimisation de protocoles individuels, de population et bayésiens pour des modèles à une ou plusieurs réponses, en solution analytique ou sous forme d'équations différentielles ordinaires, provenant de la librairie de modèles ou définis par l'utilisateur. Certaines fonctionnalités qui n'étaient pas incluses dans PFIM 4.0. sont disponibles. Il est désormais possible d'utiliser des distributions des effets aléatoires différentes (normale et log-normale) au sein du même projet. L'algorithme multiplicatif (section 2.2.1.) est intégré à PFIM 5.0. en langage C (plus performant que R, ce qui facilitera la résolution des problèmes d'optimisation les plus complexes). L'utilisation de modèles PK de perfusion et des administrations multiples a été simplifiée par rapport à PFIM 4.0. Les graphiques présentés ont aussi été revus, avec utilisation de la librairie R *ggplot2*. Certaines fonctionnalités présentes dans PFIM 4.0. doivent encore être ajoutées, par exemple, l'optimisation de protocoles par l'algorithme de Fedorov-Wynn et l'inclusion des effets de covariables catégorielles. La nouvelle version de PFIM ainsi développée sera prochainement disponible sur le CRAN (*Comprehensive R Archive Network*. Des mises à jour futures pourront aussi permettre d'inclure les récents développements en optimisation de protocoles tels que l'évaluation de la FIM par des méthodes adaptées aux modèles de données discrètes ou encore des méthodes de planification tenant compte de l'incertitude des paramètres et du modèle (section 2.1.).
En tant que membre actif du groupe PFIM, mes rôles incluent le traitement des requêtes d'utilisateurs de PFIM 4.0. (version R et interface). Concernant l'élaboration de PFIM 5.0., mes activités principales sont la participation aux réunions hebdomadaires et l'élaboration d'exemples dans le but de vérifier l'ajout de nouvelles méthodes au logiciel ou pour être intégrés au futur guide d'utilisateur. Un poster à propos de la nouvelle version de PFIM a aussi été produit (voir Annexes) et notamment présenté au PAGE (Population Approach Group of Europe), un congrès de pharmacométrie (Seurat et al. 2019).

Une version alpha de PFIM 5.0. est actuellement disponible en interne, une version beta sera disponible pour un nombre plus large d'utilisateurs durant le premier semestre de l'année 2021, dans le but de soumettre une version finale de PFIM 5.0. au cours du deuxième semestre de l'année 2021.

## 2.3. Comparaison de protocoles pour évaluer la décroissance tumorale précoce

### 2.3.1. Résumé

Le choix de protocole est complexe dans les études cliniques en oncologie, particulièrement dans les phases cliniques précoces du développement des anticancéreux. Dans ces études, il est de plus en plus recommandé de comprendre la relation entre la concentration des médicaments et la réponse clinique (Smoragiewicz et al. 2020). Cette dernière peut être suivie par des biomarqueurs tels que la décroissance tumorale précoce, c'est-à-dire l'effet des premières semaines de thérapie sur la taille tumorale. Nous avons donc évalué et comparé différents protocoles à partir de simulations, dans le cadre de l'étude d'une combinaison de deux anticorps monoclonaux (associé au traitement standard, la chimiothérapie) pour le traitement du cancer colorectal en phase précoce clinique.

Des jeux de données de tailles tumorales (la taille tumorale est en fait la somme des plus longs diamètres de chaque lésion) ont été simulés à partir de modèles d'inhibition de la croissance tumorale, proches du modèle de Claret (Claret et al. 2009). Trois différents scénarios ont été étudiés concernant l'interaction entre les deux médicaments étudiés : le cetuximab (C, traitement de référence) et une nouvelle molécule (M). Dans le premier scenario, nous avons considéré un effet nul de M dans le but d'évaluer l'erreur de type 1. Dans les deux autres scenarios, nous avons considéré respectivement une simple additivité ou une synergie d'effet entre C et M. Pour chacun de ces scenarios, trois protocoles d'études ont été évalués : à un bras,
deux bras et quatre bras de traitement, avec un même nombre total d'individus (N=60). Pour chaque protocole et chaque scenario, 500 jeux de données ont été générés et des mesures de taille tumorale à l'initiation de traitement (semaine 0) puis toutes les deux semaines jusqu'à la semaine 8. Ces jeux de données ont ensuite été analysés par modèles non-linéaires à effets mixtes, avec une estimation des paramètres par l'algorithme SAEM de Monolix 2018R2. Le critère d'évaluation de la combinaison était fondé sur la décroissance de taille tumorale à 8 semaines (ETS8), c'est-à-dire l'écart relatif entre la taille tumorale à la semaine 8 par rapport à celle de la semaine 0. Les tests statistiques, présentés dans le paragraphe suivant, ont été réalisés pour comparer, entre les bras, les distributions d'ETS8, soit observées, soit prédites à partir des *Empirical Bayes estimates* (EBE) obtenus par modélisation.

Pour le design à un bras, les 60 individus ont été traités par la combinaison C+M aux doses standard, les ETS8 de ces individus ont été comparées par test de Wilcoxon à une valeur de référence obtenue lors d'une étude antérieure. Pour le design à deux bras, 30 individus ont été traités par la combinaison C+M aux doses standard et 30 par le traitement de référence (C). Les ETS8 ont été comparés entre les deux groupes par test de Wilcoxon. Pour le design à quatre bras (15 individus par bras), en plus des deux groupes C et C+M aux doses standards, des groupes de combinaisons C+M avec des doses inférieures de C ont été étudiés, les ETS8 ont été comparés entre les groupes de manière globale par test de Kruskal-Wallis, puis les bras de combinaison C+M ont été comparés au bras de référence C par test de Dunnett.

Pour un nombre identique d'individus, le protocole à un bras présente l'avantage d'obtenir une meilleure puissance statistique qu'un autre protocole, à condition de se comparer à la bonne valeur de référence. Par exemple, avec les ETS8 observés et le scénario d'additivité, la puissance est de 89 % avec ce protocole contre 67 % pour un design à deux bras. Néanmoins, la valeur historique à laquelle on compare la distribution d'ETS8 peut être incorrecte parce que les conditions d'études ou les critères d'inclusion sont différents de l'étude qui a permis d'obtenir la valeur historique. Cela peut conduire à une inflation importante de l'erreur de type 1. Par exemple, si la valeur de référence d'ETS8 choisie était 20 % alors que la valeur correcte était 29 %, cela conduisait à une erreur de type 1 de 57 % (hors de l'intervalle de prédiction à 95 % [3,3 ; 7,3] %). De plus, si la valeur historique est plus élevée que la valeur réelle, cela pouvait conduire à une perte de puissance, par exemple si la valeur de référence choisie était 35 % au lieu de 29 %, la puissance était seulement de 45 % (vs. 89 % avec la valeur correcte). L'essai randomisé à deux bras présente l'avantage important de ne pas avoir à faire d'hypothèse sur la valeur de référence lors de l'analyse des données. Ce design permet de contrôler l'erreur
de type 1 qui était de 6 % dans notre exemple. De plus, dans le scenario de synergie, ce protocole permettait d'obtenir une puissance de 97 % (vs. 100 % pour un design à un bras). Enfin, le protocole à 4 bras peut être intéressant lorsque l'on recherche à sélectionner la dose la plus adéquate, par exemple pour des raisons de toxicité. En revanche, étant donné le plus faible nombre d'individus par bras avec ce protocole et la procédure statistique plus conservatrice, dans le scénario d'additivité, les puissances pour détecter un effet de la combinaison avec la dose forte de C (standard), moyenne et faible étaient respectivement de 20, 13 et 3 %. Dans le cas d'une synergie, les puissances respectives atteignaient 60, 52 et 26 %.

Par ailleurs, certains éléments de ce travail ont montré l'intérêt des modèles mixtes pour l'analyse des données longitudinales de tailles tumorales. En effet, ce type de données est sujet à une erreur de mesure pouvant être importante. Les modèles permettent de prendre en compte ces erreurs de mesure, et ceci permet de réduire la perte de puissance qui leur sont dues. Nous avons observé cette tendance, par exemple, pour le design à quatre bras et le scenario de synergie, la puissance était de 65 % à partir des prédictions basées sur le modèle contre 60 % à partir des données observées.

En conclusion, ce travail a notamment permis de mettre en évidence les avantages et inconvénients de différents designs pour évaluer l'ETS pour un contexte de combinaison de molécule dont l'une est le traitement de référence. Il a fait l'objet d'un article publié dans *Clinical Pharmacology and Therapeutics: Pharmacometrics & Systems Pharmacology*.

### 2.3.2. Article 3 (publié)

*Comparison of various phase I combination therapy designs in oncology for evaluation of early tumor shrinkage using simulations*, publié dans *Clinical Pharmacology and Therapeutics: Pharmacometrics & Systems Pharmacology*, 2020

Jérémy Seurat, Pascal Girard, Kosalaram Goteti, France Mentré
<figure><img src="image_17.png" /></figure>

ARTICLE

# Comparison of Various Phase I Combination Therapy Designs in Oncology for Evaluation of Early Tumor Shrinkage Using Simulations

Jérémy Seurat<fnref n="1" />*, Pascal Girard², Kosalaram Goteti³ and France Mentré<fnref n="1" />

There is still a lack of efficient designs for identifying the dose response in oncology combination therapies in early clinical trials. The concentration response relationship can be identified using the early tumor shrinkage time course, which has been shown to be a good early response marker of clinical efficacy. The performance of various designs using an exposure-tumor growth inhibition model was explored using simulations. Different combination effects of new drug M and cetuximab (reference therapy) were explored first assuming no effect of M on cetuximab (to investigate the type I error (α)), and subsequently assuming additivity or synergy between cetuximab and M. One-arm, two-arm, and four-arm designs were evaluated. In the one-arm design, 60 patients received cetuximab + M. In the two-arm design, 30 patients received cetuximab and 30 received cetuximab + M. In the four-arm design, in addition to cetuximab and cetuximab + M as standard doses, combination arms with lower doses of cetuximab were evaluated (15 patients/arm). Model-based predictions or "simulated observations" of early tumor shrinkage at week 8 (ETS8) were compared between the different arms. With the same number of individuals, the one-arm design showed better statistical power than other designs but led to strong inflation of α in case of misestimated reference for ETS8 value. The two-arm design protected against this misestimation and, with the same total number of subjects, would provide higher statistical power than a four-arm design. However, a four-arm design would be helpful for exploring more doses of cetuximab in combination with M to better understand the interaction.

## Study Highlights

**WHAT IS THE CURRENT KNOWLEDGE ON THE TOPIC?**

☑ Currently, there is no standard for design choice in early clinical studies of combinations treatments in immuno-oncology.

**WHAT QUESTION DID THIS STUDY ADDRESS?**

☑ These clinical trial simulations compared designs for the first time in the context of a fixed dose of one agent assumed to be optimal while testing different doses of another based on early tumor shrinkage.

**WHAT DOES THIS STUDY ADD TO OUR KNOWLEDGE?**

☑ Statistical tests performed on individual model-based early tumor shrinkage at week 8 values and the one-arm design indicate a better power than two-arm or four-arm designs, but imply strong assumptions about the historical reference value, leading to strong inflation of type I error in the case of underestimating the reference. Choosing a

two-arm or a four-arm design depends on the objective of the study: a two-arm design is preferable to a four-arm design to achieve a good statistical power, but a four-arm design allows better exploration of the combination and better dose selection.

**HOW MIGHT THIS CHANGE DRUG DISCOVERY, DEVELOPMENT, AND/OR THERAPEUTICS?**

☑ These simulations show that, in early-phase oncology combination studies (with early tumor shrinkage as the outcome), the choice of a one-arm design (i.e., without a comparator arm) is not appropriate unless the conditions of the clinical trial are similar to those of the reference treatment or because of ethical reasons. This implies strong assumptions about the historical value, and the trial could lead to a wrong conclusion if the historical value is no longer relevant.

Early-phase studies are critical in oncology drug development. These trials are frequently designed to evaluate the overall response as well as toxicity and pharmacokinetics. Go/no go decisions and dose recommendations for further studies are based on the results of phase I studies.

Therefore, the design choice is crucial in the early phase.<fnref n="1" /> There is growing interest in the use of combination therapies in oncology, especially immuno-oncology.² The combination of several agents could have a synergistic effect and enhance the antitumor activity. However, toxicities are

<footnote n="1">¹Université de Paris, INSERM, IAME, F-75006 Paris, France; ²Merck Institute for Pharmacometrics, Merck Serono S.A, Lausanne, Switzerland; ³EMD Serono R&D Institute, Billerica, Massachusetts, USA. *Correspondence: Jérémy Seurat (jeremy.seurat@inserm.fr) Received: June 26, 2020; accepted: September 21, 2020. doi:10.1002/psp4.12564</footnote>
frequently observed in these studies, especially when chemotherapies are involved.³ In the past two decades, several new targeted therapies, such as monoclonal antibodies, with efficacy in cancer subpopulations and limited adverse effects compared with traditional therapies have reached the market.⁴ These drugs have recently been tested in many combination trials in addition of chemotherapy/radiotherapy.⁵,⁶ However, the designs of clinical trials investigating combinations of immune-oncology agents are often based on empiricism without clear assumptions about the drug interactions. Therefore, the design choice in combination trials could be challenging, especially when clinical information is only available for single agent given as monotherapy. Recommendations have been made concerning the design of early studies with a combination of targeted anticancer therapies.⁷ The importance of the assumed pharmacodynamic drug interaction and the usefulness of biomarkers are highlighted in these recommendations. More recent recommendations insist on the importance of including drug concentrations in the analysis and understanding the pharmacokinetic/pharmacodynamic relationships.⁸ The importance of including biomarker-driven objectives in phase I trials has also been shown.⁹

Among clinical trial end points widely used in oncology, overall survival is universally accepted but involves large studies with long follow-up, which is not compatible with early-phase studies. End points based on tumor assessment can be used. The Response Evaluation Criteria in Solid Tumors (RECIST) criteria, widely adopted in oncology clinical trials including on targeted therapies, consist of four categories that depend on the evolution of tumor size (TS).¹⁰,¹¹ However, in the field of modeling and simulation it seems preferable to use the continuous TS data when available rather than RECIST as this would result in a loss of information.¹²,¹³ Furthermore, the link between the early tumor shrinkage (ETS) and overall survival or progression free-survival has been proposed as an early efficacy marker based on several studies.¹⁴⁻¹⁶

Nonlinear mixed effect models (NLMEM) are increasingly used to support drug decision making, especially in early-phase trials.¹⁷ In oncology, these models can be used to capture the pharmacokinetics of antitumoral agents or biomarker kinetics.¹⁸ NLMEM can describe tumor dynamics, and several tumor growth models are found in the literature such as Gompertz¹⁹ or Simeoni.²⁰ Tumor growth inhibition (TGI) models including a drug induced decrease in TS are used for various types of cancer and antitumoral agents, including targeted therapies.¹⁵,²¹ TS (or other biomarkers) can also be linked to survival by a joint model.²² However, before conducting a trial and modeling the data, it is crucial to choose an appropriate design to obtain good precision of the parameter estimates and thus precise results. For this purpose, model-based clinical trial simulation (CTS) can be performed to evaluate and compare designs.

In simulations, we used a combination of candidate drug M given on top of cetuximab, a monoclonal antibody targeting the epidermal growth factor receptor of tumor cells²³ that has been approved for metastatic colorectal cancer with wild-type RAS and usually given in combination with

chemotherapy. Several ongoing or completed studies indicate that new combinations including cetuximab could be promising for different cancers.²⁴,²⁵ The aim of this study was to compare by CTS several combination designs in the treatment of refractory colorectal cancer with different assumptions about the interaction between cetuximab and M. The performances, type I error (α), and power of several designs were compared to test the superiority of the combination treatment to cetuximab alone based on ETS at week 8 (ETS8). Another objective was to compare, by modeling and simulation of exposure-TGI, the test performances on predicted ETS8 obtained from the model to the observed ETS8.

## METHODS

The different models used to perform the simulations were the model of pharmacokinetic exposure during the first 8 weeks of therapy for each drug and the TGI model. We studied three different designs (one arm, two arms, and four arms) in three different scenarios: no effect of M, additivity, and synergy. Individual ETS8 values were computed as predictions after global model fitting, observations, or true values.

### Models

**Cetuximab and drug M pharmacokinetic exposures.** Cetuximab has been described by various pharmacokinetic models.²⁶,²⁷ A significant influence of body surface area (BSA) on clearance of cetuximab has been shown.²⁸ Therefore, it is recommended to adjust the dose of cetuximab (Dose<sub>C</sub>) on BSA.²⁹ We assumed a simple model with independent clearance with dose, as found in a previous study in patients with metastatic colorectal cancer.³⁰ Therefore, cetuximab exposure expressed as area under the curve (AUC<sub>C</sub>) integrated over 2 weeks (between two administrations), in mg × week/L, was simulated as in Eq. 1.

$$AUC_{C:2 weeks} = \frac{Dose_{E}(mg/m^{2}) \times BSA(m^{2})}{CL_{C}(L/week) \times (BSA/1.85)^{0.75}}, \quad (1)$$

where CL<sub>C</sub> is the cetuximab clearance value (CL<sub>C</sub>) presented in Table 1, including BSA effect standardized by 1.85 (i.e., the median BSA that was expected before the clinical trial), and BSA was simulated according to a lognormal distribution with a median of 1.75 and m² and an interindividual variability of 10% according to derived BSA distribution in previous cetuximab studies.³¹⁻³³

The current maximal dosing regimen for cetuximab given in combination with another treatment is 500 mg per m² every 2 weeks. This dosing regimen has a comparable activity and safety to the standard dosing regimen.³⁴ A dosing regimen of 400 and 200 mg per m² every 2 weeks were also considered in case of combination with drug M.

Concerning the drug M, its area under the curve (AUC<sub>M</sub>) was integrated over 2 weeks (between two administrations) and simulated as in Eq. 2:

$$AUC_{M:2 weeks} = \frac{Dose_{M}}{CL_{M}}, \quad (2)$$
<table><caption>Table 1 Parameter values for area under the curve models of cetuximab and drug M and for the tumor growth inhibition model</caption><thead><tr><th>Parameter (unit)</th><th>μ</th><th>ω<sup>2</sup></th></tr></thead><tbody><tr><td>Pharmacokinetics models</td><td></td><td></td></tr><tr><td>CL<sub>C</sub> (L week<sup>-1</sup>)</td><td>3.9</td><td>0.0025</td></tr><tr><td>CL<sub>M</sub> (L week<sup>-1</sup>)</td><td>5.0</td><td>0.01</td></tr><tr><td>Tumor growth inhibition model</td><td></td><td></td></tr><tr><td>TS0 (mm)</td><td>100<sup>a</sup></td><td>0.5<sup>a</sup></td></tr><tr><td>KS (week<sup>-1</sup>)</td><td>0.001</td><td>—</td></tr><tr><td>KD<sub>50C</sub> (week<sup>-1</sup>)</td><td>0.015</td><td>1.5</td></tr><tr><td>KD<sub>C</sub> (L × mg<sup>-1</sup>week<sup>-2</sup>)</td><td>0.00025<sup>a</sup></td><td>1.5<sup>a</sup></td></tr><tr><td>KD<sub>M</sub> (L × mg<sup>-1</sup>week<sup>-2</sup>)</td><td>0.00025<sup>a</sup></td><td>1<sup>a</sup></td></tr><tr><td>KR (week<sup>-1</sup>)</td><td>0.2</td><td>1</td></tr><tr><td>Interaction parameters</td><td>Additivity</td><td>Synergy</td><td>Additivity</td><td>Synergy</td></tr><tr><td>INT</td><td>0</td><td>2<sup>a</sup></td><td>—</td><td>0.16<sup>a</sup></td></tr><tr><td>INT<sub>50</sub> (mg × week/L)</td><td>—</td><td>75</td><td>—</td><td>—</td></tr><tr><td>Correlations</td><td></td><td></td><td></td><td></td></tr><tr><td>Cov(η<sub>KR</sub>, η<sub>KDC</sub>)</td><td></td><td>1<sup>a</sup> (ρ = 0.82)</td><td></td><td></td></tr><tr><td>Residual error</td><td></td><td></td><td></td><td></td></tr><tr><td>σ<sup>2</sup><sub>prop</sub></td><td></td><td>0.023<sup>a</sup></td><td></td><td></td></tr></tbody></table>

μ are the fixed effects, ω² the variance of random effects η, and σ² the variance of residuals.

CL<sub>C</sub>, cetuximab clearance value; CL<sub>M</sub>, clearance of drug M; Cov, covariance; INT, interaction parameter; INT<sub>50</sub>, exposure needed for 50% of interaction maximal effect; KD<sub>C</sub>, effect of drug cetuximab; KD<sub>M</sub>, effect of drug M; KD<sub>50C</sub>, standard of care effect; KR, constant rate of resistance development; KS, natural growth of TS; L, week<sup>-1</sup>, Litre, week<sup>-1</sup>; η, random effects; TS0, tumor size at baseline.

<sup>a</sup>Estimated parameter.

where CL<sub>M</sub> is the clearance of drug M (CL<sub>M</sub>) presented in Table 1.

**Tumor size.** Tumor size is defined as the sum of the longest diameters of target lesions that are simulated through a TGI model.¹⁵ The concentration of anticancer antibodies such as cetuximab is associated with the response rate.²⁷ A previous TGI model was developed based on clinical data from patients affected by metastatic colorectal cancer.²⁸ In this analysis, the TS dynamics was described by the ordinary differential equation (Eq. 3) linking TS and both drug exposures (AUC<sub>C</sub> and AUC<sub>M</sub>).

$$ \frac{d\text{TS}}{dt} = \text{KS} \times \text{TS} - \text{KD} \times e^{-\text{KR} \times t} \times \text{TS} $$

$$ \text{KD} = \text{KD}_{\text{50C}} + \text{KD}_{\text{C}} \times \text{AUC}_{\text{C}} + \text{KD}_{\text{M}} \times \text{AUC}_{\text{M}} \times \left( 1 + \frac{\text{INT} \times \text{AUC}_{\text{C}}}{\text{INT}_{\text{50}} + \text{AUC}_{\text{C}}} \right) (3) $$

$$ \text{TS}(0) = \text{TS0} $$

In this model (see Supplementary file), KS is the natural growth of TS, KD includes the different pooled drug effects on tumor shrinkage, KR represents the constant rate of resistance development, and TS0 is the TS at baseline. KD is the sum of KD<sub>50C</sub> (for standard of care, chemotherapy), KD<sub>C</sub> × AUC<sub>C</sub> (for cetuximab, with KD<sub>C</sub> the effect of cetuximab), and KD<sub>M</sub> × AUC<sub>M</sub> × (1 + INT×AUC<sub>C</sub>/(INT<sub>50</sub>+AUC<sub>C</sub>)) (for drug M, with KD<sub>M</sub> the effect of M).

According to Eq. 3, the effect of drug M is the sum of its proper effect and its interaction with cetuximab. This interaction model is inspired by the global pharmacodynamics interaction model.³⁵ Depending on the sign of the interaction parameter (INT), this model is able to capture different types of interactions. A positive INT represents synergism, a negative INT represents antagonism, and if this parameter is null, the effects of both drugs are additive. Parameter

values of this model are given in Table 1 based on previous reported analysis on cetuximab (excepted KD<sub>M</sub> and interaction parameters) in the treatment of colorectal cancer.²⁸ The standard of care effect (KD<sub>50C</sub>) includes the placebo response. KD<sub>M</sub> is selected from a sensitivity analysis to obtain plausible tumor shrinkage. Two interaction scenarios are investigated: additivity in which the two compounds act independently and synergy in which the cetuximab exposure increases the drug M effect. The difference between these two scenarios in terms of median TS evolution from 0 to 8 weeks is illustrated in Figure S1.

Let *i* denote the *i*-th individual (*i* = 1, ..., *N*) and *j* the *j*-th TS measurement of an individual. The statistical model for the TS observation TS<sub>ij</sub> in individual *i* at time t<sub>ij</sub> is given by Eq. 4.

$$ \text{TS}_{ij} = f(t_{ij}, \theta_i) + \varepsilon_{ij}, \quad (4) $$

where *f* is the TGI model described in Eq. 3, θ<sub>j</sub> is the vector of individual parameters, and ε<sub>ij</sub> is the residual errors of mean 0 and of variance proportional to the TS: Var(ε<sub>ij</sub>) = *f*<sup>2</sup>(t<sub>ij</sub>, θ<sub>j</sub>) × σ<sup>2</sup><sub>prop</sub>. The random effect model is exponential as θ<sub>j</sub> = μ × exp(η<sub>j</sub>), where μ is the fixed effect and η<sub>j</sub> the random effect, except for INT in which the random effect model is additive: θ<sub>j</sub> = μ + η<sub>j</sub>. The AUCs of both cetuximab and M are simulated assuming a constant clearance and used as regressors in the TGI model.

Simulation settings

**Designs.** We studied three types of design—one arm, two arms, or four arms—considering the same total number of individuals (*N* = 60). Simulated TS are performed at predose (baseline) and weeks 2, 4, 6, and 8. The inclusion criterion was a baseline TS > 20 mm. Dropout, incomplete adherence, missing TS measurements, and protocol
deviations were not considered. The therapy for all patients included a standard of care, which is omitted next for simplification. Simulations were performed using R software (R Foundation for Statistical Computing, Vienna, Austria) version 3.4.3., specifically, packages mvtnorm to generate individual parameters³⁶ and deSolve to simulate longitudinal data.³⁷ Simulated TS < 10 mm were censored and considered to be below the limit of quantification. The dosing regimens of cetuximab and drug M for each possible design are described in the following paragraphs. The dose of M was assumed to be fixed. For each design and each interaction scenario, 500 data sets were simulated.

**One-arm: M + C500.** In the one-arm design, the N = 60 patients were allocated to a combination M + C500: drug M and cetuximab. The dosing regimen was fixed for each drug: cetuximab at 500 mg/m² every 2 weeks (C500) and drug M at its standard dose.

**Two-arm: C500/M + C500.** In the two-arm design, the N = 60 patients were randomized among the two possible arms: cetuximab “alone” (C500) or M + C500, resulting in 30 patients in each arm. Typical simulated data for patients in the two arms for the additivity and synergy scenarios are presented in **Figure S2**.

**Four-arm: C500/M ± C200/M ± C400/M + C500.** In the four-arm design, two combination arms with lower doses of cetuximab were considered (M + C400 and M + C200) in addition to C500 and M + C500. In the lower dose arms, the dose of C was 400 mg/m² (C400) or 200 mg/m² every 2 weeks. Each arm included 15 patients.

**Scenarios of the effect of M + C.** Three scenarios concerning the effect of M and the interaction between M and C were considered. To investigate the type I error of tests (superiority of M + cetuximab vs. cetuximab alone), a first scenario without an effect of drug M (i.e., with KD<sub>M</sub> = 0) was studied. In the other two scenarios, an effect of M (μ<sub>KDM</sub> = 0.00025) was implemented. These two scenarios differ in the interaction between M and cetuximab: additivity (INT = 0) or synergy (μ<sub>INT</sub> = 2).

## Analysis

**Predictions, observations, and true individual ETS8 values.** For each individual *i*, ETS8 is expressed as a percentage and calculated as in Eq. 5.

$$ETS8_i = \frac{(TS_i(0) - TS_i(8))}{TS_i(0)} \times 100, \quad (5)$$

where TS<sub>i</sub>(0) and TS<sub>i</sub>(8) are individual TS at baseline and after 8 weeks of treatment, respectively. Predicted (IPRED), observed (OBS), and true (TRUE) ETS8 are calculated using predictions, observations, and true individual TS at both week 0 and week 8, respectively, as described in the next paragraphs.

Individual predictions (IPRED) were obtained from model-based parameters. Parameters with a superscript letter in **Table 1** are estimated based on the TGI model (Eq. 1) by fitting

the simulated data using the SAEM (Stochastic Approximation Expectation-Maximization) algorithm³⁸ in MONOLIX 2018R2 software (Lixoft, Antony, France). MONOLIX runs were performed with five Markov chains, a minimum of 500 iterations during the exploratory phase and a minimum of 200 iterations during the smoothing phase. An adequate method is used to account for below the limit of quantification data.³⁹ TS individual predictions at baseline and week 8 were obtained using the mode of empirical Bayes estimates of the individual parameters. In the synergy model, the μ<sub>INT</sub> and its interindividual variability ω<sub>INT</sub> were estimated, whereas in the additivity model they were fixed to 0. Simulated data without an effect of M were refitted for both the additivity and synergy models to evaluate the type I error for each assumed interaction.

Individual simulated TS observations (OBS) were directly calculated from the simulated TS values plus residual error while true individual TSs (TRUE) were simply derived from individual parameters without residual error. It is important to realize that true values are not observed in a real clinical trial. An example of OBS, IPRED, and TRUE TSs in one individual is illustrated in **Figure S3**.

**Design evaluation.** The performances of each design were evaluated and compared in terms of type I error and power to detect the superiority of the combination (drug M + cetuximab) to cetuximab based on the ETS8 values (OBS, IPRED, and TRUE). The statistical tests described in the next paragraphs were two-sided using α = 5%.

For the one-arm design, as all the individuals were allocated to the combination arm, individual ETS8 values were compared with a reference ETS8 for cetuximab alone from historical data using a one-sample Wilcoxon test. Different scenarios were investigated concerning the reference in which the historical value is either not valid (lower or higher than true, due, e.g., to standard of care changes or differences in patient cohorts) or true for this reference ETS8. In this study, the reference value was 29%, (i.e., the median of all individual simulated ETS8 values with C500).

For the randomized trial with two arms, C500 and cetuximab + M, a two-sample Wilcoxon test was performed to compare the ETS8 values of individuals from each arm.

In the four-arm design, a global Kruskal-Wallis test was performed first to compare the ETS8 values between the four arms. If a difference between at least two arms was detected by this global test, a Dunnett test was performed to compare each combination arm to the reference C500 alone using a Benjamini-Hochberg procedure to account for multiplicity of tests.⁴⁰

From the data sets simulated without an effect of drug M (KD<sub>M</sub> = 0, i.e., under H<sub>0</sub>), the type I error was evaluated as the proportion of trials for which H<sub>0</sub> is rejected and compared to the nominal value of 5%. The 95% prediction interval for proportion π of H<sub>0</sub> rejection (type I error or power) is given by binomial law<sup>41</sup> ((3.3%, 7.3%) with 500 data sets). From the data sets simulated with an effect of drug M (additivity or synergy scenario, i.e., under H<sub>1</sub>), the power of the tests was evaluated as the proportion of trials for which H<sub>0</sub> is rejected.
## RESULTS

First, we checked whether model-based ETS8 values (IPRED) were overpredicting or underpredicting the true values. Overall, ETS8 are distributed accurately around the identity line for each design and under both additivity and synergy scenarios. However, a trend of overpredicted values for ETS8 < 25% is shown, although it seems to not be influenced by the scenario, design, or treatment arm (Figure S4).

### One-arm: M + C500

For the one-arm design, the ETS8 values for patients were compared to an ETS8 value for current reference therapy obtained from historical data. Distributions of individual ETS8 (OBS, IPRED, and TRUE) are presented in **Figure 1** in the cases of additivity and synergy. The median of OBS, IPRED, and TRUE ETS8 values were close to each other in each scenario. The observed median ETS8 was 43% for additivity and 55% for synergy, whereas the reference ETS8 for cetuximab alone was 29%. According to our model and simulation settings, assuming additivity, almost three of four patients had greater tumor shrinkage at week 8 than those of reference treatment. These boxplots enable us to see the interindividual variability in ETS8. As predictions account for residual error, the 5th to 95th percentile of ETS8 is narrower for IPRED (from 23% to 78% for additivity and from 28% to 92% for synergy) than OBS (from 3% to 78% for additivity and from 13% to 89% for synergy). In both scenarios, the distribution of predictions seems to be more representative of true values than observations.

With this one-arm design, as the ETS8 distribution is only known for M + C500, a one-sample Wilcoxon test was performed (Table 2). Using observed data and assuming a reliable ETS8 reference value under C500 (in our case 29%), the type I error is controlled (5.6%) as it lies in the prediction interval around the nominal value $\alpha = 5\%$. However, the type I error is not controlled using model-based predictions (12.6% for additivity and 13.2% for synergy), which is partly attributed to a lack of identifiability between the M and cetuximab effect. Overall, assuming a reliable reference, good power (at least 89%) is reached using either OBS or IPRED.

In this kind of design, for multiple possible reasons, the assumed reference value from historical data could be

wrong. Therefore, we investigated the type I error and power with lower (25% and 20%) and higher (35% and 40%) ETS8 reference values. Performing a two-sided Wilcoxon test, the type I error was inflated in cases using either OBS or TRUE. For example, if we compare the ETS8 distribution to 20% instead of 29%, we wrongly conclude to superiority of the combination in more than the half of the CTS using observations or predictions. Moreover, assuming a higher reference value than true leads to a lack of power. For example, in the additivity scenario, assuming a reference value of 35% leads to a power of 45.2% using observations (instead of 89.0% with the real reference value).

### Two-arms: C500/M + C500

For the randomized two-arm design, the distribution of individual ETS8 (OBS, IPRED, and TRUE) values in each arm is shown in **Figure 2** for additivity and synergy. As with the one-arm design, the median OBS, IPRED, and TRUE ETS8 values are close to each other in each arm and each scenario. For C500, the median ETS8 is ~ 29% (i.e., the reference value in the one-arm design) and for M + C500, the median ETS8 is ~ 43% with additivity and ~ 55% with synergy. Notably, with C500, more than 1 of 10 patients presented a negative observed ETS8 (i.e., increasing TS between weeks 0 and 8), whereas none of the individuals had a true negative ETS8 under this therapy with our simulation settings, which is captured well using predicted ETS8. The IPRED distribution of C500 seems to be similar between the two scenarios. This indicates that the estimated cetuximab specific effect is not influenced by the model or the scenario, as in the simulated data, and allows for comparison between the scenarios. As with the individuals given M + C500 in the one-arm design, the variability in ETS8 is considerable with C500 in the two-arm designs, especially in OBS (the 5th to 95th percentile extends from -12% to 65%). Therefore, the difference in ETS8 distribution between the two arms is difficult to distinguish in the case of additivity. The IPRED distribution is narrower and closer to TRUE distribution than OBS.

With this design, comparing ETS8 distributions between the two arms by performing two-sample Wilcoxon tests, no assumption on ETS8 reference value is required (Table 3). Therefore, choosing a two-arm design instead of a one-arm design avoids a potential inflation of the type I error. In both

<figure><img src="image_18.png" /><figcaption>**Figure 1** Boxplots of early tumor shrinkage at week 8 (ETS8) from observations (OBS; solid boxes, left subgroup for each treatment), predictions (IPRED; dotted boxes, middle subgroup for each treatment), and true values (TRUE; dashed boxes, right subgroup for each treatment) in the one-arm design (combination M + C500). The blue horizontal line is the true median ETS8 (29%) of the historical reference treatment. Whiskers indicate 5th to 95th percentiles.</figcaption></figure>
<table><caption>Table 2 Type I error and power of tests in one-arm design</caption><thead><tr><th rowspan="3">Theoretical median of ETS8 from historical data (%)</th><th colspan="6">Additivity</th><th colspan="6">Synergy</th></tr><tr><th colspan="3">Type I error (%)</th><th colspan="3">Power (%)</th><th colspan="3">Type I error (%)</th><th colspan="3">Power (%)</th></tr><tr><th>OBS</th><th>IPRED</th><th>TRUE</th><th>OBS</th><th>IPRED</th><th>TRUE</th><th>OBS</th><th>IPRED</th><th>TRUE</th><th>OBS</th><th>IPRED</th><th>TRUE</th></tr></thead><tbody><tr><td>Lower than true</td><td>20</td><td>56.6</td><td>79.4</td><td>68.6</td><td>100.0</td><td>100.0</td><td>100.0</td><td>56.6</td><td>60.2</td><td>68.6</td><td>100.0</td><td>100.0</td><td>100.0</td></tr><tr><td></td><td>25</td><td>13.2</td><td>34.2</td><td>18.0</td><td>98.6</td><td>99.6</td><td>99.8</td><td>13.2</td><td>34.2</td><td>18.0</td><td>100.0</td><td>100.0</td><td>100.0</td></tr><tr><td>True reference</td><td>29</td><td>5.6</td><td>12.6</td><td>6.0</td><td>89.0</td><td>97.6</td><td>93.6</td><td>5.6</td><td>13.2</td><td>6.0</td><td>100.0</td><td>100.0</td><td>100.0</td></tr><tr><td>Higher than true</td><td>35</td><td>27.4</td><td>31.6</td><td>37.4</td><td>45.2</td><td>70.6</td><td>51.0</td><td>27.4</td><td>32.2</td><td>37.4</td><td>99.4</td><td>100.0</td><td>99.8</td></tr><tr><td></td><td>40</td><td>74.8</td><td>75.0</td><td>82.6</td><td>9.8</td><td>25.0</td><td>10.6</td><td>74.8</td><td>74.8</td><td>82.6</td><td>90.2</td><td>98.0</td><td>93.0</td></tr></tbody></table>

The type I error and power from clinical trial simulations are given as the observed proportion of significant tests performed on early tumor shrinkage at week 8 (ETS8) from 500 simulated clinical trials. Under $H_0$, the 95% prediction interval from binomial distribution is 3.3%–7.3% for 500 simulations. For the one-arm design, the number of individuals by data set is $N = 60$, and all individuals are allocated to the combination arm. The ETS8 distribution (of observations (OBS), predictions (IPRED), or true values (TRUE)) is compared with a theoretical median of ETS8 under cetuximab alone, performing a one-sample Wilcoxon test (two-sided, $\alpha = 5\%$). Different possible theoretical median of ETS8 values (true or not) are investigated. Two possible drug interactions are considered: additivity and synergy (see Models section).

<figure><img src="image_19.png" /><figcaption>**Figure 2** Boxplots of early tumor shrinkage at week 8 (ETS8) from observations (OBS; solid boxes, left subgroup for each treatment), predictions (IPRED; dotted boxes, middle subgroup for each treatment), and true values (TRUE; dashed boxes, right subgroup for each treatment) in the two-arm design. Whiskers indicate 5th to 95th percentiles. C500 is cetuximab alone vs. the combination treatment (M + C500).</figcaption></figure>

additivity and synergy scenarios and using either observed or predicted ETS8, the type I error is well controlled as it lies in the 95% prediction interval (3.3%–7.3%). Compared with the previous design, there is a decrease in power in the additivity scenario (using OBS: 66.6% vs. 89.0% in one-arm design). This result is explained by having only 30 patients per arm (vs. 60 in the one-arm design). Nevertheless, this design remains powerful in the case of synergy (using OBS: 97.4% vs. 100% in the one-arm design). Moreover, in each scenario, using the model-based ETS8 value leads to a slight increase in power compared with using the raw data. For example, the power in the additivity scenario was 69.2% using IPRED.

### Four-arms: C500/M + C500/M + C400/M + C200

For the four-arm design, the distributions of individual ETS8 values (OBS, IPRED, and TRUE) in each of the arms are presented in **Figure 3** for additivity and synergy. The IPRED distribution seems to be closer to TRUE than OBS, especially for extreme percentiles, except in the C500 arm. As expected, the distributions of C500 and M + C200 are close to each other for additivity. Moreover, in synergy, the ETS8 distribution is close between the two combination arms with the higher doses of cetuximab.

Type I error and power are presented in **Table 4**. Performing a global Kruskal–Wallis to detect a difference between the ETS8 distributions of the four arms, more power was

<table><caption>Table 3 Type I error and power of tests in two-arm design</caption><thead><tr><th rowspan="2"></th><th colspan="2">Additivity</th><th colspan="2">Synergy</th></tr><tr><th>Type I error (%)</th><th>Power (%)</th><th>Type I error (%)</th><th>Power (%)</th></tr></thead><tbody><tr><td>OBS</td><td>6.2</td><td>66.6</td><td>6.2</td><td>97.4</td></tr><tr><td>IPRED</td><td>6.2</td><td>69.2</td><td>7.0</td><td>98.6</td></tr><tr><td>TRUE</td><td>6.2</td><td>89.0</td><td>6.2</td><td>99.8</td></tr></tbody></table>

The type I error and power from clinical trial simulations are given as the observed proportion of significant tests performed on early tumor shrinkage at week 8 (ETS8) from 500 simulated clinical trials. Under $H_0$, the 95% prediction interval from binomial distribution is 3.3%–7.3% for 500 simulations. For the two-arms design, 30 patients are allocated to the combination arm, and 30 are allocated to cetuximab alone. The ETS8 distributions (of observations (OBS), predictions (IPRED), or true values (TRUE)) are compared between the two arms, performing a two-sample Wilcoxon test (two-sided, $\alpha = 5\%$). Two possible drug interactions are considered: additivity and synergy (see Models section).
<figure><img src="image_20.png" /><figcaption>**Figure 3** Boxplots of early tumor shrinkage at week 8 (ETS8) from observations (OBS; solid boxes, left subgroup for each treatment), predictions (IPRED; dotted boxes, middle subgroup for each treatment), and true values (TRUE; dashed boxes, right subgroup for each treatment) in the four-arm design. Whiskers indicate 5th to 95th percentiles. C500 is cetuximab alone, and M + C200, M + C400, and M + C500 are combination treatments.</figcaption></figure>

obtained in the synergy scenario than the additivity scenario. For example, using observed ETS8, the power was 69.4% in synergy vs. 30.2% in additivity. If the global test was significant, a Dunnett test was performed, comparing the ETS8 distribution in each of the three combination arms to those of the reference arm C500. Under the null hypothesis, if we compare the distribution of ETS8 for C500 and M + C500, the Dunnett test was conservative compared with a single Wilcoxon test, the type I error using observations was 2% with this procedure in the four-arm design vs. 6.2% for the Wilcoxon test in the two-arm design. Therefore, in the additivity scenario, the power to detect the superiority of M + C500 vs. C500 decreased to 19.8% using observations and 25.6%

using predictions (vs. 66.6% and 69.2% in two-arm design, respectively), which is partly explained by the lesser number of patients per arm (15 vs. 30). Nevertheless, in synergy and using model-based ETS8, we achieved 65.4% power to detect the superiority for M + C500 and 52% for M + C400. Thus, under certain conditions, this design is interesting for exploring other doses in the combination of two drugs.

## DISCUSSION

This work compares designs for early clinical combination trials in oncology in the context of a fixed dose of one agent and different doses of another with the same number of

<table><caption>**Table 4** Type I error and power of tests in four-arm design</caption><thead><tr><th></th><th colspan="2">Additivity</th><th colspan="2">Synergy</th></tr></thead><tbody><tr><td colspan="5">Power of Kruskal-Wallis test (%)</td></tr><tr><td>OBS</td><td colspan="2">30.2</td><td colspan="2">69.4</td></tr><tr><td>IPRED</td><td colspan="2">52.0</td><td colspan="2">76.6</td></tr><tr><td>TRUE</td><td colspan="2">44.0</td><td colspan="2">87.8</td></tr><tr><td colspan="5">Type I error of Dunnett test (%) C500 vs. M + C500</td></tr><tr><td>OBS</td><td colspan="2">1.6</td><td colspan="2">1.6</td></tr><tr><td>IPRED</td><td colspan="2">2.0</td><td colspan="2">2.0</td></tr><tr><td>TRUE</td><td colspan="2">2.4</td><td colspan="2">2.4</td></tr><tr><td></td><td>C500 vs. M + C200</td><td>C500 vs. M + C400</td><td>C500 vs. M + C500</td><td>C500 vs. M + C200</td><td>C500 vs. M + C400</td><td>C500 vs. M + C500</td></tr><tr><td colspan="7">Power of Dunnett test (%)</td></tr><tr><td>OBS</td><td>3.0</td><td>13.0</td><td>19.8</td><td>26.4</td><td>51.8</td><td>60.2</td></tr><tr><td>IPRED</td><td>5.4</td><td>13.2</td><td>25.6</td><td>18.6</td><td>52.0</td><td>65.4</td></tr><tr><td>TRUE</td><td>4.0</td><td>19.8</td><td>29.2</td><td>34.0</td><td>70.4</td><td>79.0</td></tr></tbody></table>

The type I error and power from clinical trial simulations are given as the observed proportion of significant tests performed on early tumor shrinkage at week 8 (ETS8) from 500 simulated clinical trials. For the four-arms design, 15 patients are allocated to each arm. The ETS8 distributions (of observations (OBS), predictions (IPRED), or true values (TRUE)) are compared between the arms, performing first a global Kruskal-Wallis test (two-sided, α = 5%), then, if the global test is significant, a Dunnett test where distribution of each combination arm (M + C200, M + C400, and M + C500) is compared with those of cetuximab alone (C500). Two possible drug interactions are considered: additivity and synergy (see Models section).
subjects and based on ETS. The one-arm design leads to inflation of type I error in the case of a misestimated reference value, as reported previously.⁴² Therefore, before the analysis of data from a one-arm design, we should ensure that the conditions of the study are similar to the reference study. With the same number of subjects and an accurate reference value, the power of a one-arm design is higher than that of a two-arm design. However, to avoid the risk of an inaccurate reference value, a two-arm randomized design including an arm with the reference therapy is preferable. Such a design conducts to control the type I error. Another possibility is to choose a four-arm design with different possible dosing regimens. Although with the same total number of patients a two-arm design leads to a greater power of tests, a four-arm design offers more accurate dose selection.

Our work aimed to show simulations in an ideal setting in order to facilitate the interpretation of the results. In these simulations, we did not account for possible incomplete adherence, missing data, protocol deviations (e.g., timing of TS measurements), or dropout. It would be interesting to incorporate these design departures in future simulations and analyses to mimic a real clinical study.⁴³ Nevertheless, most of these departures are not preponderant in the first 8 weeks of a trial and therefore should not change drastically our simulation results. In this study, exposure to cetuximab was simulated from a previous pharmacokinetics model²⁸ that included a BSA effect on clearance but other covariates, such as normal fat mass, would be more physiological.⁴⁴ Moreover, the allometric coefficient of 0.75 (used for the BSA effect in the model) would be more appropriate with the normal fat mass. In our study, the clearance of M was not influenced by any covariate, but the effect of body weight or the normal fat mass could have been considered. In the model that we used, we assumed linear effects in the exposure ranges of cetuximab and M based on previous models in which the maximal effect (Emax) could not be estimated. However, it would be relevant to explore more mechanistic models and assume nonlinear effects, as is the case for the interaction between cetuximab and M in the model. The TGI model involves a global resistance parameter KR. It would be relevant to associate different resistance parameters with each drug and its exposure.

Using NLMEM to analyze longitudinal data is a promising tool, especially when the outcome is based on TS. Tumor measurements are notably influenced by scan variability and the individual reviewing the image, which may result in noise in the data⁴⁵ and a loss of power. The NLMEM accounts for the information from the longitudinal data, which were five TS measurements in our study, and allows estimation of the residual error, including the noise in the data. Therefore, performing tests based on model-based predictions can lead to better power than observations. For example, in our two-arm design with the additivity scenario, the power using model-based ETS8 was 69.2% vs. 66.6% using observations.

In this study, ETS was chosen as outcome. We could also dichotomize ETS as responder or nonresponder outcome. Different thresholds can be chosen to separate the patients depending on their ETS.⁴⁶ We can dichotomize the

subjects as responders if ETS8 ≥ 30% or nonresponder if ETS8 < 30%, which corresponds roughly to a complete or partial response vs. stable or progressive disease according to the RECIST criteria. With this binary outcome, in our two-arm design and additivity scenario, the power of Fisher's exact test was 46.4% instead of 69.2% with the Wilcoxon test of continuous model-based ETS8 values.

Guidelines highlight the importance of design choice for an early-phase combination trial in oncology and that assumptions and knowledge about pharmacokinetic/pharmacodynamic interactions should be taken into account to select the most effective trial design.⁴⁷ There is no standard design for combination studies, and designs that include randomization helps to define the effect of the combination vs. one or both of the individual components.¹

We plan to extend this work to perform model-based adaptive two-stage designs⁴⁸,⁴⁹ using the Fisher information matrix to optimize the second stage of the study⁵⁰ in which different dose arms could also be added or dropped at the end of the first stage.

**Supporting Information.** Supplementary information accompanies this paper on the CPT: Pharmacometrics & Systems Pharmacology website (www.psp-journal.com).

**Acknowledgments.** The authors would like to thank Thu Thuy Nguyen and Vishnu Dutt Sharma who contributed to this project. They also thank Hervé Le Nagard and Lionel de la Tribouille for the use of the CATIBioMed calculus facility.

**Funding.** This work was supported by Merck KGaA, Darmstadt, Germany.

**Conflict of Interest.** Pascal Girard and Kosalaram Goteti are employed by Merck Institute for Pharmacometrics and EMD Serono R&D Institute, Inc., respectively, an affiliate of Merck KGaA, Darmstadt, Germany. All other authors declared no competing interests for this work.

**Author Contributions.** J.S., P.G., K.G., and F.M. wrote the manuscript. J.S., P.G., K.G., and F.M. designed the research. J.S., P.G., K.G., and F.M. performed the research. J.S., P.G., K.G., and F.M. analyzed the data.

**Disclaimer.** As Editor-in-Chief of CPT: Pharmacometrics and Systems Pharmacology, France Mentré was not involved in the review or decision process for this article.

1. Lo Russo, P.M., Boerner, S.A. & Seymour, L. An overview of the optimal planning, design, and conduct of phase I studies of new therapeutics. *Clin. Cancer Res.* **16**, 1710–1718 (2010).

2. Ledford, H. Cocktails for cancer with a measure of immunotherapy. *Nature* **532**, 162–164 (2016).

3. Sandler, A. et al. Paclitaxel-carboplatin alone or with bevacizumab for non-small-cell lung cancer. *N. Engl. J. Med.* **355**, 2542–2550 (2006).

4. Huang, J. et al. Emerging trends in us oncological approvals: a 13-year review (1999–2011). *Drug Inf. J.* **46**, 344–357 (2012).

5. Cunningham, D. et al. Cetuximab monotherapy and cetuximab plus irinotecan in irinotecan-refractory metastatic colorectal cancer. *N. Engl. J. Med.* **351**, 337–345 (2004).

6. Sullivan, R.J. et al. Atezolizumab plus cobimetinib and vemurafenib in BRAF-mutated melanoma patients. *Nat. Med.* **25**, 929–935 (2019).

7. Seymour, L.K. et al. Design and conduct of early clinical studies of two or more targeted anticancer therapies: recommendations from the task force on methodology
for the development of innovative cancer therapies. *Eur. J. Cancer Oxf. Engl. 1990* **49**, 1808–1814 (2013).

8. Smoragiewicz, M. et al. Design and conduct of early clinical studies of immunotherapy: recommendations from the task force on Methodology for the Development of Innovative Cancer Therapies 2019 (MDICT). *Clin. Cancer Res.* **26**, 2461–2465 (2020).

9. Paller, C.J. et al. Factors affecting combination trial success (FACTS): investigator survey results on early-phase combination trials. *Front. Med.* **6** (2019). https://doi.org/10.3389/fmed.2019.00122.

10. Eisenhauer, E.A. et al. New response evaluation criteria in solid tumours: Revised RECIST guideline (version 1.1). *Eur. J. Cancer* **45**, 228–247 (2009).

11. Litière, S. et al. RECIST 1.1 for response evaluation apply not only to chemothera- py-treated patients but also to targeted cancer agents: a pooled database analysis. *J. Clin. Oncol.* **37**, 1102–1110 (2019).

12. Sharma, M.R., Maitland, M.L. & Ratain, M.J. RECIST: no longer the sharpest tool in the oncology clinical trials toolbox–point. *Cancer Res.* **72**, 5145–5149 (2012), discussion 5150.

13. Vera-Yunca, D. et al. Machine learning analysis of individual tumor lesions in four metastatic colorectal cancer clinical studies: linking tumor heterogeneity to overall survival. *AAPS J.* **22**, 58 (2020).

14. Wang, Y. et al. Elucidation of relationship between tumor size and survival in non- small-cell lung cancer patients can aid early decision making in clinical drug devel- opment. *Clin. Pharmacol. Ther.* **86**, 167–174 (2009).

15. Claret, L. et al. Model-based prediction of phase III overall survival in colorectal can- cer on the basis of phase II tumor dynamics. *J. Clin. Oncol.* **27**, 4103–4108 (2009).

16. Colloca, G.A., Venturino, A. & Guarneri, D. Early tumor shrinkage after first-line medical treatment of metastatic colorectal cancer: a meta-analysis. *Int. J. Clin. Oncol.* **24**, 231–240 (2019).

17. Kimko, H. & Pinheiro, J. Model-based clinical drug development in the past, present and future: a commentary. *Br. J. Clin. Pharmacol.* **79**, 108–116 (2015).

18. Petrylak, D.P. et al. Evaluation of prostate-specific antigen declines for surrogacy in patients treated on SWOG 99–16. *J. Natl. Cancer Inst.* **98**, 516–521 (2006).

19. Gompertz, B. On the nature of the function expressive of the law of human mortal- ity, and on a new mode of determining the value of life contingencies. *Philos. Trans. R. Soc. Lond.* **115**, 513–583 (1825).

20. Simeoni, M. et al. Predictive pharmacokinetic-pharmacodynamic modeling of tumor growth kinetics in xenograft models after administration of anticancer agents. *Cancer Res.* **64**, 1094–1101 (2004).

21. Bonate, P.L. & Suttle, A.B. Modeling tumor growth kinetics after treatment with pazopanib or placebo in patients with renal cell carcinoma. *Cancer Chemother. Pharmacol.* **72**, 231–240 (2013).

22. Tardivon, C. et al. Association between tumor size kinetics and survival in patients with urothelial carcinoma treated with atezolizumab: implication for patient fol- low-up. *Clin. Pharmacol. Ther.* **106**, 810–820 (2019).

23. Jonker, D.J. et al. Cetuximab for the treatment of colorectal cancer. *N. Engl. J. Med.* **357**, 2040–2048 (2007).

24. Garrett, C.R. et al. Phase I dose-escalation study to determine the safety, pharma- cokinetics and pharmacodynamics of brivanib alaninate in combination with full- dose cetuximab in patients with advanced gastrointestinal malignancies who have failed prior therapy. *Br. J. Cancer* **105**, 44–52 (2011).

25. Van Cutsem, E. et al. Binimetinib, encorafenib, and cetuximab triplet therapy for patients with BRAF V600E-mutant metastatic colorectal cancer: safety lead-in results from the phase III beacon colorectal cancer study. *J. Clin. Oncol.* **37**, 1460–1469 (2019).

26. Le Louedec, F. et al. Cetuximab pharmacokinetic/pharmacodynamics relationships in advanced head and neck carcinoma patients. *Br. J. Clin. Pharmacol.* **85**, 1357– 1366 (2019).

27. Azzopardi, N. et al. Cetuximab pharmacokinetics influences progression-free survival of metastatic colorectal cancer patients. *Clin. Cancer Res.* **17**, 6329–6337 (2011).

28. Girard, P. et al. Drug-disease model of tumor size and overall survival in metastatic colorectal cancer patients treated with cetuximab administered weekly or every second week. European Cancer Congress, Amsterdam, Netherlands. (2013).

29. Tan, A.R. et al. Pharmacokinetics of cetuximab after administration of escalating single dosing and weekly fixed dosing in patients with solid tumors. *Clin. Cancer Res.* **12**, 6517–6522 (2006).

30. Tabernero, J. et al. Cetuximab administered once every second week to patients with metastatic colorectal cancer: a two-part pharmacokinetic/pharmacodynamic phase I dose-escalation study. *Ann. Oncol.* **21**, 1537–1545 (2010).

31. Bokemeyer, C. et al. Fluorouracil, leucovorin, and oxaliplatin with and without cetuximab in the first-line treatment of metastatic colorectal cancer. *J. Clin. Oncol.* **27**, 663–671 (2009).

32. Van Cutsem, E. et al. Cetuximab and chemotherapy as initial treatment for meta- static colorectal cancer. *N. Engl. J. Med.* **360**, 1408–1417 (2009).

33. Cheng, A.-L. et al. Efficacy, tolerability, and biomarker analyses of once-every-2- weeks cetuximab plus first-line FOLFOX or FOLFIRI in patients with KRAS or all RAS wild-type metastatic colorectal cancer: the phase 2 APEC study. *Clin. Colorectal Cancer* **16**, e73–e88 (2017).

34. Fernandez-Plana, J. et al. Biweekly cetuximab in combination with FOLFOX-4 in the first-line treatment of wild-type KRAS metastatic colorectal cancer: final results of a phase II, open-label, clinical trial (OPTIMIX-ACROSS Study). *BMC Cancer* **14**, 865 (2014).

35. Wicha, S.G., Chen, C., Clewe, O. & Simonsson, U.S.H. A general pharmacodynamic interaction model identifies perpetrators and victims in drug interactions. *Nat. Commun.* **8**, 2129 (2017).

36. Genz, A. Numerical computation of multivariate normal probabilities. *J. Comput. Graph. Stat.* **1**, 141–149 (1992).

37. Soetaert, K., Petzoldt, T. & Setzer, R.W. Solving differential equations in *R*: package desolve. *J. Stat. Softw.* **33**, 9 (2010).

38. Kuhn, E. & Lavielle, M. Maximum likelihood estimation in nonlinear mixed effects models. *Comput. Stat. Data Anal.* **49**, 1020–1038 (2005).

39. Ahn, J.E., Karlsson, M.O., Dunne, A. & Ludden, T.M. Likelihood based approaches to handling data below the quantification limit using NONMEM VI. *J. Pharmacokinet. Pharmacodyn.* **35**, 401–421 (2008).

40. Benjamini, Y. & Hochberg, Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. *J. R. Stat. Soc. Ser. B Methodol.* **57**, 289–300 (1995).

41. Clopper, C.J. & Pearson, E.S. The use of confidence or fiducial limits illustrated in the case of the binomial. *Biometrika* **26**, 404–413 (1934).

42. Viele, K. et al. Use of historical control data for assessing treatment effects in clinical trials. *Pharm. Stat.* **13**, 41–54 (2014).

43. Sheiner, L.B. Is intent-to-treat analysis always (ever) enough? *Br. J. Clin. Pharmacol.* **54**, 203–211 (2002).

44. Holford, N.H.G. & Anderson, B.J. Allometric size: The scientific theory and exten- sion to normal fat mass. *Eur. J. Pharm. Sci.* **109S**, S59–S64 (2017).

45. Sullivan, D.C., Schwartz, L.H. & Zhao, B. The imaging viewpoint: how imaging af- fects determination of progression-free survival. *Clin. Cancer Res.* **19**, 2621–2628 (2013).

46. Mandrekar, S.J. et al. Evaluation of alternate categorical tumor metrics and cut points for response categorization using the RECIST 1.1 data warehouse. *J. Clin. Oncol.* **32**, 841–850 (2014).

47. Paller, C.J. et al. Design of phase I combination trials: recommendations of the clinical trial design task force of the nci investigational drug steering committee. *Clin. Cancer Res.* **20**, 4210–4217 (2014).

48. Foo, L.K. & Duffull, S. Adaptive optimal design for bridging studies with an ap- plication to population pharmacokinetic studies. *Pharm. Res.* **29**, 1530–1543 (2012).

49. Pierrillas, P.B. et al. Model-based adaptive optimal design (MBAOD) improves combination dose finding designs: an example in oncology. *AAPS J.* **20**, 39 (2018).

50. Mentré, F., Mallet, A. & Baccar, D. Optimal design in random-effects regression models. *Biometrika* **84**, 429–442 (1997).

© 2020 The Authors. CPT: Pharmacometrics & Systems Pharmacology published by Wiley Periodicals LLC on behalf of the American Society for Clinical Pharmacology and Therapeutics. This is an open access article under the terms of the Creative Commons Attribution-NonCommercial-NoDerivs License, which permits use and distribution in any medium, provided the original work is properly cited, the use is non-commercial and no modifications or adaptations are made.
<figure><img src="image_21.png" /></figure>

**Figure S1.** Median tumor size for the additivity (left) and synergy (right) scenario with different therapies. Cetuximab alone at 500 mg/m² every 2 weeks (C500) in blue and solid line, or combination treatment with different doses of cetuximab: 200, 400, or 500 mg/m² every 2 weeks (respectively M+C200, M+C400 and M+C500, respectively) in different shades of grey and purple for additivity and synergy, respectively.
<figure><img src="image_22.png" /><figcaption>**Figure S2.** Simulated tumor size data for 30 individuals under cetuximab alone (C500, in blue) or combination treatment (M+C500) for the additivity (in black) or synergy scenario (in purple).

(A) The red horizontal line represents the lower limit of quantification of 10 mm.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(B) Without residual error and an initial tumor size of 100 mm.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>**Figure S3.** Simulated observations (OBS, black dots), model based predictions (IPRED, green line and dots), and true values (TRUE, red line and dots) for tumor sizes in one individual.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>TRUE ETS8 (%)</figcaption></figure>

**Figure S4.** Individual model-based (IPRED) vs. true ETS8 (TRUE) values with different designs by row and under different interaction scenarios by column.
# Troisième partie : Modélisation de l'interaction bactériophages/bactéries pour le traitement de la pneumonie

Dans le contexte de l'inquiétante augmentation de l'émergence et de la propagation des bactéries multi-résistantes, il est nécessaire de mener des actions contre la résistance aux antibiotiques avec une forte incitation à mettre au point de nouvelles approches thérapeutiques (O'Neill 2014). Parmi les solutions les plus prometteuses et radicalement différentes des molécules classiques, la phagothérapie suscite de plus en plus l'intérêt des communautés scientifiques et médicales (Caflisch et al. 2019). Cette approche utilise des bactériophages *i.e.* des virus qui n'infectent que les bactéries. Cette phagothérapie fût rapidement adoptée dans un premier temps, mais le manque de connaissance scientifique sur les bactériophages, ainsi que le rapide développement de la production industrielle des antibiotiques a sonné le glas de ce traitement, à l'exception de certains pays de l'Europe de l'Est où la phagothérapie s'est fortement développée. Les résultats prometteurs des études animales suggèrent que la phagothérapie pourrait être envisagée comme une alternative sérieuse aux antibiotiques. Néanmoins, il reste encore beaucoup à apprendre avant que l'on puisse proposer un traitement à base de phages. En particulier, il est nécessaire de comprendre et de quantifier les interactions phages/bactéries afin de définir les modes d'administration de phages les plus appropriés.

## 3.1. Résumé

Les travaux présentés s'inscrivent dans un projet global impliquant l'Institut Pasteur (équipe « Bactériophage, bactérie, hôte ») et IAME (équipes EVRest et BIPID) qui a pour but de mieux appréhender la phagothérapie. La série d'expériences présentée a été menée sur les bactériophages 536P1, dont les propriétés thérapeutiques ont été précédemment validées dans un modèle d'infection pulmonaire aiguë à E.coli 536 (Dufour et al. 2015). Nous avons construit un modèle mathématique à partir de ces données afin de caractériser les aspects essentiels concernant les interactions entre phages et bactéries. Au cours de cette nouvelle étude, des expériences in vitro, puis chez la souris, ont été conduites pour le couple 536/536_P1, et l'analyse successive de ces données nous a permis de construire un modèle mathématique décrivant l'évolution de l'infection bactérienne et l'efficacité de la phagothérapie.
Les données *in vitro* consistent en trois types d'expériences. Pour l'étude de *one-step growth*, la multiplication des phages est observée pour un unique cycle lytique. L'étude de co-évolution a permis de suivre les concentrations bactérienne et virale de manière longitudinale. Enfin, les cinétiques de lyses permettent de suivre la concentration bactérienne selon différentes doses connues de phages. Ces expériences ont notamment permis d'estimer les paramètres du cycle lytique du bactériophage 536_P1 en présence de la bactérie 536.

Les études chez la souris saine, où la concentration pulmonaire de phages a été mesurée lors du sacrifice (à différents temps selon les souris), ont permis de caractériser la pharmacocinétique des bactériophages, après administration intraveineuse (IV, N = 37) ou intratrachéale (IT, N = 16).

Enfin, des souris ont été infectées puis traitées ou non par bactériophages à différentes doses et voies d'administration (IV et IT). Grâce à une technique originale basée sur l'imagerie des bactéries bioluminescentes, les charges bactériennes pulmonaires ont pu être quantifiées à différents temps chez les mêmes animaux. La concentration pulmonaire de phages a été relevée uniquement au temps de sacrifice (si la souris survivait jusqu'au temps prévu).

Certains paramètres du modèle final décrivant les données de souris infectées (N = 183) ont été fixés à partir des estimations obtenues à partir de modèles construits en amont, notamment celui décrivant les données *in vitro* et celui décrivant la pharmacocinétique des phages chez la souris saine. L'ensemble des paramètres ont été estimés par MNLEM, avec l'algorithme SAEM implémenté dans le logiciel Monolix (version 2018R2), à l'exception des paramètres de la pharmacocinétique des phages estimés par naïve pooling of data. Le modèle conjoint final a été évalué par VPC corrigées (Bergstrand et al., 2011) de la charge bactérienne et par VPC de survie. A partir du modèle final, il a été étudié l'impact de la dose et de la voie d'administration des phages, de l'inoculum bactérien et des paramètres du cycle lytique sur la dynamique médiane des phages et des bactéries.

Les principaux paramètres du cycle lytique ont été déterminés à partir des expériences *in vitro* et il a notamment été estimé que 570 phages sont relargués lors de la lyse d'une bactérie 536 (*burst size*). A partir des données chez la souris non infectée, c'est-à-dire sans amplification possible des phages, l'administration de phages par voie IV conduisait à une exposition 5000 fois moindre qu'après administration locale (IT). L'analyse des données chez la souris infectée, afin de caractériser les paramètres clés de l'interaction phage-bactérie, a nécessité l'utilisation d'un système d'équations différentielles ordinaires combinant un modèle de cinétique bactérienne et un modèle de dynamique virale. Le modèle final consiste en un modèle proie-
prédateur (bactérie-bactériophage, respectivement) qui intègre le rôle du système immunitaire sur la décroissance bactérienne. Le modèle a montré que pour un faible inoculum bactérien, par exemple avec une concentration bactérienne initiale de 6 log CFUeq/g, le système immunitaire seul permet de contrôler l'infection. Pour un inoculum modéré (aux alentours de 7 log CFUeq/g), les phages et le système immunitaire agissent de concert pour parvenir à réduire la charge bactérienne. Enfin, lorsque l'inoculum bactérien est très important (par exemple supérieur à 8 log CFUeq/g), la phagothérapie ne permet pas de contrôler la charge bactérienne. Le modèle montre aussi un lien significatif (p < 10⁻³ par test de Wald) entre la charge bactérienne des souris et leur survie, par le biais d'un modèle conjoint. Le risque instantané de décès était multiplié par 18 pour une charge bactérienne augmentée de 2 log CFUeq/g, ce qui correspond à un *hazard ratio* de 4,2. Ce modèle était correctement évalué par les différents graphiques diagnostiques et a permis de décrire les différents profils individuels malgré une forte variabilité entre ces derniers.

Le modèle final a permis de quantifier l'efficacité de la phagothérapie en fonction de la dose et de la voie d'administration : pour un inoculum bactérien létal (7 log CFUeq/g) pour 100 % des souris non traitées, l'administration de phages, permettait un contrôle de l'infection bactérienne chez la plupart des individus, pour tous les groupes de doses et de voie d'administration. En médiane, la charge bactérienne commençait à décroître moins d'une demi-heure après administration IT à haute dose (IT Hi), environ une heure et demie après IT à faible dose (IT Lo, dose 1000 fois plus importante que la dose faible), et moins de 4 h après administration IV à faible dose (IV Lo). Ces différences, notamment entre les voies d'administrations, pouvaient devenir plus importantes pour un phage présentant un *burst size* 10 fois moins important. Les charges bactériennes commençaient en médiane à décroître à 0,3, 2,7 et 11,3 h après administration, chez les IT Hi, IT Lo et IV Lo, respectivement. Ces travaux montrent que la voie locale est à privilégier par rapport à la voie IV, en particulier lorsque les phages présentent une mauvaise capacité d'amplification. Le modèle construit permet donc bien d'évaluer l'efficacité de la phagothérapie en fonction de la dose et du mode d'administration. Ce travail permettra d'aider au choix de protocole lors de futures études sur les bactériophages. Un article présentant ce travail est soumis.

## 3.2. Article 4 (soumis)
*Combination of in vivo phage therapy data with in silico model highlights key parameters for treatment efficacy, soumis pour publication*

Raphaëlle Delattre*, Jérémy Seurat*, Feyrouz Haddad, Thu Thuy Nguyen, Baptiste Gaborieau, Rokhaya Kane, Nicolas Dufour, Jean-Damien Ricard, Jérémie Guedj, Laurent Debarbieux
**Title**

Combination of *in vivo* phage therapy data with *in silico* model highlights key parameters for treatment efficacy

**Authors**

Raphaëlle Delattre<sup>1,2,3*</sup>, Jérémy Seurat<sup>2*</sup>, Feyrouz Haddad<sup>1,2</sup>, Thu-Thuy Nguyen<sup>2</sup>, Baptiste Gaborieau<sup>1,2,3</sup>, Rokhaya Kane<sup>1</sup>, Nicolas Dufour<sup>4</sup>, Jean-Damien Ricard<sup>2</sup>, Jérémie Guedj<sup>2**</sup> et Laurent Debarbieux<sup>1**</sup>

**Affiliations**

1 Institut Pasteur, Bacteriophage, Bacterium, Host Laboratory, Paris F-75015 France

2 Université de Paris, INSERM, IAME, Paris F-75006, France

3 Université de Paris, Cellule Pasteur, Paris F-75013, France

4 Centre Hospitalier René Dubos, Cergy Pontoise F-95503, France

Lead contact: Laurent Debarbieux (laurent.debarbieux@pasteur.fr)

<footnote>* : co-first authors</footnote>

<footnote>** : co-corresponding authors</footnote>
**Summary**

The clinical (re)development of phage therapy to treat antibiotic resistant infections requires grasping specific biological properties of bacteriophages (phages) as antibacterial. However, identification of optimal dosing regimens is hampered by the poor understanding of phage-bacteria interactions *in vivo*. Here we developed a general strategy coupling *in vitro* and *in vivo* experiments with a mathematical model to characterize the interplay between phage and bacterial dynamics during pneumonia induced by a pathogenic strain of *Escherichia coli*. The model estimates some key parameters for phage therapeutic efficacy, in particular the impact of dose and route of administration on phage dynamics and the synergism of phage and the innate immune response on the bacterial clearance rate. Simulations predict a low impact of the intrinsic phage characteristics in agreement with the current semi-empirical choices of phages for compassionate treatments. Model-based approaches will foster the deployment of future phage therapy clinical trials.

**Keywords**

Bacteriophage, antibiotic resistance, bacterial infection, pneumonia, mathematical model, innate immunity, simulation
## Introduction

The pronounced slowdown in the discovery of new antibiotics to treat bacterial infections caused by multidrug resistant (MDR) pathogens has revived the interest for bacteriophages (phages), viruses infecting bacteria (Gordillo Altamirano and Barr, 2019; Kortright et al., 2019). Pre-clinical experiments and several clinical cases reports have demonstrated the potential of phages to successfully treat infections caused by MDR bacteria (Aslam et al., 2020; Corbellino et al., 2020; Jennes et al., 2017; Melo et al., 2020; Schooley et al., 2017). These success are nonetheless balanced by clinical studies that did not confirm these positive results (Jault et al., 2019; Leitner et al., 2020; Sarker et al., 2016), showing that translation from experimental results and promising individual treatments to clinical trials remains challenging.

A major milestone for a widespread use in human populations is the determination of the optimal dose, route of administration and treatment duration. This is a particularly complex endeavor for phages, as standard assessments of clinical pharmacology used to determine the processes of Administration, Distribution, Metabolism and Excretion (ADME) of drugs are not adapted to phages. Indeed, phages massively replicate on their target bacteria with kinetic parameters governing this viral amplification (receptor recognition, genome amplification, particles assembly and bacterial lysis) being specific to each phage-bacteria pair (Dion et al., 2020). Furthermore, phages are naturally present in human-associated microbiotas and their rate and route of elimination do not follow usual metabolic pathways through kidney or liver, as for xenobiotic products (Dabrowska and Abedon, 2019).

Here we coupled *in vitro* and *in vivo* experiments with mathematical modeling to approach the key parameters of phage-bacteria-host interactions in order to guide treatment strategies. This approach has been successful in drug development for various chronic (HIV, HCV) and acute (Ebola, Zika, Influenza) viral infections (Friberg and Guedj, 2020; Perelson and Guedj, 2015). Unlike previous models that focused on specific aspects of phage-bacteria, phage-host immune system or bacteria-host immune system interactions (Cairns et al., 2009; Roach et al., 2017; Smith et al., 2011; Weld et al., 2004), we aimed to design a synthetic, semi-mechanistic model that could be used to understand the tripartite phage-bacteria-host interactions *in vivo*. We purposely designed a series of specific experiments (Figure 1) to estimate the central parameters from a previously validated experimental mice model of pulmonary infection using the bioluminescent *Escherichia coli* pathogenic strain 536 and the virulent phage 536_P1 (Dufour et al., 2015). Dynamic and direct microbiological data were collected from groups of animals that were either (bacterial) uninfected and (phage) treated, infected and untreated, or infected and treated, using different inoculum doses and routes of administration. These data led us to
develop an original mathematical model recapitulating the tripartite interactions over time. This
model identified key parameters affecting bacterial kinetics in vivo. It was then exploited to test
scenarios including variations of phage intrinsic properties. The use of such model will
rationalize the phage choice, the dose and route of administration in order to optimize the
efficacy of phage therapy.

<figure><img src="image_26.png" /><figcaption>**Figure 1. Schematic representation of the experimental design of the study.**</figcaption></figure>

*In vivo* experiments were conducted in mice that were either i) uninfected and treated (green box, with phages administered either intratracheally, IT, or intravenously, IV); ii) infected and untreated (magenta box); iii) infected and IT- or IV-treated at various doses (yellow box). The intensity of light emitted by *E. coli* strain 536-lux was recorded along all experiments and translated to bacterial lung density, while phage 536_P1 density was measured only once (upon sacrifice) in each animal at specific time points. A mathematical compartmental model was developed to fit to these data, with S, R, I₁ and I₂ representing the susceptible, refractory, infected non-productive, productively infected bacteria, respectively, and V and Vblood the phage concentration in lungs and blood, respectively.
# Results

## Biodistribution of phage in uninfected mice

We first studied the biodistribution of phage 536_P1 in uninfected mice following a single administration of 8 log₁₀ PFU by intratracheal (IT) or intravenous (IV) route (Figure 2). After IT administration, phage concentrations in the lungs varied over time only within 1-log with observed mean values between 7.0 and 8.0 log₁₀ PFU/g. By contrast, after IV administration, the observed peak concentration was reached four hours post-administration (5.1 [4.7, 5.5] log₁₀ PFU/g), which was over 2-log below the observed concentration after IT administration (7.2 ± [6.8, 7.8] log₁₀ PFU/g). Therefore, after IV administration only a fraction *F* of the phage, estimated to 0.02%, ultimately reached the lung compartment, leading to exposure approximately 5000 times lower than after IT administration and an estimation of the absorption rate into the lung compartment, *k*<sub>a</sub>, of 0.3 h⁻¹. (Eq. 1-2). The kinetics of phage elimination was also different with an estimated half-life of 12 and 3 hours after IT and IV administration, respectively (see *k*<sub>e</sub>(IT) and *k*<sub>e</sub>(IV) values in Table 1).

<figure><img src="image_27.png" /><figcaption>**Figure 2. Distribution of phage 536_P1 in lungs of treated uninfected and treated infected mice.**</figcaption></figure>

Lungs of uninfected (left, N=53) and *E. coli* infected (right, N=138) mice were homogenized at 2, 4, 8, 24 or 48h post-administration of a single dose of phage 536_P1 (8 log₁₀ PFU in uninfected mice and 5.6, 7.6, or 8.6 log₁₀ PFU in infected mice), either by intratracheal (blue) or intravenous (red) route. Lines represent the median phage concentrations at each time point. The dashed line represents the limit of quantification (LOQ=2 log₁₀ PFU/g).
<table><caption>**Table 1.** Parameter estimates of the final phage/bacteria dynamics model in mice</caption><thead><tr><th>Parameter</th><th>Name</th><th>Unit</th><th>Fixed/Estimated from</th><th>Fixed effect (rse, %)</th><th>Sd of random effect (rse, %)</th></tr></thead><tbody><tr><td>Bacterial growth rate</td><td>α</td><td>h<sup>-1</sup></td><td>Untreated infected</td><td>0.21 (-)</td><td>0 (-)</td></tr><tr><td>Initial number of bacteria</td><td>B<sub>0</sub></td><td>log<sub>10</sub> CFU g<sup>-1</sup></td><td>Untreated infected</td><td>7.1 (1)</td><td>0.14 (6)</td></tr><tr><td>Bacterial load plateau</td><td>B<sub>max</sub></td><td>log<sub>10</sub> CFU g<sup>-1</sup></td><td>Untreated infected</td><td>10.4 (-)</td><td>0 (-)</td></tr><tr><td>Time for immunity activation</td><td>t<sub>I</sub></td><td>h</td><td>Untreated infected</td><td>3.0 (-)</td><td>0 (-)</td></tr><tr><td>Maximal immunity effect</td><td>δ<sub>M</sub></td><td>CFU g<sup>-1</sup> h<sup>-1</sup></td><td>Untreated infected</td><td>6.7 10<sup>6</sup> (-)</td><td>0.05 (-)</td></tr><tr><td>Bacterial load for 50% of δ</td><td>B<sub>50</sub></td><td>CFU g<sup>-1</sup></td><td>Untreated infected</td><td>2.3 10<sup>7</sup> (-)</td><td>0 (-)</td></tr><tr><td>Infectivity rate</td><td>β</td><td>log<sub>10</sub> g CFU<sup>-1</sup> h<sup>-1</sup></td><td>Treated infected</td><td>6.9 (7)</td><td>0.11 (43)</td></tr><tr><td>Transition rate</td><td>e</td><td>h<sup>-1</sup></td><td>in vitro</td><td>1.2 (-)</td><td>0 (-)</td></tr><tr><td>Lysis rate</td><td>δ<sub>V</sub></td><td>h<sup>-1</sup></td><td>Treated infected</td><td>0.09 (25)</td><td>0 (-)</td></tr><tr><td>Burst size</td><td>b</td><td>PFU CFU<sup>-1</sup></td><td>in vitro</td><td>570 (-)</td><td>0 (-)</td></tr><tr><td>Relative bioavailability (IV/IT)</td><td>F</td><td>-</td><td>Treated uninfected</td><td>0.0002 (-)</td><td>0 (-)</td></tr><tr><td>Phage absorption rate</td><td>k<sub>α</sub></td><td>h<sup>-1</sup></td><td>Treated uninfected</td><td>3 (-)</td><td>0 (-)</td></tr><tr><td>Phage elimination rate after IT</td><td>k<sub>e(IT)</sub></td><td>h<sup>-1</sup></td><td>Treated uninfected</td><td>0.058 (-)</td><td>0 (-)</td></tr><tr><td>Phage elimination rate after IV</td><td>k<sub>e(IV)</sub></td><td>h<sup>-1</sup></td><td>Treated uninfected</td><td>0.23 (-)</td><td>0 (-)</td></tr><tr><td>Phage distribution weight</td><td>W</td><td>g</td><td>Treated uninfected</td><td>1 (-)</td><td>0 (-)</td></tr><tr><td>Initial refractory proportion</td><td>f</td><td>-</td><td>Treated infected</td><td>0.105 (45)</td><td>2.8 (19)</td></tr><tr><td>Weibull shape</td><td>κ</td><td>-</td><td>Treated infected</td><td>1.8 (11)</td><td>0 (-)</td></tr><tr><td>Weibull scale</td><td>λ</td><td>h<sup>-1</sup></td><td>Treated infected</td><td>2.4 10<sup>4</sup> (98)</td><td>0 (-)</td></tr><tr><td>Impact of bacterial kinetics on death</td><td>γ</td><td>g CFU<sup>-1</sup></td><td>Treated infected</td><td>1.44 (13)</td><td>0 (-)</td></tr><tr><td>Additive error on bacterial load</td><td>σ<sub>B</sub></td><td>log<sub>10</sub> CFU g<sup>-1</sup></td><td>Treated infected</td><td>-</td><td>0.47 (4)</td></tr><tr><td>Additive error on phage concentration</td><td>σ<sub>V</sub></td><td>log<sub>10</sub> PFU g<sup>-1</sup></td><td>Treated infected</td><td>-</td><td>2.6 (11)</td></tr></tbody></table>

Sd: standard deviation, rse: relative standard error

## Bacterial kinetics in infected-untreated mice

IT administration of the *E. coli* strain 536 led to a large variability in the lung bacterial load of untreated animals, as measured by emitted luminescence and CFUeq/g (see methods), with observed values ranging between 5.0 and 9.9 log₁₀ CFUeq/g at two hours post infection (Figure 3). This heterogeneity in the effective bacterial inoculum mimics the diversity of bacterial loads in human patients with pneumonia (Jiang et al., 2014), and allowed us to consider each mouse as an individual when developing the mathematical model (see below and Eq. 3). This heterogeneity also informed on the threshold of bacterial density above which the immune response of the host cannot control the infection. Indeed, when the bacterial load was greater
than 6.8 log₁₀ CFUeq/g at two hours post-infection, the bacterial density increased continuously, while, below this threshold, we observed a rapid reduction of the bacterial density caused by the immune response to the infection (Figure 3).

<figure><img src="image_28.png" /><figcaption>**Figure 3. Kinetics of bacterial densities in untreated and phage-treated animals.**</figcaption></figure>

Bioluminescence data from infected animals were stratified by treatment groups: untreated (black), treated with phage 536_P1 administered by IV at a low dose (5.6 log₁₀ PFU, yellow) or medium dose (7.6 log₁₀ PFU, red), or by IT at either low dose (cyan), or medium dose (blue) or high dose (8.6 log₁₀ PFU, dark blue). The dotted horizontal line at 6.8 log₁₀ CFUeq/g corresponds to the threshold value of bacterial load at two hours post-infection below which all animals controlled the infection. The dashed line represents the limit of quantification (LOQ=4.7 log₁₀ CFUeq/g). Circles correspond to live animals until the time of sacrifice while triangles correspond to dead animals before the scheduled time of sacrifice.

## Bacterial and phage kinetics in infected treated mice

Next, using the same infection procedure as above, mice (N=138) were treated two hours post-infection with various doses of phage 536_P1 (low dose: 5.6; medium dose: 6.6; high dose: 7.6 log₁₀ PFU), administered either by IT or IV route.

Phages rapidly replicated and their density in the lungs reached values 2 to 3-log₁₀ higher than in uninfected animals, showing that they encountered bacteria and efficiently replicated (Figure 2). At four hours post-phage administration (six hours post-bacterial infection), the median titer of phage in the lungs of IT-treated mice (N=59) was 9.6 [9.3, 9.9] log₁₀ PFU/g, as compared to 7.1 [7.1, 7.4] log₁₀ PFU/g in uninfected IT-treated mice (N=16). A similar difference was also observed in IV-treated animals (N=79), with phages concentration of 8.4 [7.0, 9.1] log₁₀ PFU/g as compared to 5.1 [4.7, 5.5] log₁₀ PFU/g in uninfected IV-treated animals (N=37). A dose
effect was noticed in IV-treated animals, with a median concentration of 7.6 [6.4, 8.5] log₁₀ PFU/g in the lungs of mice receiving a low dose (N=50), as compared to 10.6 [10.0,10.6] log₁₀ PFU/g in mice treated with a medium dose (N=29).

The initial effective median bacterial inoculum was similar in each of the three groups with levels of 6.9 [6.1, 8.0], 7.2 [6.3, 8.5] and 6.7 [6.3, 7.8] log₁₀ CFUeq/g in the lungs of untreated, IT- and IV-treated animals, respectively (Figure 3). Afterwards, two hours post-phage administration, the bacterial loads tended to be lower in treated (N=138) compared to untreated (N=45) animals, with levels of 7.5 [6.6, 8.8] log₁₀ CFUeq/g and 8.4 [7.6, 8.8] log₁₀ CFUeq/g, respectively (p=0.17, Wilcoxon-Mann-Whitney test). However, at this time point there was no effect of the route of administration, with values of 7.5 [6.5, 8.8] log₁₀ CFUeq/g and 7.2 [6.6, 8.5] log₁₀ CFUeq/g in IT and IV-treated groups, respectively (p=0.27 and p=0.14 to untreated, respectively, and p=0.75 when comparing IT and IV-treated groups, Wilcoxon-Mann-Whitney test).

The mortality rates at 48 hours were similar in all groups (respectively 42%, 48% and 52% in IT, IV and untreated animals) (Figure 4). Mortality rate was not dependent on the dose of administered phages (at 48 hours equal to 51%, 25% and 50% in animals receiving a low, medium or high dose of phages, respectively; p=0.92, 0.14 and 0.68 to untreated, logrank test), (Figure S4). However, it is important to note that in the subgroup of animals with a high effective bacterial inoculum (>6.8 log₁₀ CFUeq/g), treatment tended to increase the survival rate (Figure 4, p=0.07 in IT and p=0.21 in IV, as compared to untreated, respectively, logrank test), whereas all animals with a low effective inoculum (<6.8 log₁₀ CFUeq/g) survived at 48h, regardless of treatment (Figure 4).

<figure><img src="image_29.png" /><figcaption>**Figure 4.** Survival of infected and treated animals in function of the pulmonary bacterial load at two hours post-infection</figcaption></figure>
The survival curves were obtained from animals that were infected untreated (black, N=45), infected treated by IT route (blue, N=59) or infected treated by IV route (red, N=79). The left panel displays all infected animals (N=183). The middle panel includes only animals (N=78) with a pulmonary bacterial load at two hours post-infection greater than 6.8 log₁₀ CFUeq/g. The right panel includes only animals (N=99) with a pulmonary bacterial load at two hours post-infection lower than 6.8 log₁₀ CFUeq/g. Crosses represent sacrifice time. Of note, 6 animals sacrificed at two hours post infection were not represented in both middle and right panel.

Unraveling bacteria-phage interaction using a mathematical model

Next, we used these data in both infected and uninfected animals together with parameters
determined *in vitro* (see methods) to build a mathematical model of phage-bacteria-host
interaction (Eq. 4-10, Figure 1). Fitting the model to bacterial load data, we first determined
that a model accounting for the control of bacterial growth by the innate immune response was
needed to recapitulate the threshold effect associated with the effective bacterial inoculum. This
phenomenon could be well captured using our model (Eq. 3), with estimates of δ<sub>M</sub>, the maximal
decay induced by the immune response, of 6.8 log<sub>10</sub> CFUeq/g/h and B<sub>50</sub>, the bacterial load
needed to reach δ<sub>M</sub>/2, of 7.4 log<sub>10</sub> CFUeq/g. This corresponds to a threshold for bacterial
control equal to 6.94 log<sub>10</sub> CFUeq/g. Thus, if the administered bacterial inoculum, B<sub>0</sub>, was
below this value, bacteria were rapidly eliminated by the immune system with a rate ~ δ<sub>M</sub>/B<sub>50</sub>
= 0.3 h<sup>-1</sup>, corresponding to a half-life of 2.3 hours. However, when the bacterial inoculum was
above this value, the model predicts a saturation of the immune response, enabling bacterial
growth with a bacterium doubling time estimated at 3.5 hours.

We performed *in vitro* experiments from which we estimated two key parameters of phage
dynamics that could not be identified *in vivo*, namely the median eclipse phase and the burst
size (see methods and Supplementary Figure S2). Assuming that these parameters had the same
value *in vivo*, we could estimate the infection rate of bacteria by phages, noted *β*, to 6.9 log₁₀
g/CFU/h, and the elimination rate of productively infected bacteria, *δ*_v, to 0.09 h⁻¹ (eg, a half-
life of 7.7 h).

Finally, we also used the model to assess the impact of bacterial kinetics on survival. Assuming
a log-linear effect of bacterial load on the instantaneous risk of death, h(t), the bacterial load
was found to be highly associated with death (γ=1.44, p<10⁻³, Wald test). This corresponds to
a hazard ratio of 4.2; in other words, an increase of 2 log of the bacterial load raised the
instantaneous risk of death by ~18. The joint model developed could simultaneously
recapitulate bacterial kinetics, phage counts and survival rate (Supplementary Figures S5-S6),
despite the heterogeneity of the individual profiles experimentally observed (Figure 5).

<figure><img src="image_30.png" /><figcaption>**Figure 5.** The mathematical model recapitulates the individual variations of the pulmonary bacterial loads and estimates the survival of infected animals.</figcaption></figure>

Data and model predictions for three mice (lines) randomly chosen from each experimental
groups (columns) are shown. The circles correspond to the pulmonary bacterial loads, the color
lines correspond to the individual predictions from the model, and the green line represents the
predicted survival probability. Lines are solid between t=0 and the last observation time point
(death or pre-specified sacrifice) and are dotted after the pre-specified sacrifice. Interrupted
lines indicate the death of animals before the pre-specified time point. The horizontal dashed
line represents the limit of quantification (LOQ=4.7 log₁₀ CFUeq/g).

Effects of bacterial inoculum, dose, route of administration and phage characteristics on
the bacterial kinetics

The model was then used to identify the key factors governing bacterial kinetics.
Unsurprisingly, the bacterial inoculum, B₀, was the most important parameter. In animals with
B₀<6 log₁₀ CFUeq/g, bacterial loads rapidly decreased below the limit of quantification within
48h, regardless of treatment. In animals with large bacterial inoculum B₀>8 log₁₀ CFUeq/g,
phage therapy was not sufficient to control the continuous growth of bacteria, regardless of the
route and administered dose (Figure 5, Figure 6). The benefit of phage therapy was visible
mostly for animals having an intermediate effective inoculum size, e.g. 6<B₀<8 log₁₀ CFUeq/g,
in which phages were critical to prevent bacteria growth and to maintain the bacterial load under
control by the immune system. In those animals, the decline of bacterial load was biphasic, with a rapid decline in the first four hours of treatment, followed by a slower decline afterwards. This slower second phase was attributed in our model to the unrestricted replication of a sub-population of bacteria that were not eliminated by phages, representing both phage-resistant and phage-inaccessible (spatially out of reach) bacteria. The initial proportion of this phage-refractory population was estimated to 10% (Table 1), but could be variable across individuals (from 0.03% to 76%). The dynamics of refractory bacteria was then individual dependent, leading either to a slower but continuous bacterial decline or to a rebound of fully refractory bacteria within 24-48 hours (see individual examples in Figure 5). Importantly, the need to incorporate a compartment of refractory bacteria in the model was strongly supported by data and led to a significant improvement of the data fitting (p<10⁻³, Likelihood Ratio Test, not shown).

<figure><img src="image_31.png" /><figcaption>**Figure 6. Simulations of different bacterial inoculums inform on the impact of refractory bacteria**</figcaption></figure>

Median profile assuming initial bacterial inoculums of 6, 7 and 8 log₁₀ CFUeq/g (top, middle and bottom lines, respectively). Predictions were calculated for each bacterial population (S, R, I₁ and I₂) and phages (V). The colors are indicative of both the route of administration (IT or IV) and the dose of phage received (Lo: 5.6 log₁₀ PFU; Med: 7.6 log₁₀ PFU; Hi: 8.6 log₁₀ PFU): untreated (black), IV-Lo (orange), IV-Med (red), IT-Lo (cyan), IT-med (blue), IT-high (dark blue).
The raw data did not allow to observe a role of route of administration in treatment outcome,
except in the case of a low bacterial inoculum ($B_0=6 \log_{10} \text{CFUeq/g}$) in which the rapid
clearance of bacteria limited the encounter rate with phage. Our simulations indicated that IT
and IV otherwise reached similarly high phages concentrations in lungs (Figure 6). IT
administration was however associated with a substantially faster time to peak viral load, and
this difference was particularly obvious with a low or an intermediate bacterial inoculum. For
instance, with an intermediate bacterial inoculum ($B_0=7 \log_{10} \text{CFUeq/g}$), the time to peak viral
load varied between 1.6 h (IT high dose) and 6.0 h (IV low dose). As a consequence of this
delay in phage maximal replication, IV route had a less efficiency to stop bacterial growth, with
a time to peak bacterial load that could be as large as 4.0 h (IV low dose), to less than 1.1 h with
IT route, and as low as 0.3 h in case of IT high dose.

Finally, we also used the model to explore the effect of intrinsic phage characteristics, namely
burst size and lysis rate, on the bacterial kinetics. A 10-fold reduction in burst size (i.e. 57
PFU/CFU) predicted a larger impact of the route of administration, with a time to bacterial peak
of 10 h post-IV compared to less than 3 h post-IT administration, even when using a low dose
(Figure S7). However, a higher lysis rate value counter-balanced the effects of a low burst size
value (Figure S8). These showed the sensitivity of the model to phage characteristics.
## Discussion

Despite many encouraging data from both ancient literature and recent compassionate treatments, phage therapy remains underexploited, notably because of ambiguous results from clinical trials (Luong et al., 2020). Indeed, the translation from *in vitro* to clinical conditions remains unmet, in particular because using phages as a drug requires a specific method to elaborate a rational scheme for safe and optimal administration. While some studies performed *in silico* or in laboratory conditions (Cairns et al., 2009; Kasman et al., 2002; Leung and Weitz, 2017; Payne and Jansen, 2001, 2003) led to propose several mathematical models of phage-bacteria interactions, very few studies including the host immune response were conducted *in vivo* (Dabrowska and Abedon, 2019; Roach et al., 2017). Here, we developed a pharmacometric model, exploiting *in vivo* data from a murine pneumonia treated by a single virulent phage administered at three doses and via two administration routes, which robustness could support further the (re)deployment of phage therapy.

Our experimental approach was inspired from the work performed during the development of new drugs, where intensive pharmacokinetic and pharmacodynamic data are required to optimize the dosing regimens of a treatment (Marshall et al., 2019). Here, we collected an unprecedented biological set of data using an *in vivo* imaging tool to longitudinally follow the same animals over time and from which we developed a semi-mechanistic model to characterize the dynamics of a phage treatment.

The biodistribution of phages in the lungs of uninfected mice showed that phage loads were much higher and their exposition lasted longer following local (IT) compared to systemic (IV) administration. Consequently, IT administration led to a rapid decline of bacterial loads within the next two hours post-administration. Nevertheless, 12 hours post local and systemic phage administrations, all dosing regimens led to comparable levels of phage exposure, due to the large burst size of phage 536_P1. Despite being in agreement with previous literature (Morello et al., 2011; Oechslin et al., 2017), these biodistribution data cannot be yet generalized to all phages given the very large diversity of these viruses. Also, we did not study repeated IV administered doses that could allow for longer and higher phage exposition in lungs.

The heterogeneity of the initial bacterial loads observed two hours post-infection explains the lack of statistical differences in Kaplan-Meir survival estimations when pooling mice to compare routes of administration and doses. The source of this heterogeneity was the
consequence of choices made to maximize the probability of interactions between phages and bacteria, compared to our previous work using the same couple of phage-bacteria (Dufour et al., 2015). Here, we increased the inoculum and we used the direct intratracheal route (instead of the intranasal route) to reach the lung, in order to reduce spatial dispersion along the upper respiratory tract. The large range of situations observed was instrumental to enrich the mathematical model with multiple events (N=236 analyzed mice) to improve the assessments of phage-bacteria parameters, doses and routes of administration. Moreover, a direct consequence of the above conditions was the precise quantification of the threshold value (6.94 $log_{10}$ CFUeq/g) above which the immune response was overwhelmed, which was very close to the value (6.80 $log_{10}$ CFUeq/g) deduced from our model.

The modeling of phage-bacteria interactions in the host required the definition of a refractory bacterial population. This compartment could not be phenotypically characterized but encompasses bacteria that are escaping phage predation because they are phage-resistant and/or they remain spatially inaccessible to phages (Leung and Weitz, 2017). The inclusion of this compartment is supported by both the emergence of phage resistance observed *in vitro* after 10 h of incubation with phage 536_P1 (Supplementary Figure S2), and previous experimental and theoretical studies showing that the spatial distribution of bacteria in organs can prevent bacteria to be infected by phages (Lourenco et al., 2020; Sousa and Rocha, 2019). This was also in agreement with the modelling of phage therapy in immunocompromised hosts (Roach et al., 2017). Despite the large variation of the bacterial loads that directly affects the initial density of this refractory population, our model could appropriately recapitulate the kinetics of phage-bacteria interactions for each individual mouse, including the survival probability, testifying of its robustness. In addition, simulations performed with larger or lower burst sizes and lysis rates, showed that the model predicts outputs that are in agreement with results from animal studies in which different phages were used to treat bacterial infections (Alemayehu et al., 2012; Carmody et al., 2010; Henry et al., 2013; Morello et al., 2011).

To ensure the model recapitulates kinetics of phage-bacteria interactions during therapy, it was necessary to use as much as possible parameters deduced from animal experiments instead of *in vitro* conditions. For instance, growth rate of bacteria, infectivity rate and lysis rate, were all deduced from experimental data. Nevertheless, other parameters should be more deeply investigated to improve the fidelity of the model. For instance, growth of refractory population was approximated as similar to the susceptible bacteria, neglecting the possibility that the
fitness of mutants might be affected in the host. Indeed, such mutants could be, either or both,
less virulent and strongly affected in their growth, allowing immune cells to eliminate them
more effectively (Oechslin et al., 2017; Pouillot et al., 2012; Smith et al., 1987). The
involvement of the immune response in the success of respiratory phage therapy was previously
demonstrated and modelled, but remains mechanistically poorly explored (Roach et al., 2017).
In this work we chose to account for the immune response in a synthetic way which
encompasses mechanical factors, such as mucins and surfactant proteins, as well as various
cellular actors of the innate immune system (Bergeron et al., 1998; Marriott and Dockrell,
2007).

Predictions from the mathematical model showed the importance of viral load, burst size and lysis rate to fully capture bacterial kinetics and animal outcome. Therefore, in a clinical context, the administration route, and to a lesser degree the dose, will play a critical role in the rate at which bacterial growth will be controlled, supporting the determination of key parameters of phage-bacteria interactions as a prerequisite for optimizing phage therapy. However, the model predicted that substantial variations of the intrinsic characteristics of phages (burst size and lysis rate) were not jeopardizing treatment output, which fits with the current semi-empirical choice of phages for compassionate treatments. Indeed, currently phages are chosen exclusively for their capacity to infect bacteria on *in vitro* conditions where phage-bacteria interactions are artificially optimized. The model is therefore robust, which argues in favor of its relevance for many other phages.

To conclude, the model developed during this work could be used to predict the efficacy of
virtually any phage for which a minimum set of *in vitro* and *in vivo* data must be obtained,
which will considerably lower the number of experiments needed to validate such a phage
during pre-clinical development. Beyond phage therapy, the model could also be implemented
to test combined anti-infectious therapies such as the association of phages with antibiotics (Lin
et al., 2020).
# Methods

## Determination of in vitro infection parameters of phage 536_P1

Parameters of the interactions between strain 536 and phage 536_P1 have been determined *in vitro*, notably the rate of transition *e* from infected non-productive bacteria to productively infected bacteria (eg, log(2)/e is the median eclipse phase duration), and the burst size *b*, i.e. the number of virions released by a lysed bacteria. The median eclipse phase time and the burst size were estimated from non-linear mixed effect models performed with the SAEM algorithm implemented in Monolix software (version 2018R2) (Kuhn and Lavielle, 2005) using data from one-step-growth experiments (see Supplementary material).

## Presentation of the murine model of pulmonary phage therapy

The *E. coli* strain 536 is a virulent clinical strain belonging to the B2 phylogenetic group (Brzuszkiewicz et al., 2006). The *E. coli* strain 536-lux was obtained by the transformation of wild-type 536 strain with plasmid pCM17, which carries the *luxCDABE* operon from *Photorhabdus luminescens* under *OmpC* promoter (Morin and Kaper, 2009). This plasmid is stable over 7 days without adding antibiotic in the medium (Dufour et al., 2015).

Phage 536_P1 (Dufour et al., 2015) is a virulent *Myoviridae* (149.4 Kb) that was amplified on strain 536 and purified according to molecular biology protocols including an endotoxin removal step (Endotrap Blue, Hyglos, Germany). The stock solution of 8.2 log₁₀ PFU/mL has an endotoxin concentration of 2.15 EU/mL.

Pneumonia was induced by intratracheal (IT) administration of strain 536-lux in eight-week-old BALB/cJRj anesthetized male mice (mixture of Ketamine and Xylazine administered intraperitoneally, IP). Exponentially growing bacteria in Luria Bertoni Broth (Lennox) at 37°C were centrifuged and washed twice and finally resuspended in PBS to obtain inoculums ranging from 7.6 to 8.2 log₁₀ CFU in 20µL.

Phage treatment was applied under general anesthesia (isoflurane, inspired fraction=2%) using IT or intravenous (IV, retro-orbital injection) administration two hours after infection. A volume of 20µL of phage 536_P1 was used, containing either a low (Lo; 6.6 log₁₀ PFU), medium (Me; 7.6 log₁₀ PFU) or high (Hi; 8.6 log₁₀ PFU) amount of phages.

At 2, 4, 8, 24 and 48h after phage administration, mice were sacrificed, and blood (in heparin coated tubes) and lungs were recovered and put on ice. Lungs were weighted, homogenized in PBS with protease inhibitor (EDTA free micro Tabs, SigmaAldrich). Lungs and blood were
then centrifuged (5000g, 15 min, 4°C). The supernatant was collected, serially diluted and spotted (in triplicate) on agar plates covered with strain 536 to titer phage 536_P1.

## Bioluminescence and longitudinal follow up of bacterial load in lungs

In all our experiments photon emission of the luminescent bacteria was recorded at regular intervals under general anesthesia (isoflurane, inspired fraction=2%) using IVIS Spectrum and analyzed with Living Image software (PerkinElmer, Richmond, CA) as previously described (Debarbieux et al., 2010). Photons emitted within the chest area of each mouse were collected as photon/s/cm²/steradian.

To fully exploit the *in vivo* imaging tool we established the correlation between bioluminescence and CFU count in lungs using data from 17 untreated animals that were euthanized at different time points and for which lung bacteria was quantified by plating (Figure S1). The relationship obtained allowed us to derive from the luminescence signal an equivalent counts of CFU per gram of lung (noted CFUeq/g) that was used to obtain in all infected animals (untreated or treated by phages) the bacterial densities.

## Phage distribution

A dose of 8 log₁₀ PFU of phage 536_P1 was administered in uninfected mice (N=53) using IT (N=16) or IV (N=37) route. At 2, 4, 8, 24 and 48h after phage administration, mice were sacrificed, and blood (in heparin coated tubes) and lungs were recovered and put on ice. Lungs were weighted, homogenized in PBS with protease inhibitor (EDTA free micro Tabs, SigmaAldrich). Lungs and blood were then centrifuged (5000g, 15 min, 4°C). The supernatant was collected, serially diluted and spotted (in triplicate) on agar plates covered with strain 536 to titer phage 536_P1.

We used the following standard pharmacokinetic model to estimate phage biodistribution:

$$ \frac{dV_B}{dt} = -k_a V_B \qquad (1) $$

$$ \frac{dV}{dt} = F k_a V_B - k_e V \qquad (2) $$

Where V and V_B denote phage quantity in lung and blood compartments, respectively, k_e is the elimination rate from lungs, k_a is the transfer rate from blood to lungs, F is the bioavailability in lungs after IV administration. The initial condition of these equations are V_B(0) = dose and
$V(0) = 0$ in case of IV administration and $V_B(0) = 0$ and $V(0) = dose$ in case of IT administration. $V$ is divided by the phage distribution weight ($W$) to obtain the phage concentration in lung compartment. The parameter $W$ was deduced from the extrapolation of the median phage concentration at baseline ($8 \log_{10} \text{PFU/g}$) following the IT administration of a single dose of $8 \log_{10} \text{PFU}$, and $F$ was obtained by performing the ratio of IV/IT area under the curves in the lung. Estimations of the other parameters were performed by naïve pooling of data (i.e. considering all the measures for each group (IV and IT) as if they came from a single individual and performing a nonlinear regression).

## Bacterial kinetics in untreated mice

We characterize the natural bacterial kinetics (i.e., in absence of treatment) in animals infected with doses of $7.6 \log_{10} \text{CFU}$ (34 mice) or $8.2 \log_{10} \text{CFU}$ (12 mice) and the bacterial kinetics obtained by bioluminescence (see above) using the following model (one mouse was not analyzed due to death during bacterial inoculation, i.e. N=45):

$$ \frac{dB}{dt} = \alpha \left(1 - \frac{B}{B_{max}}\right) B - \delta_M \frac{B}{B_{50} + B} \quad (3) $$

Where $\alpha$ is the exponential growth rate of bacteria (eg, $\log(2)/\alpha$ is the doubling time) and $B_{max}$ is the maximal carrying capacity. The initial condition of $B$ was estimated by the parameter $B_0$ ($B(0) = 10^{B_0}$). Consistent with previous work (Dufour et al., 2019), we assumed that in the first 3 hours of infection the immune system had no activity (eg, bacteria could grow exponentially during this period). After this period the elimination of bacteria by the host was modeled by a non-linear saturable term, given by $\delta_M \frac{s}{B_{50}+B}$. The model parameters were estimated in non-linear mixed effect models performed with the SAEM algorithm implemented in Monolix software (version 2018R2) (Kuhn and Lavielle, 2005).

## Phage-bacteria interaction in treated mice

We characterized phage/bacterial interactions using data obtained in 142 mice infected by $7.6 \log_{10} \text{CFU}$ and treated subsequently (four mice were not analyzed due to death during bacterial inoculation, i.e. N=138), using the protocol described above, with a single dose (low, medium or high) of phages administered using IT or IV route (detail of the groups in Figure 1). Bacterial kinetics was derived from bioluminescence (see above) and lungs were collected at 2, 4, 8, 24 or 48h after phage administration to count PFU. We also considered data in untreated mice (see previous section) for the analysis.
The model describing interactions between phage 536_P1 and strain 536 relied on a system of ordinary differential equations (ODE) (Figure 1, Equations 4 to 10, presented below). We combined standard model of viral (Perelson and Guedj, 2015) and bacterial kinetics (Friberg and Guedj, 2020) to tease out the main parameters governing phage-bacteria interaction in mice. In lungs, phages can infect susceptible bacteria (S) following a mass-action term, with rate $\beta$. We also considered a population of refractory bacteria (of initial proportion $f$), noted $R$, that includes both phage-resistant and phage-inaccessible bacteria.

$$ \frac{dS}{dt} = \alpha \left(1 - \frac{S + I_1 + I_2 + R}{B_{max}}\right)S - \delta_M \frac{S}{B_{50} + S + I_1 + I_2 + R} - \beta VS \quad (4) $$

$$ \frac{dI_1}{dt} = \beta VS - eI_1 - \delta_M \frac{I_1}{B_{50} + S + I_1 + I_2 + R} \quad (5) $$

$$ \frac{dI_2}{dt} = eI_1 - \delta_V I_2 - \delta_M \frac{I_2}{B_{50} + S + I_1 + I_2 + R} \quad (6) $$

$$ \frac{dR}{dt} = \alpha \left(1 - \frac{S + I_1 + I_2 + R}{B_{max}}\right)R - \delta_M \frac{R}{B_{50} + S + I_1 + I_2 + R} \quad (7) $$

$$ \frac{dV_B}{dt} = -k_a V_B \quad (8) $$

$$ \frac{dV}{dt} = - \beta VS + b \delta_V I_2 + F k_a V_B - k_e V \quad (9) $$

We considered the following initial condition for this system: $S(0) = 10^{B_0}(1-f)$, $I_1(0) = 0$, $I_2(0) = 0$, $R(0) = 10^{B_0} \times f$, $V_B(0) = 0$, $V(0) = 0$. The dose of phages is added to $V$ or $V_B$ at the phage administration time by IT or IV route, respectively.

The model parameters in bold were estimated in non-linear mixed effect models performed with the SAEM algorithm implemented in Monolix software using data from all infected animals (version 2018R2) (Kuhn and Lavielle, 2005). All other parameters were fixed to the values determined above, either in vitro (e, b), in uninfected treated mice ($k_a$, $F$, $k_e$) or in infected untreated mice ($\alpha$, $B_{max}$, $\delta_M$, $B_{50}$).

### Impact of bacterial kinetics on survival

Finally, we aimed to take into account the potential bias due to the fact that more severely infected animals died earlier, making the population study less and less representative as time goes by. In this approach, called joint models (Desmee et al., 2015), the probability to survive up to time $t$, $S(t)$, is modeled as: $S(t) = e^{-\int_{0}^{t} h(t)dt}$ where
$$h(t) = \frac{\kappa}{\lambda} \left(\frac{t}{\lambda}\right)^{\kappa-1} e^{\gamma \times \log(S(t)+I_1(t)+I_2(t)+R(t))} \quad (10)$$

$h(t)$ is the instantaneous hazard function, $\kappa$ and $\lambda$ are scale parameters characterizing the risk of death in absence of infection, and $\gamma$ quantifies the impact of bacterial load on mortality. Thus, in this model $\gamma = 0$ is indicative of an absence of association between bacterial load and mortality, while $\gamma > 0$ indicates that increased bacterial load is associated with an increased risk of death.

## Model evaluation

Individual fits were investigated to evaluate the different models. For the final model, we also evaluated the model predictions i) for the bacterial kinetics by using prediction corrected-visual predictive checks (Bergstrand et al., 2011), accounting for the risk of death and ii) visual predictive checks of the survival. Simulations for validation methods were made using mlxR package.

## Model based simulations for the impact of treatment procedures and phage parameter values

Model based predictions were performed to investigate various scenarios as follows. For each simulation, the median profile of different model compartments ($S, R, I_1, I_2$ and $V$) was investigated using the parameters determined in Table 1 and different experimental conditions regarding the route of administration (IV or IT), the dose of phage administered (Lo: 5.6, Med: 7.6 or Hi: 8.6 $log_{10}$ PFU) and initial bacterial load (6, 7 or 8 $log_{10}$ CFUeq/g). We also evaluated the effects of different burst size (10 fold lower, 10 fold larger than the estimated value in vitro), using a single initial bacterial load of 7 $log_{10}$ CFUeq/g. Furthermore, when setting a burst size 10 times lower than the estimated value (and a fixed initial bacterial load of 7 $log_{10}$ CFUeq/g), we investigated the effect of three different lysis rates values: the value estimated in vivo (0.09 $h^{-1}$); a value three times higher (0.27 $h^{-1}$); a value estimated in vitro (1.39 $h^{-1}$).

## Ethical statement

A total of 241 eight-week-old BALB/cJRj male mice (Janvier, France) were housed in animal facility in accordance with French and European regulations on the care and protection of laboratory animals. Food and drink were provided ad libitum. Protocol was approved by the veterinary staff of the Institut Pasteur animal facility (approval number 20.173) and the National
Ethics Committee regulating animal experimentation (APAFIS#26874-2020081309052574 v1).
**Authors contribution**

JDR, JG, LD and RD, conceived the study and secured funding. BG, FH, RD and RK performed experiments. JDR, JG, JS, LD and RD analyzed the data. JG, JS and TTN performed modeling. ND provided feedback. JS and RD wrote the first draft. JG and LD reviewed and edited final draft.

**Declaration of interests**

Authors declare no competing interest.
**Supplementary figures**

<figure><img src="image_32.png" /><figcaption>**Supplementary figure S1.** Correlation between bioluminescence and CFU from lungs of infected and untreated mice. The regression line equation is $\log_{10}(\text{CFU/g}) = 1.44 + 1.07 \log_{10}(\text{Bioluminescence})$. From 17 mice, 14 were infected with an inoculum of 7.6 $\log_{10}$ CFU (grey dots) and 3 with an inoculum of 8.2 $\log_{10}$ CFU (black dots).</figcaption></figure>
<figure><img src="image_33.png" /><figcaption>**Supplementary figure S2.** *In vitro* experiments describing the dynamic interactions between strain 536 and phage 536_P1.</figcaption></figure>

The one step growth data (3 replicates) are presented on top left, the dynamic quantification of strain 536 and phage 536_P1 (3 replicates at MOI 0 and at MOI 0.001) are presented on bottom left, and the kinetics of lysis (3 or 4 replicates at MOI 0, MOI 0.0001 and MOI 10) are presented on the right (see Supplementary methods).
<figure><img src="image_34.png" /></figure>

**Supplementary figure S3.** Longitudinal bacterial load data and individual fits from all infected mice (N=183). Solid lines represent individual predictions from the infection time to the last observation time. The dotted lines are individual predictions starting from the sacrifice time point. Interrupted solid lines correspond to animal that die before the scheduled sacrifice time. Circles correspond to live animals while triangles correspond to animals at time of sacrifice after recording the bioluminescence.
<figure><img src="image_35.png" /></figure>

**Supplementary figure S4.** Survival curves over time of *E. coli* strain 536 infected mice and untreated (black, N=45), or treated by phage 536_P1 administered intravenously (IV) at low dose of 5.6 log₁₀ PFU (orange, N=50), medium dose of 7.6 log₁₀ PFU (red, N=29), intratracheally (IT) at low dose (cyan, N=24), medium dose (blue, N=25) or high dose of 8.6 log₁₀ PFU (dark blue, N=10). Crosses represent sacrifice time points.
<figure><img src="image_36.png" /><figcaption>**Supplementary figure S5.** Prediction-corrected visual predictive check plot for the bacterial load, using the final model. The 90% confidence intervals around the 10th (bottom), 50th (middle) and 90th (top) model-predicted percentiles, obtained from 1000 simulated data sets from the model, are represented as blue (10th and 90th percentiles) and pink (50th percentile) areas. The blue dots represent the 10th, 50th and 90th percentiles of the observations from the original data set.</figcaption></figure>
<figure><img src="image_37.png" /></figure>

**Supplementary figure S6.** Visual predictive check plot for the survival in each group, using the final model. Untreated mice are in black, phage 536_P1 treated mice (2 h after infection) groups are in colors, with different phage administration scheme: intravenously (IV) at low dose (5.6 log₁₀ PFU) in orange or medium dose (7.6 log₁₀ PFU) in red, intratracheally (IT) at low dose in cyan, medium dose in blue or high dose (8.6 log₁₀ PFU) in dark blue (IT Hi). The dashed lines around the solid line represent the 10th (bottom), 50th (middle) and 90th (top) percentiles of survival obtained from 1000 simulated data sets from the model. The stepped line is the observed survival (Kaplan-Meier curve), with crosses representing censoring.
<figure><img src="image_38.png" /></figure>

**Supplementary figure S7.** Median profile of model predictions for different burst sizes: 57, 570 (i.e. the estimated value) and 5700 PFU/CFU in the top, middle, and bottom line, respectively), with a fixed inoculum of 7 log₁₀ CFUeq/g. Each column represents different model compartments: bacteria (S, R, I₁ and I₂) and phages in lungs. The colors are indicative of both the route of administration (IT or IV) and the dose of phage received (Lo: 5.6 log₁₀ PFU; Med: 7.6 log₁₀ PFU; Hi: 7.6 log₁₀ PFU): untreated (black), IV-Lo (orange), IV-Med (red), IT-Lo (cyan), IT-med (blue), IT-high (dark blue).
# Supplementary Material

## Purification of bacteriophages

Phage 536_P1 was amplified on strain 536 and purified by first, filtration of lysate to 0.1µm followed by tangential ultrafiltration (Vivaflow™, Sartorius, 100kDa), two ultracentrifugations on cesium chloride, dialysed against TN (Tris-HCl pH8, NaCl 150 mM), passed three times through an endotoxin removal column (Endotrap Blue, Hyglos, Germany) and filtered-sterilized at 0.22µm.

## In vivo experiments

Intratracheal administration (20µL) was performed in anesthetized mice using elongated tips guided by the glottis visualization with a cold light. The correct position of the tip was checked by the liquid movement inside the tip caused by the ventilatory movement of the animal. The liquid was then expulsed from the tip with a pipette, and 100 µL of air was then pushed to aerosolize the solution to the bronchoalveolar space. Mice were maintained in prone position during at least 5 minutes.

## In vitro experiments

### One-Step-Growth experiments

The one-step growth experiment was performed in triplicate, using Lysogeny Broth (LB, Difco Bacto-Tryptone: 10 g/L, Difco Yeast extract Difco 5 g/L, NaCl 5 g/L). As previously described (Hyman and Abedon, 2009), an exponential-phase culture (9 ml) of strain 536 (optical density at 600 nm 0,25, ie 9 log₁₀ CFU) was mixed with 1mL of 8 log₁₀ PFU of phage 536_P1 (MOI of 0.1). The mixture was placed under constant shaking (100 rpm) at 37°C. Duplicated samples were collected at 2-min intervals during 30 minutes. Six minutes after mixing phages and bacteria, the mixture was 1000 fold diluted in LB medium in order to lower bacterial concentration and avoid any reinfection. Out of the duplicated samples, one was immediately diluted and plated for phage titration, while the second was treated with 10% (vol/vol) chloroform to release intracellular phages in order to determine the eclipse period before phage titration.

### Dynamic quantification of bacteria and phages over time

We first established the correlation between optical density and CFU for strain 536 by plating on LB agar dilutions of independent aliquots of exponentially growing cells from which OD600 nm was recorded before dilutions. Then, an exponential bacterial culture (50 mL of 1.10⁷