# Marginally Stable Current Sheets in Collisionless Magnetic Reconnection

C. Granier<EMAIL_ADDRESS>Université Côte d’Azur, CNRS, Observatoire
de la Côte d’Azur, Laboratoire J. L. Lagrange, Boulevard de l’Observatoire, CS
34229, 06304 Nice Cedex 4, France  Istituto dei Sistemi Complessi - CNR and
Dipartimento di Energia, Politecnico di Torino, Torino 10129, Italy  D.
Borgogno  Istituto dei Sistemi Complessi - CNR and Dipartimento di Energia,
Politecnico di Torino, Torino 10129, Italy  L. Comisso  Department of
Astronomy and Columbia Astrophysics Laboratory, Columbia University, New York,
NY 10027, USA  D. Grasso  Istituto dei Sistemi Complessi - CNR and
Dipartimento di Energia, Politecnico di Torino, Torino 10129, Italy  E. Tassi
Université Côte d’Azur, CNRS, Observatoire de la Côte d’Azur, Laboratoire J.
L. Lagrange, Boulevard de l’Observatoire, CS 34229, 06304 Nice Cedex 4, France
R. Numata  Graduate School of Information Science, University of Hyogo, Kobe
650-0047, Japan

###### Abstract

Non-collisional current sheets that form during the nonlinear development of
magnetic reconnection are characterized by a small thickness, of the order of
the electron skin depth. They can become unstable to the formation of
plasmoids, which allows the magnetic reconnection process to reach high
reconnection rates. In this work, we investigate the marginal stability
conditions for the development of plasmoids when the forming current sheet is
purely collisionless and in the presence of a strong guide field. We analyze
the geometry that characterizes the reconnecting current sheet, and what
promotes its elongation. Once the reconnecting current sheet is formed, we
identify the regimes for which it is plasmoid unstable. Our study shows that
plasmoids can be obtained, in this context, from current sheets with an aspect
ratio much smaller than in the collisional regime, and that the plasma flow
channel of the marginally stable current layers maintains an inverse aspect
ratio of $0.1$.

††preprint: APS/123-QED

Magnetic reconnection is a fundamental plasma process that involves a rapid
topological change of the magnetic field leading to an efficient magnetic
energy conversion. Magnetic reconnection typically occurs via current sheets
(CS), where non-ideal plasma effects become important allowing the change of
magnetic field line connectivity Zweibel and Yamada (2009); Yamada  _et al._
(2010); Ji  _et al._ (2022).

It is well established that the instabilities of thin CS, that lead to the
formation of plasmoids, have a fundamental impact on the reconnection rate
Daughton  _et al._ (2009); Bhattacharjee  _et al._ (2009). Indeed, even in the
resistive magnetohydrodynamics (MHD) framework, the development of plasmoids
in the reconnection layer induces a fast magnetic reconnection regime
characterized by a reconnection rate that can exceed the estimates based on
the Sweet-Parker (SP) theory Sweet (1958); Parker (1957) by several orders of
magnitude. In collisional CS, it has been shown that plasmoids develop when
the Lundquist number $S=L_{\rm{cs}}v_{A}/(\eta c^{2}/4\pi)$ exceeds the
threshold value $S_{\star}\sim 10^{4}$ Biskamp (1986). Here, the Lundquist
number is defined with the length of the CS, $L_{\rm{cs}}$. The other
quantities are the plasma resistivity $\eta$, the speed of light $c$, and the
Alfvén speed $v_{A}$. The threshold value on the Lundquist number,
$S_{\star}$, separates the Sweet-Parker regime from the plasmoid-mediated
regime of collisional reconnection. In addition, it controls the reconnection
rate $R_{\rm{rec}}$ in the plasmoid-mediated regime, $R_{\rm{rec}}\sim
S_{\star}^{-1/2}v_{A}B_{\rm{up}}$ Huang and Bhattacharjee (2010); Uzdensky
_et al._ (2010); Comisso  _et al._ (2015a); Comisso and Grasso (2016), where
$B_{\rm{up}}$ is the reconnecting magnetic field. The extension of the
resistive reconnection regime with the inclusion of the ion dynamics
associated with the ion sound Larmor radius, $\rho_{s}$, or the ion inertial
length, $d_{i}$, complicates the picture. Indeed, when the thickness of the
reconnecting CS shrinks below $\rho_{s}$ (for reconnection with a guide field)
or $d_{i}$ (for reconnection without guide field), the process becomes even
faster and approaches $R_{\rm{rec}}\sim 0.1v_{A}B_{\rm{up}}$ Comisso and
Bhattacharjee (2016). Given its importance, the transition between the
different regimes of reconnection has been thoroughly investigated, and the
current understanding of reconnection driven by plasma resistivity has been
summarized in the form of parameter space diagrams Ji and Daughton (2011);
Daughton and Roytershteyn (2012); Huang and Bhattacharjee (2013); Karimabadi
and Lazarian (2013); Comisso  _et al._ (2015b); Le  _et al._ (2015); Loureiro
and Uzdensky (2015); Bhat and Loureiro (2018).

In contrast, the marginal stability of reconnecting CS in the collisionless
regime has seen relatively little investigation. This subject was approached
in Ji and Daughton (2011), in which it is argued that, below the scales
$d_{i}$ or $\rho_{s}$, no plasmoids were formed. Yet, it is acknowledged that
reconnection in nature is often driven by collisionless effects beyond the
resistive MHD description.

In this Letter, we investigate a phase space described by the two kinetic
scales $d_{e}$ (electron inertial length) and $\rho_{s}$, compared to the
current length $L_{\rm{cs}}$. We show how the aspect ratio of the marginally
stable reconnection layer depends on these relevant kinetic scales. We believe
this study might also be useful to support observational and experimental
results. In particular, recent observations revealed many reconnection onsets
driven by electrons, in the presence of a strong guide field, close to the
dayside magnetopause and magnetosheath Burch  _et al._ (2016); Phan  _et al._
(2018). Moreover, in Ref. Phan  _et al._ (2018), current sheets having a
thickness of the order of the electron inertial length were identified. A
study also gave direct experimental proof of plasmoid formation at the X point
and at the electron scale in a regime where no plasmoids were predicted by the
theory Olson  _et al._ (2016).

We assume a plasma immersed in a strong (guide) magnetic field of amplitude
$B_{0}$, resulting in low plasma $\beta$ (the ratio of plasma pressure to
magnetic pressure). In order to reduce the problem to a few essential
ingredients, in our analysis we adopt a simple two-fluid model that retains
electron inertia effects, as well as ion sound Larmor radius effects.
Specifically, the equations governing the plasma dynamics are (e.g. Cafaro
_et al._ , 1998)

$\frac{\partial n_{e}}{\partial t}+[\phi,n_{e}]=[A_{\parallel},u_{e}]\,,$ (1)
$\frac{\partial}{\partial
t}\left(A_{\parallel}-d_{e}^{2}u_{e}\right)+\left[\phi,A_{\parallel}-d_{e}^{2}u_{e}\right]=\rho_{s}^{2}[n_{e},A_{\parallel}]\,,$
(2)

where $A_{\parallel}$ and $\phi$ are the magnetic and electrostatic potential,
$n_{e}=\nabla_{\perp}^{2}\phi$ is the electron density perturbation, and
$u_{e}=\nabla_{\perp}^{2}A_{\parallel}$ is the parallel electron velocity,
also proportional to the current density. The variables are normalized as
$\\{t,x,A_{\parallel},\phi\\}=\\{{v_{A}}\hat{t}/L,{\hat{x}}/{L},\hat{A}_{\parallel}/(LB_{0}),c\hat{\phi}/(v_{A}LB_{0})\\}$,
where the caret ($\,\hat{}\,$) indicates dimensional quantities and with $L$
the characteristic equilibrium scale length set by the equilibrium magnetic
field. The normalized magnetic field and the perpendicular
$\bm{E}\times\bm{B}$ velocity are related to $A_{\parallel}$ and $\phi$ as
${\bm{B}}={\bm{\hat{z}}}+\nabla A_{\parallel}\times{\bm{\hat{z}}}$ and
${\bm{u}}_{\perp}={\bm{\hat{z}}}\times\nabla\phi$, respectively. The
parameters are $d_{e}=\sqrt{m_{e}c^{2}/{4\pi e^{2}n_{0}}}/L$ and
$\rho_{s}=\sqrt{{T_{e}m_{i}c^{2}}/{e^{2}B_{0}^{2}}}/L$, corresponding to the
normalized electron skin depth and ion sound Larmor radius, respectively. Ions
are assumed to be cold. In Eqs. (1)-(2),
$[f,g]=\partial_{x}f\partial_{y}g-\partial_{y}f\partial_{x}g$.

In order to analyze the marginal stability conditions of the plasmoid
instability in the collisionless regime, we conducted a large number ($\sim
30$) of numerical simulations of the system of Eqs. (1)-(2). The numerical
solver is pseudo-spectral and the advancement in time is done through a third
order Adams–Bashforth scheme. We considered a periodic 2D domain $2L_{x}\times
2L_{y}$, resolved with a number of grid points up to $2000\times 2400$. On the
other hand, in the figures presented in the paper, only a part of the
computational domain is shown. We set up an initial equilibrium given by
$\phi^{(0)}(x)=0$, $A_{\parallel}^{(0)}(x)=1/\cosh^{2}\left(x\right)$. The
tearing stability parameter for this equilibrium is
$\Delta^{\prime}_{\rm{box,m}}=2[(5-k_{y}^{2})(k_{y}^{2}+3)]/[k_{y}^{2}(k_{y}^{2}+4)^{1/2}]$.
This equilibrium is tearing unstable if $\Delta_{\rm{box,m}}^{\prime}>0$, thus
for a wavenumber $k_{y}=\pi m/L_{y}<\sqrt{5}$. We will always refer to
$\Delta^{\prime}_{\rm{box}}$ as being associated to the mode $m=1$, and we
change its value by taking different box lengths along the $y$ direction. With
this set up, one or several tearing modes are initially unstable. The dominant
mode generates two magnetic islands separated by a reconnection point
(X-point). During the nonlinear phase, a slowly thinning CS forms self-
consistently at the X-point location. The evolution of this current sheet will
be decisive for the formation of plasmoids.

In the following, we characterize this reconnecting CS according to the
parameters $d_{e}$, $\rho_{s}$, and $\Delta^{\prime}_{\rm{box}}$.
Specifically, we measured the length and the width of the CS at a time $t$
just before the plasmoid onset. We define the measure $L_{\rm{cs}}$, such
that, taking the variation from the highest current position $u_{e}|_{X}$
($u_{e}$ evaluated at the $X$-point), the standard deviation of the current
distribution from $y=0$ to $y=L_{\rm{cs}}/2$ equals unity, I.e.
$\sqrt{{\sum^{N}_{i=1}\left[u_{e}|_{X}-u_{e}(0,i\Delta
y,t)\right]^{2}}/{N}}=1$, where $\Delta y$ is the distance between two grid-
points along $y$ and $N$ indicates the number of points from $y=0$ to
$y=L_{\rm{cs}}/2$. This method makes it possible to account for the decrease
of the current intensity along the layer. Once $L_{\rm{cs}}$ is identified,
the half width of the CS, which we denote by $\delta_{\rm{cs}}/2$, corresponds
to the distance, along $x$, between $u_{e}|_{X}$ and the position where the
current reaches the value
$u_{e}(\delta_{\rm{cs}}/2,0)=u_{e}(0,L_{\rm{cs}}/2)$. We also measure the
width and length of the outflow velocity channel coming out from the end of
the CS. The length $L_{\rm{outf}}$ corresponds to the distance between the
upward and downward peaks in the distribution of $u_{y}=\partial_{x}\phi$.
while the width $\delta_{\rm{outf}}$ is also measured with the standard
deviation method. The aspect ratios $A_{\rm{cs}}=L_{\rm{cs}}/\delta_{\rm{cs}}$
and $A_{\rm{outf}}=L_{\rm{outf}}/\delta_{\rm{outf}}$ are also reported.

We first focus on the limit $\rho_{s}=0$ shown in Fig. 1. As discussed in Ref.
Ottaviani and Porcelli (1993), in this limit $\delta_{\rm{cs}}\propto d_{e}$.
To better show the geometry, the colored contour maps of $u_{e}$ with
superimposed contour lines of $A_{\parallel}$ in black, are shown for certain
cases on Fig. 1 as well as on Fig. 2. For low values of $d_{e}$, a high and
uniform current density allows the parallel alignment of a high density of
magnetic field lines (see color map of $u_{e}$ for $d_{e}=0.05$ and
$\rho_{s}=0$). On the other hand, for high $d_{e}$ values, the current is not
uniform enough along the layer for the magnetic field lines to line up
perfectly, since their density decreases in the region where the current is
weaker, (see $d_{e}=0.6$ and $\rho_{s}=0$). As we discuss below, this latter
case is less likely to develop plasmoids. Finally, in the limit $\rho_{s}=0$,
we obtain the approximate scalings $L_{cs}\propto d_{e}^{-1/10}$ and
$A_{\rm{cs}}\propto d_{e}^{-1}$.

Figure 1: Characteristics of the reconnecting CS as a function of $d_{e}$ for
fixed $\rho_{s}=0$ and $\Delta^{\prime}_{\rm{box}}=60$. The CS is unstable to
the formation of plasmoids in all five cases. For $d_{e}=0.05$ and
$d_{e}=0.6$, we show the color maps of $u_{e}$ with isolines of
$A_{\parallel}$ in black.

When $\rho_{s}$ is taken into account, ion sound Larmor effects can become
important and the CS changes into a cross shaped structure aligned with the
magnetic island separatrices Cafaro  _et al._ (1998). Indeed, in Fig. 2, when
$\rho_{s}$ is increased (for $\rho_{s}\sim d_{e}$), a part at the end of the
layer splits to extend along the separatrices (see $d_{e}=0.05$ and
$\rho_{s}=0.05$). Here, the measured $L_{\rm{cs}}$ still corresponds to the
length distributed symmetrically on both sides of $y=0$. We measured
$L_{\rm{cs}}\propto\rho_{s}^{-1/2}$. As for the aspect ratios, they scale as
$A_{\rm{cs}}\propto\rho_{s}^{-0.6}$ and $A_{\rm{outf}}\propto\rho_{s}^{-1/2}$.
For the series of simulations with $\Delta^{\prime}_{\rm{box}}=14.3$ and
$d_{e}=0.1$, the reconnection process occurs without forming any plasmoids
(gridded red region) until $\rho_{s}\sim 0.4$.

For $\rho_{s}\gtrsim d_{e}$ and $\Delta^{\prime}_{\rm{box}}=60$ (green
diagonally striped region), the aspect ratio $A_{\rm{cs}}$ is sufficiently
large and one plasmoid emerges from the center of the CS. This corresponds to
a low wavenumber fluctuation that develops in the CS, which is entering the
nonlinear phase.

For $\rho_{s}\gg d_{e}$ (green dotted region), the CS reaches a perfect cross
shape Cafaro  _et al._ (1998). This very different geometry can still lead to
a more complex plasmoid formation. Indeed, in the regime $\rho_{s}\gg d_{e}$,
the first plasmoids that break up the CS are symmetrically located above and
below the $X$-point. This process is detailed in Fig. 3. We observe 4 main
phases: (I) formation of the $X$-shaped current, (II) its ends meet to form a
local $Y$-shaped CS, (III) plasmoids emerge and enter the nonlinear phase,
(IV) they are expelled by the outflow and the center plasmoid emerges. This
type of plasmoid onset takes place for $\rho_{s}>0.4\gg d_{e}$ in Fig. 2.

Figure 2:  Same measures as Fig. 1 but for fixed $d_{e}$ and varying
$\rho_{s}$. For $\Delta^{\prime}_{\rm{box}}=60$ and $d_{e}=0.05$, all five
cases are plasmoid unstable. For $\Delta^{\prime}_{\rm{box}}=14.3$ and
$d_{e}=0.1$, plasmoids grow only when $\rho_{s}\geq 0.4$. For $\rho_{s}=0.05$
and $\rho_{s}=0.5$, we show the color maps of $u_{e}$ with isolines of
$A_{\parallel}$ in black. Figure 3: Time evolution of the parallel current
density. The plots show the color maps of $u_{e}$, while black lines are
contour lines of $A_{\parallel}$. For this simulation $d_{e}=0.085$,
$\rho_{s}=0.3$ and $\Delta^{\prime}_{\rm{box}}=30$.

We now discuss the dependence on the $\Delta^{\prime}_{\rm{box}}$ parameter,
for $\rho_{s}=0$ and for $\rho_{s}\gg d_{e}$. In order to clearly identify a
CS, we have considered large $\Delta^{\prime}_{\rm{box}}$ values, which vary
from $11.3$ to $240$. For $\rho_{s}=0$ (Fig. 4), the $L_{\rm{cs}}$ depends
linearly on $\Delta^{\prime}_{\rm{box}}$, as in the resistive case Waelbroeck
(1989); Jemella  _et al._ (2003, 2004); Loureiro  _et al._ (2005). We do not
obtain plasmoids for $\Delta^{\prime}_{\rm{box}}\leq 14.3$, in agreement with
Del Sarto  _et al._ (2003) where this regime is shown to be prone to the
development of the Kelvin-Helmholtz instability. In the cases with
$21\leq\Delta^{\prime}_{\rm{box}}\leq 38$, one plasmoid emerges and breaks up
the reconnecting CS. For $\Delta^{\prime}_{\rm{box}}=240$, two other plasmoids
are formed when the reconnecting CS becomes more elongated (unstable) as
$\Delta^{\prime}_{\rm{box}}$ increases. In the limit $\rho_{s}=0$, the outflow
channel follows the CS and we observe indeed the scaling
$A_{\rm{outf}}\propto\Delta^{\prime}_{\rm{box}}$ (not shown here).

For $\rho_{s}\gg d_{e}$ (Fig. 5), on the other hand, the case with
$\Delta^{\prime}_{\rm{box}}=14.3$ is plasmoid unstable. In this regime, the
small-scale, oscillating current layer pattern located inside the two magnetic
islands, identified in Refs. Grasso  _et al._ (2001); Del Sarto  _et al._
(2003), is visible on the two left panels. In the rightmost panel of Fig. 5,
we show the measured aspect ratio of the outflow velocity channel just before
the appearance of the first plasmoid. For the least unstable reconnecting CS
($\Delta^{\prime}_{\rm{box}}=14.3)$, we measured $L_{\rm{outf}}=1.21$ and
$\delta_{\rm{outf}}=0.12$, which implies a steady state reconnection rate of
$R_{\rm{rec}}\sim(\delta_{\rm{outf}}/L_{\rm{outf}})v_{A}B_{\rm{up}}\sim
0.1v_{A}B_{\rm{up}}$. The red area corresponds to stable cases. The green
striped area corresponds to the onset of only one plasmoid located at the
center of the CS. Finally, the green dotted region corresponds to the cases
where the first plasmoids emerge from a local $Y$-shaped CS (as described on
Fig. 3).

We can construct a parameter space diagram (Fig. 6), analogously to what was
done for reconnection induced by plasma resistivity Ji and Daughton (2011);
Daughton and Roytershteyn (2012); Cassak and Drake (2013), which allows one to
identify the collisionless plasmoid regimes that take place once a
reconnecting layer of a certain length is formed. According to our numerical
simulations, the critical aspect ratio above which plasmoids break up the
reconnecting CS is
$A^{(1)}_{\star}=(L_{\rm{cs}}/\delta_{\rm{cs}})^{(1)}_{\star}\sim(L_{\rm{cs}}/d_{e})^{(1)}_{\star}\sim
10$ when $d_{e}\gg\rho_{s}$. On the other hand, for $\rho_{s}\gtrsim d_{e}$,
the plasmoid formation has a different threshold
$A^{(2)}_{\star}=(L_{\rm{cs}}/\delta_{\rm{cs}})^{(2)}_{\star}$, and the
simulations indicate $A^{(2)}_{\star}<A^{(1)}_{\star}$.

Figure 4: CS characteristics for varying $\Delta^{\prime}_{\rm{box}}$ at fixed
$d_{e}=0.085$ and $\rho_{s}=0$. From left to right: color maps of $u_{e}$ with
isolines of $A_{\parallel}$ in black for $\Delta^{\prime}_{\rm{box}}=14.3$,
$\Delta^{\prime}_{\rm{box}}=38$, $\Delta^{\prime}_{\rm{box}}=240$ and plot of
the same quantities as in Fig. 1 but here as functions of
$\Delta^{\prime}_{\rm{box}}$. Figure 5: CS characteristics for varying
$\Delta^{\prime}_{\rm{box}}$ at fixed $d_{e}=0.085$ and $\rho_{s}=0.3$. From
left to right: color maps of $u_{e}$ with isolines of $A_{\parallel}$ in black
for simulations with $\Delta^{\prime}_{\rm{box}}=14.3$,
$\Delta^{\prime}_{\rm{box}}=38$, $\Delta^{\prime}_{\rm{box}}=240$, and aspect
ratio of the outflow velocity channel before the onset of plasmoids.

We can evaluate $A^{(2)}_{\star}$ by taking inspiration from the plasmoid
instability theory presented in Refs. Comisso  _et al._ (2016, 2017); Huang
_et al._ (2017). We consider a forming CS in which the amplitude of the
tearing mode grows as
$A_{\parallel}(k,t)=A_{0}\mbox{exp}(\int_{t_{0}}^{t}\gamma(k)\mbox{d}t^{\prime})$,
where $\gamma$ and $k$ are the tearing mode growth rate and wavenumber,
respectively, while $A_{0}$ is the magnetic flux amplitude at $t_{0}$. The
plasmoid half-width is given by $w(k,t)=2(A_{\parallel}a/B_{\rm{up}})^{1/2}$,
where $A_{\parallel}$ is evaluated at the resonant surface and $a$ is the
half-width of the CS. We verified that the CS profiles are well fitted by a
Harris sheet Biskamp (2000), for which
$\Delta^{\prime}_{\rm{cs}}=2[(ka)^{-1}-ka]/a$. Given that the CS is slowly
shrinking toward a finite width, we assume that, just before the plasmoid
onset, the CS is in nearly steady-state and we neglect its time dependence.
From Eqs. (1)-(2), one can derive the dispersion relation of the collisionless
tearing mode for arbitrary values of $\Delta^{\prime}_{\rm{cs}}$ Porcelli
(1991). For the marginally stable CS, one can consider the limit
$\delta_{\rm{in}}\Delta^{\prime}_{\rm{cs}}\ll 1$, with $\delta_{\rm{in}}$
indicating the width of the inner tearing layer. In this case, the full
dispersion relation Porcelli (1991) reduces to
$\gamma^{(1)}=\left[{\Gamma(1/4)}/{2\pi\Gamma(3/4)}\right]^{2}\Delta_{\rm{cs}}^{{}^{\prime}2}d_{e}^{3}k/a$
for $\rho_{s}^{2}\ll d_{e}^{2}$, while
$\gamma^{(2)}={\Delta^{\prime}_{\rm{cs}}d_{e}\rho_{s}k}/{a\pi}$ when
$\rho_{s}^{2}\gg d_{e}^{2}$.

Figure 6: Parameter space diagram that identifies where the plasmoid
instability leads to the break up of the CS. The cases that give plasmoids
(green dots) are either above the $A^{(1)}_{\star}$ or $A^{(2)}_{\star}$
thresholds. The threshold $A^{(2)}_{\star}$ is approximated by the scaling
$L_{{\rm{cs}},\star}^{(2)}/d_{e}\propto L_{\rm{cs}}/\rho_{s}$ from Eq. (4).

We denote by ${\tau}_{\star}=L_{{\rm{cs}},\star}/v_{A}$ the timescale for the
plasma to be expelled from the CS because of the Alfvénic outflow. If the
magnetic flux amplitude becomes nonlinear (with plasmoid half-width
$w_{\rm{nl}}$) in a time shorter than ${\tau}_{\star}$, the CS is broken by at
least one plasmoid. Otherwise it remains stable. Therefore, taking
$w(k,{\tau}_{\star})=2\left({A_{0}a}/{B_{\rm{up}}}\right)^{1/2}e^{\frac{1}{2}{\tau}_{\star}\gamma}$,
the threshold for the plasmoid formation can be written as

${\tau}_{\star}\gamma=2\;\mbox{ln}\left[\frac{w_{\rm{nl}}}{2}\left(\frac{B_{\rm{up}}}{A_{0}a}\right)^{1/2}\right]\,.$
(3)

Assuming that the needed amplification factor of the magnetic flux
perturbation is the same for the $\rho_{s}^{2}\ll d_{e}^{2}$ and
$\rho_{s}^{2}\gg d_{e}^{2}$ cases, requiring
${\tau}_{\star}^{(2)}\gamma^{(2)}\sim{\tau}_{\star}^{(1)}\gamma^{(1)}$, making
use of the numerical result $(L_{\rm{cs}}/d_{e})^{(1)}_{\star}\sim 10$, and
considering that for $(ka)^{2}\ll 1$ we have $\Delta^{\prime}_{\rm{cs}}k\sim
1/a^{2}\sim 1/d_{e}^{2}$, with $k\propto 1/L_{\rm{cs}}$, gives us the
threshold condition

$\frac{L_{{\rm{cs}}}}{d_{e}}=\frac{L_{{\rm{cs}},\star}^{(2)}}{d_{e}}\propto\frac{L_{\rm{cs}}}{\rho_{s}}\,.$
(4)

We identified a proportionality coefficient for which the proposed scaling,
shown by the dashed blue line in Fig. 6, correctly captures the plasmoid
formation that occurs for significantly lower values of the CS aspect ratio
when $L_{\rm{cs}}/\rho_{s}\lesssim 1$.

While the aspect ratio of the CS controls the plasmoid growth, the aspect
ratio of the plasma flow channel regulates the rate of inflowing plasma via
mass conservation. For $d_{e}\gg\rho_{s}$, the aspect ratios $A_{\rm{cs}}$ and
$A_{\rm{outf}}$ essentially coincide since the plasma behaves as a one fluid.
Therefore, for an incompressible flow in steady state, the marginal stability
threshold $A^{(1)}_{\star}\sim 10$ yields the reconnection rate
$R_{\rm{rec}}\sim 0.1v_{A}B_{\rm{up}}$. On the other hand, for
$\rho_{s}\gtrsim d_{e}$, two-fluid effects lead to a decoupling of the plasma
flow channel from the electric current density, and in this case we find that
$R_{\rm{rec}}\sim(\delta_{\rm{outf}}/L_{\rm{outf}})^{(2)}_{\star}v_{A}B_{\rm{up}}\sim
0.1v_{A}B_{\rm{up}}$ even when $A^{(2)}_{\star}\ll A^{(1)}_{\star}$. Since the
global reconnection rate is controlled by the marginally stable CS Comisso and
Bhattacharjee (2016), eventually $R_{\rm{rec}}\sim 0.1v_{A}B_{\rm{up}}$ in the
entire green and blue parameter space regions of Fig. 6.

In summary, we have identified, with two-fluid numerical simulations and
analytical arguments, the marginal stability conditions for the development of
plasmoids in collisionless reconnecting CS. We find that in the collisionless
regime, reconnecting CS are unstable to the formation of plasmoids for
critical aspect ratios that can be as small as
$L_{\rm{cs}}/\delta_{\rm{cs}}\leq 10$. For the marginally stable CS, we find
that the aspect ratio of the outflow channel is
$L_{\rm{outf}}/\delta_{\rm{outf}}\sim 10$ independent of the microscopic
plasma parameters. The space of collisionless plasma parameters
($L_{\rm{cs}}/d_{e}$ and $L_{\rm{cs}}/\rho_{s}$) for which magnetic
reconnection driven by electron inertia occurs in the plasmoid-mediated regime
is organized in a new phase space diagram for collisionless reconnection. A
new phase space diagram spanned by $L_{\rm{cs}}/d_{e}$ and
$L_{\rm{cs}}/\rho_{s}$ for collisionless reconnection is presented, which
extends the diagram of plasmoid onset for collisional plasmas Ji and Daughton
(2011). Our results allow one to separate the collisionless laminar regime of
reconnection from the collisionless plasmoid-mediated regime. The properties
of the marginally stable CS obtained in this study contribute to the
understanding of the rate of collisionless reconnection mediated by the
plasmoid instability.

###### Acknowledgements.

C.G. acknowledges the financial support provided by the Universitée Franco-
Italienne through the Vinci program. This work benefits from the support of
the CNR contract DFM.AD003.261 (IGNITOR)-Del. CIPE n.79 del 07/08/2017 and of
the grant DOE DESC0021254. The numerical simulations were performed using the
EUROfusion high performance computer Marconi Fusion and Gallileo100 hosted at
CINECA (project FUA35-FKMR and IsC86$\\_$MR-EFLRA), and the computing
facilities provided by Mesocentre SIGAMME hosted by Observatoire de la Côte
d’Azur.

## References

  * Zweibel and Yamada (2009) E. G. Zweibel and M. Yamada, Annu. Rev. Astron. Astrophys. 47, 291 (2009).
  * Yamada  _et al._ (2010) M. Yamada, R. Kulsrud, and H. Ji, Reviews of Modern Physics 82, 603 (2010).
  * Ji  _et al._ (2022) H. Ji, W. Daughton, J. Jara-Almonte, A. Le, A. Stanier, and J. Yoo, Nature Reviews Physics 4, 263 (2022), 2202.09004 .
  * Daughton  _et al._ (2009) W. Daughton, V. Roytershteyn, B. J. Albright, H. Karimabadi, L. Yin, and K. J. Bowers, Phys. Rev. Lett. 103, 065004 (2009).
  * Bhattacharjee  _et al._ (2009) A. Bhattacharjee, Y.-M. Huang, H. Yang, and B. Rogers, Physics of Plasmas 16, 112102 (2009), 0906.5599 .
  * Sweet (1958) P. A. Sweet, in  _Electromagnetic Phenomena in Cosmical Physics_ , Vol. 6, edited by B. Lehnert (1958) p. 123.
  * Parker (1957) E. N. Parker, J. Geophys. Res. 62, 509 (1957).
  * Biskamp (1986) D. Biskamp, The Physics of Fluids 29, 1520 (1986).
  * Huang and Bhattacharjee (2010) Y.-M. Huang and A. Bhattacharjee, Physics of Plasmas 17, 062104 (2010), 1003.5951 .
  * Uzdensky  _et al._ (2010) D. A. Uzdensky, N. F. Loureiro, and A. A. Schekochihin, Phys. Rev. Lett. 105, 235002 (2010).
  * Comisso  _et al._ (2015a) L. Comisso, D. Grasso, and F. L. Waelbroeck, Physics of Plasmas 22, 042109 (2015a), 1502.06408 .
  * Comisso and Grasso (2016) L. Comisso and D. Grasso, Physics of Plasmas 23, 032111 (2016), 1603.00090 .
  * Comisso and Bhattacharjee (2016) L. Comisso and A. Bhattacharjee, Journal of Plasma Physics 82, 595820601 (2016), 1609.02998 .
  * Ji and Daughton (2011) H. Ji and W. Daughton, Physics of Plasmas 18, 111207 (2011).
  * Daughton and Roytershteyn (2012) W. Daughton and V. Roytershteyn, Space Science Reviews 172, 271 (2012).
  * Huang and Bhattacharjee (2013) Y.-M. Huang and A. Bhattacharjee, Physics of Plasmas 20, 055702 (2013), 1301.0331 .
  * Karimabadi and Lazarian (2013) H. Karimabadi and A. Lazarian, Physics of Plasmas 20, 112102 (2013).
  * Comisso  _et al._ (2015b) L. Comisso, D. Grasso, and F. L. Waelbroeck, Journal of Plasma Physics 81, 495810510 (2015b), 1505.00726 .
  * Le  _et al._ (2015) A. Le, J. Egedal, W. Daughton, V. Roytershteyn, H. Karimabadi, and C. Forest, Journal of Plasma Physics 81, 305810108 (2015).
  * Loureiro and Uzdensky (2015) N. F. Loureiro and D. A. Uzdensky, Plasma Physics and Controlled Fusion 58, 014021 (2015).
  * Bhat and Loureiro (2018) P. Bhat and N. F. Loureiro, Journal of Plasma Physics 84, 905840607 (2018), 1804.05145 .
  * Burch  _et al._ (2016) J. L. Burch, R. B. Torbert, T. D. Phan, L.-J. Chen, T. E. Moore, R. E. Ergun, J. P. Eastwood, D. J. Gershman, P. A. Cassak, M. R. Argall, S. Wang, M. Hesse, C. J. Pollock, B. L. Giles, R. Nakamura, B. H. Mauk, S. A. Fuselier, C. T. Russell, R. J. Strangeway, J. F. Drake, M. A. Shay, Y. V. Khotyaintsev, P.-A. Lindqvist, G. Marklund, F. D. Wilder, D. T. Young, K. Torkar, J. Goldstein, J. C. Dorelli, L. A. Avanov, M. Oka, D. N. Baker, A. N. Jaynes, K. A. Goodrich, I. J. Cohen, D. L. Turner, J. F. Fennell, J. B. Blake, J. Clemmons, M. Goldman, D. Newman, S. M. Petrinec, K. J. Trattner, B. Lavraud, P. H. Reiff, W. Baumjohann, W. Magnes, M. Steller, W. Lewis, Y. Saito, V. Coffey, and M. Chandler, Science 352, aaf2939 (2016).
  * Phan  _et al._ (2018) T. Phan, J. Eastwood, M. Shay, J. Drake, B. Sonnerup, M. Fujimoto, P. Cassak, M. Oieroset, J. Burch, R. Torbert, A. Rager, J. Dorelli, D. Gershman, C. Pollock, P. Sharma Pyakurel, C. Haggerty, Y. Khotyaintsev, B. Lavraud, Y. Saito, and W. Magnes, Nature 557 (2018).
  * Olson  _et al._ (2016) J. Olson, J. Egedal, S. Greess, R. Myers, M. Clark, D. Endrizzi, K. Flanagan, J. Milhone, E. Peterson, J. Wallace, D. Weisberg, and C. Forest, Physical Review Letters 116 (2016).
  * Cafaro  _et al._ (1998) E. Cafaro, D. Grasso, F. Pegoraro, F. Porcelli, and A. Saluzzi, Phys. Rev. Lett. 80, 4430 (1998).
  * Ottaviani and Porcelli (1993) M. Ottaviani and F. Porcelli, Phys. Rev. Lett. 71, 3802 (1993).
  * Waelbroeck (1989) F. L. Waelbroeck, Physics of Fluids B 1, 2372 (1989).
  * Jemella  _et al._ (2003) B. D. Jemella, M. A. Shay, J. F. Drake, and B. N. Rogers, Phys. Rev. Lett. 91, 125002 (2003).
  * Jemella  _et al._ (2004) B. D. Jemella, J. F. Drake, and M. A. Shay, Physics of Plasmas 11, 5668 (2004).
  * Loureiro  _et al._ (2005) N. F. Loureiro, S. C. Cowley, W. D. Dorland, M. G. Haines, and A. A. Schekochihin, Phys. Rev. Lett. 95, 235003 (2005).
  * Del Sarto  _et al._ (2003) D. Del Sarto, F. Califano, and F. Pegoraro, Physical Review Letters 91, 235001 (2003).
  * Grasso  _et al._ (2001) D. Grasso, F. Califano, F. Pegoraro, and F. Porcelli, Phys. Rev. Lett. 86, 5051 (2001).
  * Cassak and Drake (2013) P. A. Cassak and J. F. Drake, Physics of Plasmas 20, 061207 (2013).
  * Comisso  _et al._ (2016) L. Comisso, M. Lingam, Y.-M. Huang, and A. Bhattacharjee, Physics of Plasmas 23, 100702 (2016).
  * Comisso  _et al._ (2017) L. Comisso, M. Lingam, Y.-M. Huang, and A. Bhattacharjee, The Astrophysical Journal 850, 142 (2017).
  * Huang  _et al._ (2017) Y.-M. Huang, L. Comisso, and A. Bhattacharjee, Astrophys. J. 849, 75 (2017), 1707.01863 .
  * Biskamp (2000) D. Biskamp,  _Magnetic Reconnection in Plasmas_, Cambridge Monographs on Plasma Physics (Cambridge University Press, 2000).
  * Porcelli (1991) F. Porcelli, Phys. Rev. Lett. 66, 425 (1991).


# Backfilling Cohorts in Phase I Dose-Escalation Studies

Helen Barnett1, Oliver Boix2, Dimitris Kontos3, Thomas Jaki1,4
1 MRC Biostatistics Unit, University of Cambridge
2 Bayer AG
3 ClinBAY
4 Department of Mathematics and Statistics, Lancaster University

###### Abstract

The use of ‘backfilling’, assigning additional patients to doses deemed safe,
in phase I dose-escalation studies has been used in practice to collect
additional information on the safety profile, pharmacokinetics and activity of
a drug. These additional patients help ensure that the Maximum Tolerated Dose
(MTD) is reliably estimated and give additional information in order to
determine the recommended phase II dose (RP2D). In this paper, we study the
effect of employing backfilling in a phase I trial on the estimation of the
MTD and the duration of the study. We consider the situation where only one
cycle of follow-up is used for escalation as well as the case where there may
be delayed onset toxicities. We find that, over a range of scenarios, there is
an increase in the proportion of correct selections and a notable reduction in
the trial duration at the cost of more patients required in the study.

Keywords:
Dose-Finding; Dose-Escalation; Backfilling; Phase I Trials; Model-Based; Late-
onset Toxicity.

##  1 Introduction

In Phase I dose-finding studies, the main objective is often to find the
Maximum Tolerated Dose (MTD) or the recommended Phase II dose (RP2D), the dose
recommended for further testing in Phase II. The MTD is defined as highest
dose that has an acceptable level of toxicity [11], most often corresponding
to a certain probability of occurrence of a Dose Limiting Toxicity event
(DLT). In oncology, a DLT is frequently defined as a grade 3 or higher
toxicity by the grading scale of the National Cancer Institute [7]. In the
following, we use the terms DLT and toxicity interchangeably.

In a Phase I dose-escalation study, a set of doses is investigated, patients
are recruited in cohorts and an escalation procedure is used to carefully
escalate from low doses that are expected to be very safe to dose levels that
have an acceptable level of toxicity and at the same time induce some
desirable activity in a patient. The escalation procedure can be rule-based
[e.g. 11], model-based [e.g. 13] or model-assisted [e.g. 5] and cohorts are
most commonly small [14], with a size of 3 often used in such trials.

However, with such small sample sizes, comes greater uncertainty in the
estimate of interest. It is, for example, desirable to establish the MTD
quickly and accurately but with no more patients than necessary. As with any
clinical trial, a balance must be taken between the accuracy, the duration of
the study and the trial size. Typically, the larger the trial, the higher the
accuracy, but also the longer the duration of the study.

An approach that has gained popularity in the recent years (e.g. [4, 10]), is
the use of ‘backfilling’ of cohorts on lower doses [2]. The principle is that,
once a dose is deemed safe enough to escalate to a higher dose-level,
additional patients may be allocated to lower doses to increase the
understanding of the safety, tolerability and activity of these doses. The
decision to backfill a dose may be taken solely on the criteria that the dose
is deemed ‘safe’, or it may require the additional condition that an activity
signal must be seen.

Whilst it is clear that backfilling will result in a, potentially large,
increase in patient numbers, additional insight around the safety,
tolerability and potential activity of the treatment is gained. Moreover,
there is the potential that the trial duration could be substantially reduced
due to the improved understanding of the dose-toxicity relationship. The exact
nature of this relationship between trial duration and sample size is,
however, unclear - in particular in the setting where late-onset toxicities
are of concern.

In this work we investigate the impact of backfilling on the probability of
correctly selecting the MTD and on the duration in phase I dose-escalation
trials. The rest of the paper is organized as follows. In Section 2, we give
an example of a dose-escalation trial that utilized backfilling. In Section 3,
we first introduce an example trial simulation with and without backfilling,
then describe the dose-finding algorithms used in this work. In Section 4, we
present the results from simulation studies to demonstrate the impact of
backfilling on operating characteristics. Finally, we conclude with a
discussion in Section 5.

##  2 Motivating Trial Example

The first in-human Phase I study reported in [12] investigated the safety and
activity of the activin A inhibitor, STM 434, in advanced solid tumors. The
study initally considered five doses levels, 0.25, 0.5, 1, 2, and 4 mg/kg
administered every 4 weeks. The treatment scheduled subsequently changed to
bi-weekly after the half-life was estimated to be lower than anticipated and
an additional dose of 8 mg/kg was added later due to lower than anticipated
predicted exposures. A 3+3 design [11] was used to guide dose-escalation.
Dose-limiting toxicities (DLT) were deﬁned as any grade $\geq$3 nonhematologic
toxicity, any grade $\geq$4 hematologic toxicity lasting 7 days, febrile
neutropenia, or grade 3 thrombocytopenia with active bleeding and the DLT
assessement period was 28 days. A minimum of 3 evaluable patients were
required prior to dose-escalation and ”Backfill” slots were permitted at doses
that had been declared safe.

A total of 32 patients participated in the trial of which three experienced a
DLT and of 28 patients that were evaluable according to RECIST [3], 16
achieved stable disease while the remaining had progressive disease (Table 1).
Since two of the DLTs occurred in the highest dose, the MTD was declared to be
4 mg/kg every 2 weeks but no dose-expansion was undertaken following the
safety review committees recommendation on the basis of the overall safety
profile observed.

Dose (mg/kg) | 0.25 | 0.5 | 0.5 | 0.75 | 1 | 2 | 4 | 8
---|---|---|---|---|---|---|---|---
| 4 weekly | 4 weekly | 2 weekly | 2 weekly | 2 weekly | 2 weekly | 2 weekly | 2 weekly
DLT | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 2
SD | 1 | 3 | 3 | 1 | 2 | 3 | 2 | 1
Total subjects | 4 | 6 | 4 | 3 | 4 | 4 | 4 | 3
Table 1: Number of subjects and DLT per dose-level.

From the results of the study we can see that backfilling slots were indeed
used as several of the cohorts have more than three (or six) patients required
for the 3+3 design.

Note, that the 3+3 design does, however, not utilize these additional data
when recommending subsequent doses. In the remainder of this work we will
therefore explore the impact of backfilling on the operating characteristics
of a trial that utilizes a model-based escalation approach. This exploration
is motiviated by a recent study of Thorium-227 in combination with an antibody
(NCT03507452). As this study utilized doses of Thorium-227 of 1.5 MBq in steps
of 1.0 or 1.5 MBq, with antibody doses of 10 mg, we will consider the same
dose-levels in our subsequent evaluations. Moreover, due to the radioactive
nature of Thorium-227, late-onset toxicities $-$ toxicities that occur after
the first treatment cycle $-$ are of potential concern. We therefore also
explore the impact of backfilling under a traditional CRM-type model as well
as a time-to-event model.

##  3 Methodology

In this section we outline the methods used in this work. In Section 3.1 we
illustrate and discuss the concept of backfilling in more detail. We then
describe the statistical models used for dose-escalation in Section 3.2, and
finally outline the rules used in the design implementation in Section 3.3.

###  3.1 Single Simulation Examples

To illustrate the concept of backfilling, and the difference between trials
that use backfilling and those that do not, we present an example of a single
trial simulation. Following the motivating study (NCT03507452) the following
six doses are investigated: 1.5MBq , 2.5MBq , 3.5MBq , 4.5MBq , 6.0MBq ,
7.0MBq. Cohorts of size three enter the trial, starting at the lowest dose. If
no backfilling is implemented, then cohorts are assigned to doses according to
the escalation procedure, until some stopping rule is triggered. If full
backfilling is implemented then, when a dose is considered safe enough that
the escalation continues above it, two additional backfilling cohorts are
recruited to that dose. Of course one may also choose to backfill more
cautiously, but in this illustration we use two cohorts as standard. Full
backfilling therefore means that we start backfilling from the lowest dose of
1.5MBq. Patients are followed up for three cycles of treatment, with new
cohorts assigned after each cycle, with a cycle lasting 6 weeks. Since the
follow-up period is longer than the times between assignments, the Time-to-
Event Continual Reassessment (TITE-CRM) is implemented for the dose-escalation
procedure as described in Section 3.2.

Figure 1 illustrates a trial is conducted without backfilling (left panel) and
with full backfilling (right panel). Every other aspect of the trials was
comparable (see Section 3.3 for further details). In this case, the trial
without backfilling determined the MTD to be the 4.5MBq dose while the trial
using backfilling recommended 3.5MBq. Figure 1 shows that the escalation in
both settings is the same - they both reach the 4.5MBq dose, but in the trial
without backfilling two further cohorts are allocated to this dose before
being recommended as the MTD. When backfilled cohorts were used, more
information was available at the lower doses allowing the precision stopping
rule to be reached quickly and hence the lower dose is recommended. In
addition to the difference in recommended dose one can also observe that the
duration of the study is notably shorter when backfilling is used.

Figure 1: A single trial simulation using the TITE-CRM, with and without
backfilling. Note that the backfilled cohorts in grey all enter the trial at
the same time as the escalation cohort above in black, the horizontal offset
is for ease of reading.

###  3.2 BLRM and TITE-BLRM

In order to investigate the impact of backfilling on the estimation of the MTD
and the duration of the trial, we consider two settings. The first setting
assumes that patients are followed up for one cycle of treatment only (i.e.
DLT period 1 cycle), and the next cohort of patients is assigned once the
previous cohort’s follow-up period has been fully observed. The second setting
assumes that there may be late onset toxicities. Patients are therefore
followed up for three cycles of treatment. A new cohort of patients is
admitted every cycle, so that only partial information is available for the
previous two cohorts, as their full follow-up period has not yet been
observed. In the first setting we use a Baysian Logistic regression model
(BLRM) [8], and in the second setting a Time-To-Event version of the BLRM
(TITE-CRM) [1].

In each setting, a set of $J$ doses labelled $d_{j}$ for $j=1,\ldots,J$ are
investigated, with patients labelled $i=1,\ldots,n$.

The BLRM is conducted in the following way. Starting at the lowest dose
$d_{1}$, cohorts of patients enter the trial. After each cohort has been fully
observed, the dose assignment of the next cohort is decided. A two-parameter
logistic model is used to describe the dose-response relationship:

$F(d,\bm{\beta})=\frac{\exp(\beta_{0}+\beta_{1}d)}{1+\exp(\beta_{0}+\beta_{1}d)},$

where $F(d,\bm{\beta})$ is the probability of DLT at dose $d$, and
$\bm{\beta}=(\beta_{0},\beta_{1})$ is parameter vector with prior:

$\begin{pmatrix}\beta_{0}\\\ \log(\beta_{1})\end{pmatrix}\sim
N_{2}\left(\begin{pmatrix}c_{1}\\\
c_{2}\end{pmatrix},\begin{pmatrix}v_{1}&0\\\ 0&v_{2}\end{pmatrix}\right).$

The posterior distribution of $\bm{\beta}$ is updated after each cycle using
the likelihood

${\mathcal{L}(\bm{\beta})=\prod_{i=1}^{n}F(d,\bm{\beta})^{y_{i}}\\{1-F(d,\bm{\beta})\\}^{1-y_{i}}},$

where $y_{i}$ is an indicator taking the value 1 if patient $i$ had observed a
DLT response and 0 otherwise. The updated posterior for $\bm{\beta}$ is then
used to estimate the probability of DLT at each dose. The dose assignment of
the next cohort is then dose $d_{j}$ that minimizes
$|F(d_{j},\hat{\bm{\beta}})-\tau_{1}|$, where $\tau_{1}$ is the target DLT
rate for one cycle of follow-up, subject to certain rules - see Section 3.3.
The final dose recommendation is then the dose $d_{j}$ that
minimizes$|F(d_{j},\bm{\beta})-\tau_{1}|$ once a stopping rule has been
triggered.

In the second setting, the TITE-CRM [1] is used. Here, patients are followed
up for three cycles of treatment. However, if it was required to wait until
the entire follow-up for the previous cohort had been completely observed to
assign the dose for the next patient, the trial length would be very
undesirably long. Therefore, each new cohort is assigned their dose once the
previous cohort has been observed for one cycle of treatment. There is
therefore only partial information available for the previous two cohorts. The
TITE-CRM takes this into account by weighting the observations in the
following way.

The dose response model $F(d,\bm{\beta})$ is identical to (3.2), and is
weighted to form $G(d,w,\bm{\beta})$:

${G(d,w,\bm{\beta})=wF(d,\bm{\beta)}},$

where the weights $0\leq w\leq 1$ are a function of time-to-event of a patient
response. The posterior distribution of $\bm{\beta}$ is updated after each
cycle using likelihood

${\mathcal{L}(\bm{\beta})=\prod_{i=1}^{n}G(d_{[i]},w_{i,n},\bm{\beta})^{y_{i,n}}\\{1-G(d_{[i]},w_{i,n},\bm{\beta})\\}^{1-y_{i,n}}},$

where $n$ is the number of patients that have been treated so far, $d_{[i]}$
is the dose assigned to patient $i$ and $y_{i,n}$ is an indicator which takes
the value 1 if patient $i$ has observed a DLT after the $n$ patients treated
so far have been observed for at least one cycle, and 0 otherwise.

The updated posterior for $\bm{\beta}$ is used to estimate the probability of
DLT at each dose, as in the one cycle setting. The dose assignment of the next
cohort is then the dose $d_{j}$ that minimizes
$|F(d_{j},\bm{\beta})-\tau_{3}|$, where $\tau_{3}$ is the target DLT rate for
three cycles of follow-up, again subject to certain rules - see Section 3.3.

In this implementation, we use the simple specification of weights suggested
by Cheung and Chappell [1]: $w_{i,n}=u_{i,n}/S$, where $u_{i,n}$ is the
current number of cycles for which patient $i$ has been observed and $S$ is
the total number of cycles in the follow-up period. If patient $i$ observes a
DLT response, then $w_{i,n}=1$. The final dose recommendation is the dose
level that minimises $|F(d_{j},\hat{\beta})-\tau_{3}|$ once a stopping rule
has been implemented and the follow-up for all enrolled patients has been
completed.

Due to the added complexity of the method, the TITE-CRM uses a start-up period
such that the dose assignment is escalated one level at a time until a DLT
response is observed. Once a DLT response is observed, then the TITE-CRM model
is used. Although in the original description of the methodology, this start-
up period requires each patient to be followed up for their entire follow-up
time before the next patient’s dose is assigned, in our implementation, only
one cycle is required for follow-up before the next is assigned, which is in
line with the rest of the trial.

###  3.3 Rules

To evaluate the performance in the setting of a Phase I dose-escalation trial,
we use enforcement and stopping rules that could be used in such a trial. We
define $p_{s,d_{j}}$ as the $P(DLT)$ in cycles up to and including cycle $s$
for any given dose $d_{j}$.

Enforcement Rules:

  1. 1.

Hard Safety: If there is a high probability that the toxicity of an
experimented dose exceeds the target toxicity, this and all higher doses are
excluded from further experimentation (i.e. dose $d_{j}$ and all above are
excluded when $P(p_{1,d_{j}}>\tau_{1})>\psi$ for some threshold $\psi$). In
this implementation we use a threshold for excessive toxicity of $\psi=0.95$,
with a $Beta(1,1)$ prior for the DLT outcome. For example for $\tau_{1}=0.3$,
if there are at least 3 DLT responses out of 3 patients, or at least 4 DLT
responses out of 6 patients, or at least 5 DLT responses out of 9 patients,
then all dose assignments must be lower than that dose for the rest of the
study. If the lowest dose is excluded then the trial stops with no dose
recommendation made.

  2. 2.

K-fold Skipping Doses: No more than a 2-fold-rise in dose value for the next
dose assignment based on the highest experimented dose so far.

Stopping Rules:

  1. 1.

Sufficient Information: If a dose is recommended for the next cohort on which
three cohorts have already been assigned in the escalation (excluding
backfilling cohorts) the trial is stopped.

  2. 2.

Lowest Dose Deemed Unsafe: If $P(p_{1,d_{1}}>30\%)>0.80$ according to the
escalation model and at least one cohort of patients has been assigned to dose
$d_{1}$, the trial is stopped.

  3. 3.

Highest Dose Deemed Very Safe: If $P(p_{1,d_{J}}\leq 30\%)>0.80$ according to
the model and at least one cohort of patients has been assigned to dose
$d_{J}$, the trial is stopped.

  4. 4.

Precision: If the MTD is estimated precisely enough, the trial is stopped.
This precision is defined as $CV(MTD)<30\%$, with the coefficient of variation
calculated as an adjusted median absolute deviation divided by the median.
This stopping rule is only used once at least three cohorts of patients that
are not part of a backfilling cohort have had at least one cycle of treatment
in the escalation. The data of the backfilling cohorts are included in the
estimation of the precision, however.

  5. 5.

Hard Safety: If the lowest dose is considered unsafe according to the hard
safety enforcement rule, the trial is stopped.

  6. 6.

Maximum Patients: If the maximum number of patients ($n=n_{\mbox{max}}$) have
been recruited.

##  4 Simulations

###  4.1 Set-up

In order to investigate the effect of backfilling on the operating
characteristics of the dose-finding designs, we conduct a simulation study.
The same set of six doses as used in the example in Section 3.1 are used:
1.5MBq , 2.5MBq , 3.5MBq , 4.5MBq , 6.0MBq , 7.0MBq. Seventeen scenarios are
considered, to cover a wide range of potential dose responses. Table 2 gives
the probability of a DLT in the first six week cycle for each of the
considered scenarios, with the MTD associated with the target, $\tau_{1}=0.3$,
highlighted in boldface. Scenarios 1 – 6 represent cases where each level in
turn is the MTD, with higher and lower doses equidistant in terms of
probability of DLT. Scenarios 7 and 8 are non-linear around the MTD. In
scenarios 9 – 12, no dose is exactly on target, with all doses unsafe in
scenario 9. Scenarios 13 – 17 are a set of varying scenarios, typically used
to test performance of a dose-finding algorithm.

Scenario | 1.5MBq | 2.5MBq | 3.5MBq | 4.5MBq | 6.0MBq | 7.0MBq
---|---|---|---|---|---|---
1 | 0.30 | 0.40 | 0.50 | 0.60 | 0.70 | 0.80
2 | 0.20 | 0.30 | 0.40 | 0.50 | 0.60 | 0.70
3 | 0.10 | 0.20 | 0.30 | 0.40 | 0.50 | 0.60
4 | 0.05 | 0.10 | 0.20 | 0.30 | 0.40 | 0.50
5 | 0.05 | 0.10 | 0.15 | 0.20 | 0.30 | 0.40
6 | 0.02 | 0.05 | 0.10 | 0.15 | 0.20 | 0.30
7 | 0.15 | 0.20 | 0.25 | 0.30 | 0.45 | 0.60
8 | 0.05 | 0.15 | 0.30 | 0.35 | 0.40 | 0.45
9 | 0.40 | 0.45 | 0.50 | 0.55 | 0.60 | 0.65
10 | 0.05 | 0.15 | 0.25 | 0.35 | 0.45 | 0.55
11 | 0.15 | 0.20 | 0.35 | 0.40 | 0.45 | 0.50
12 | 0.05 | 0.10 | 0.15 | 0.20 | 0.25 | 0.40
13 | 0.06 | 0.07 | 0.08 | 0.09 | 0.11 | 0.12
14 | 0.10 | 0.14 | 0.21 | 0.30 | 0.46 | 0.58
15 | 0.16 | 0.30 | 0.50 | 0.70 | 0.89 | 0.95
16 | 0.55 | 0.91 | 0.99 | 1.00 | 1.00 | 1.00
17 | 0.05 | 0.05 | 0.05 | 0.80 | 0.80 | 0.80
Table 2: $p_{1,d_{j}}=P(DLT)$ in cycle 1 for each dose $d_{j}$ across 17
defined scenarios. MTD in bold.

For the setting where the follow-up period is three cycles, the conditional
probability of DLT response in subsequent cycles is multiplied by a factor of
1/3. So that

$p_{2,d_{j}}=p_{1,d_{j}}+(1-p_{1,d_{j}})\frac{p_{1,d_{j}}}{3},$

and

$p_{3,d_{j}}=p_{1,d_{j}}+(1-p_{1,d_{j}})\frac{p_{1,d_{j}}}{3}+(1-p_{1,d_{j}})(1-\frac{p_{1,d_{j}}}{3})\frac{p_{1,d_{j}}}{9}.$

For example if $p_{1,d_{j}}=0.3$, then $p_{2,d_{j}}=0.37$ and
$p_{3,d_{j}}=0.391$, hence we use $\tau_{3}=0.391$ as the target toxicity for
three cycles of follow-up.

The maximum sample size is chosen to be $n_{\mbox{max}}=54$ which is
relatively large for a Phase I trial, but has been chosen to allow escalation
to dose $d_{6}$ with backfilling implemented.

We conduct 5,000 simulations for each scenario, and compare the performance to
the non-parametric benchmark [9]. This benchmark gives an indication of the
‘difficulty’ of a scenario, so that we can quantify any differences between
the approaches accordingly. Note that this uses the maximum sample size in
every simulation, and is not subject to any stopping rules.

In these implementations, as well as comparing backfilling all doses
considered safe (Fully Backfilled), and not (Not Backfilled), we also consider
the setting where a dose level is only backfilled once an activity signal has
been observed at that dose or any dose below (Partially Backfilled). For
example, if an activity signal is first seen in the third dose level and not
the first or second, then the decision to escalate to the fourth dose level
would mean the additional backfilling cohorts are assigned to dose three but
not doses one and two. In this partially backfilling setting, the underlying
probability of observing an activity signal (at least one complete response in
a cohort) is 0.00, 0.15, 0.30, 0.45, 0.60, 0.75 at each of the six doses
respectively.

###  4.2 Prior Specification

We consider two options for the values of the hyper-parameters of the prior
distribution. In the first option, the hyper-parameters for the prior used
are: $c_{1}=\log(0.5)$, $c_{2}=0$, $v_{1}=4$ and $v_{2}=1$, chosen so that the
prior is relatively vague, with a mean effective prior sample size of 1.3
patients per dose level. These choices are in line with those used by
Neuenschwander et al. [8]. Importantly, these hyper-parameters are the same
for the implementation with and without backfilling, and for the follow-up of
one and three cycles. The decision to choose the same prior for both one cycle
and three cycles is to ensure a fair comparison between the two approaches. It
is important, however, to note that in practice, different prior distributions
are likely suitable for the different settings, as explored further below.

An alternative option considered here is to calibrate the values of these
hyper-parameters over a small range of scenarios, to choose the values
yielding the best performance across a wide range of settings (see for example
Mozgunov et al. [6]). We have chosen scenarios 1, 3, 4, 6, 9 and 13 to
represent a diverse set of dose-response relationships. The calibration is
done separately for the four settings resulting in the combination of 1
cycle/3 cycles and backfilling/no backfilling, in order to give each approach
the maximum chance of successful performance. It is important to note that the
prior used for both fully and partially backfilling is the same in each case,
and has been calibrated under the assumption of full backfilling. The values
of these hyper-parameters resulting from the calibration procedure are
provided in Table 3.

We have chosen to implement both approaches to the prior specification here as
we wish to demonstrate a direct comparison between backfilling and not
backfilling using the same prior for both, as well as demonstrating the
potential of both approaches with a suitably matched prior. In practice it may
be the case that the decision to backfill is made during the design stage, in
which case we would use the corresponding calibrated prior. However, it may
instead be the case that backfilling is introduced later in the trial (per
amendment), and the same specification of prior that was intended to be used
without backfilling must still be used. Hence we highlight the importance of
investigating the performance for both a calibrated and non-calibrated prior.

| 1 cycle | 3 cycles
---|---|---
No Backfilling |  $c_{1}=\log(1/2)$, $c_{2}=\log(1/4)$, |  $c_{1}=\log(1/16)$, $c_{2}=\log(1/16)$,
|  $v_{1}=4$ and $v_{2}=2$ |  $v_{1}=1$ and $v_{2}=1$
Backfilling |  $c_{1}=\log(1/2)$, $c_{2}=\log(1/4)$, |  $c_{1}=\log(1/16)$, $c_{2}=\log(1/16)$,
(Partial & Full) |  $v_{1}=2$ and $v_{2}=1$ |  $v_{1}=1$ and $v_{2}=1$
Table 3: Calibrated values for prior hyper-parameters.

###  4.3 Results

The most often considered metrics of performance in dose-finding trials are
the proportion of correct selections (PCS) and the proportion of acceptable
selections (PAS). A correct selection is defined as selecting the MTD as
defined in Table 2, or as a safety stopping rule being correctly triggered. An
acceptable selection is defined as a dose whose true probability of toxicity
during the first cycle is between 0.18 and 0.33. Why this range? Any reference
that could be given to why?

####  4.3.1 Non-Calibrated Prior

Figure 2: Proportion of correct and acceptable selections across scenarios for
one cycle of follow-up Figure 3: Proportion of correct and acceptable
selections across scenarios for three cycles of follow-up

Figures 2 and 3 show the PCS and PAS in the setting with one cycle and three
cycles respectively. From Figure 2, it can be seen that, when one cycle of
follow-up is considered, employing backfilling increases the PCS and PAS in
almost all scenarios, as expected. The only scenario which did not see an
increase in PCS is scenario 7, where there is a marginal decrease in PCS, but
a noticeable increase in PAS. This is a particularly challenging scenario,
where the dose below the MTD has true probability of DLT only 5% below target.
The largest increase is seen in scenario 13, where full backfilling gives a 9%
increase in correct selections. In nearly all scenarios, partial backfilling
gives a performance between no backfilling and fully backfilling as expected.
Interestingly in scenario 17, the performance of partial backfilling is worse
than both no and full backfilling. With such an ‘easy’ scenario, the PCS is
already very high, and the decrease is due to a shift to recommending the
unsafe fourth dose. The average increase in correct selections from no
backfilling to fully backfilling is 4%, and is 5% for acceptable selections.
The approach to partially backfilling we have taken means the probability of
backfilling any given level dose increases with dose. Therefore when the lower
doses are more toxic, fewer doses will be backfilled and partial backfilling
gives more similar results to no backfilling. Although in all cases it is
noticeable that the performance is well below that of the benchmark, it is
worth noting that the benchmark has the advantage of a much larger sample size
on average.

To compare the measures of trial sizes of the two approaches, Figure 4 shows
the relationship between the mean total sample size and the mean trial
duration for each scenario. The diamonds represent the setting with one cycle
of follow-up, with blue indicating that full backfilling was used, purple
indicating partial backfilling and black indicating no backfilling. It is
clear to see that the use of full backfilling substantially increases the
total sample size in scenarios where the MTD is at least the second dose
level. It is notable, however, that this increase in sample size owed to the
use of backfilling decreases the average trial duration notably. On average
across scenarios, full backfilling increases sample size by 12 patients on
average but reduces the trial duration by 6 weeks. It appears that on average
each additional patient reduces the duration by half a week. For partial
backfilling, although the overall increase in sample size is smaller, as is
the decrease in trial duration, the trade-off is still on average each
additional patient reduces the duration by half a week.

The other important metrics used in dose-finding trials concern the safety of
the patients within the trial. Figure 5 shows the distribution of mean
percentage of DLT responses and the mean number of patients assigned to overly
toxic doses. Although backfilling increases the number of DLT responses
observed in every scenario, the percentage of patients observing a DLT
response decreases slightly when full backfilling is used, and is largely
unchanged for partial backfilling. For the most part, the increase in numbers
of patients exposed to overly toxic doses are in scenarios where all doses are
unsafe. Hence one mistaken escalation increases this exposure by three cohorts
when backfilling instead of one when not. The benefit in reducing the number
of patients exposed to overly toxic doses comes in scenarios where the MTD is
in the middle of the dose range, and the backfilling expansion provides more
information on the lower doses, resulting in more cautious escalation.

Figure 4: Mean total sample size and mean trial duration across scenarios.
Darker colours indicate scenarios where a higher number of doses are unsafe.
Figure 5: Mean number of DLTs and mean number of patients treated at overly
toxic doses across scenarios. Darker colours indicate scenarios where a higher
number of doses are unsafe.

Interestingly, the patterns observed in the setting with one cycle of follow-
up are not all duplicated in the setting with three cycles of follow-up. It is
not the case here that the PCS and PAS increase across all scenarios when
backfilling is employed. In fact in some scenarios, there is a noticeable
decrease in PCS when backfilling is employed. For example, in scenarios 5 – 7,
where the PCS is low in both settings, it is lower when backfilling is used.
This reflects further the point raised earlier that the backfilling leads to a
more cautious escalation. In most cases, the PAS is higher when backfilling is
implemented.

In terms of mean sample size and trial duration, the circles on Figure 4
display this relationship. This relationship is similar to that observed with
one cycle. For every additional patient, the trial duration decreases by 0.3
weeks, for both full and partial backfilling. Again, the magnitude of this
overall increase or decrease is larger for fully backfilling. Likewise, the
comparisons of the safety aspects of the trials are similar. With three cycles
of observation, fully backfilling increases the total number of DLTs in all
scenarios, decreases the percentage of DLT responses overall, and decreases
the number of patients exposed to overly toxic doses in some scenarios, and
increases in others. On average across scenarios, there is a slight decrease
in the number of patients exposed to overly toxic doses.

####  4.3.2 Calibrated Prior

Now that we have investigated the behaviour of each approach when the same
prior specification is used, we move to look at the potential of each approach
when the prior is calibrated for the design. Note that for the case of partial
backfilling, the calibrated prior with full backfilling is used.

Figures 6 and 7 illustrate the proportion of correct and acceptable selections
in the settings with a toxicity assessment period of one and three cycles,
respectively. All but scenarios 6 and 7 see an increase in the proportion of
correct selections when full backfilling is employed in the setting using one
cycle. In most scenarios the partial backfilling gives a level of performance
in between the full backfilling and no backfilling as before. However,
interestingly in some scenarios, such as scenarios 3, 8 and 10, partially
backfilling gives a higher proportion of correct and acceptable selections
than both when using full and no backfilling. In these scenarios, the MTD is
the third dose level. When the MTD is the fourth dose level, partially
backfilling performs worse than the other two approaches. This can be
explained by the calibration which assumes full backfilling. But since full
and partial backfilling result in different levels of information at lower
doses, the calibration is no longer optimal for partial backfilling. From this
we can deduce that calibrating the prior for the setting at hand has distinct
advantages.

When the DLT period is three cycles, a larger difference between backfilling
and not backfilling can be seen. In most cases, full backfilling increases the
proportion of correct selections (up to 21% in scenario 10). In some cases,
however, backfilling again does not increase the proportion of correct
selections (e.g. scenarios 5 and 6 where the MTD is in the higher dose range).
Partial backfilling once more gives a performance level in between no and full
backfilling as previously.

Figure 8 shows the percentage of patients who experience DLT responses, and
the number of patients treated at overly toxic doses. When backfilling is
employed, more patients are on average treated at overly toxic doses. The
exceptions are scenarios 5 and 12, where the fifth highest dose is the MTD and
therefore there is only one overly toxic dose. In terms of percentage of
patients observing DLTs, this is less than or equal in every scenario for the
use of backfilling. The percentage of patients observing a DLT response is
less when backfilling is employed.

The mean trial duration is drastically reduced with the use of full
backfilling, on average by seven weeks. The largest reduction of 10 weeks is
seen in scenarios 4 and 10, when the MTD is in the middle of the dose range,
for both a DLT period of one and three cycles. These are also scenarios where
the mean total sample size increases by a large amount. The mean total sample
size increases on average 16 patients when one cycle of follow-up is used, and
18 patients when three is used. For every additional patient, over the 17
considered scenarios, the trial duration reduces by 0.4 weeks for both one and
three cycles of follow-up.

Interestingly, the use of partial backfilling on average slightly increases
the trial duration over the setting with no backfilling when one cycle of
follow-up is used, driven by an increase in duration in scenarios where the
MTD is mid dose-range. This is once again a consequence of the prior that is
calibrated for full backfilling.

Figure 6: Proportion of correct and acceptable selections across scenarios for
one cycle of follow-up, using calibrated priors. Figure 7: Proportion of
correct and acceptable selections across scenarios for three cycles of follow-
up, using calibrated priors. Figure 8: Mean number of DLTs and mean number of
patients treated at overly toxic doses across scenarios. Darker colours
indicate scenarios where a higher number of doses are unsafe, using calibrated
priors. Figure 9: Mean total sample size and mean trial duration across
scenarios. Darker colours indicate scenarios where a higher number of doses
are unsafe, using calibrated priors.

##  5 Discussion

In this work, we have investigated the effect of backfilling on the operating
characteristics of dose-escalation studies. The main reasons to utilize
backfilling include to gain better understanding of the safety, tolerability
and activity of the treatment under investigation and to aid determining the
recommended Phase II dose. In this work we have focused on the implications of
backfilling on the estimation of the MTD and the duration of the study. We
found that backfilling increases the chance of identifying the MTD while
reducing the duration of the study. This comes at the cost of an increased
number of patients required in studies that use backfilling. The impact of
backfilling on the accuracy is larger in the setting with a DLT assessment
period of one cycle then when three cycles of follow-up are used.

In our evaluations, patients for backfilling are available for recruitment
immediately which clearly is an optimistic assumption. Moreover the use of
backfilling does increase the number of patients in the study. Nevertheless we
did see consistent benefits of backfilling. Specifically, one additional
patient did yield a reduction of the study duration of approximately half a
week.

Our investigations explore two different settings with respect to the prior
distributions used. In the first setting the prior is chosen with a study in
mind that does not plan to use backfilling, while the second considers
backfilling as an option at the outset. Encouragingly we find that,
irrespective of the setting, the benefits of backfilling on estimation of the
MTD and duration of the study are fairly consistent.

## Data Availability Statement

All data is simulated according to the specifications described.

## Acknowledgements

This report is supported by the NIHR Cambridge Biomedical Research Centre
(BRC-1215-20014). The views expressed in this publication are those of the
authors and not necessarily those of the NHS, the National Institute for
Health Research or the Department of Health and Social Care (DHSC). T Jaki and
H Barnett received funding from UK Medical Research Council (MC_UU_00002/14).

## References

  * [1] Ying Kuen Cheung and R. Chappell.  Sequential designs for phase I clinical trials with late-onset toxicities.  Biometrics, 56(4):1177–1182, 2000.
  * [2] Hakim-Moulay Dehbi, John O’Quigley, and Alexia Iasonos.  Controlled backfill in oncology dose-finding trials.  Contemporary clinical trials, page 106605, 2021.
  * [3] Elizabeth A Eisenhauer, Patrick Therasse, Jan Bogaerts, Lawrence H Schwartz, Danielle Sargent, Robert Ford, Janet Dancey, S Arbuck, Steve Gwyther, Margaret Mooney, et al.  New response evaluation criteria in solid tumours: revised recist guideline (version 1.1).  European journal of cancer, 45(2):228–247, 2009.
  * [4] Erika Paige Hamilton, Elizabeth Claire Dees, Judy Sing-Zan Wang, Amy Kim, Manav Korpal, Vicki Rimkunas, Nathalie Rioux, Joanne Schindler, and Dejan Juric.  Phase I dose escalation of H3B-6545, a first-in-class highly Selective ER$\alpha$ Covalent Antagonist (SERCA), in women with ER-positive, HER2-negative breast cancer (HR+ BC).  Journal of Clinical Oncology, 37(15_suppl):1059–1059, 2019.
  * [5] Suyu Liu and Ying Yuan.  Bayesian optimal interval designs for phase i clinical trials.  Journal of the Royal Statistical Society: Series C: Applied Statistics, pages 507–523, 2015.
  * [6] Pavel Mozgunov, Rochelle Knight, Helen Barnett, and Thomas Jaki.  Using an interaction parameter in model-based phase i trials for combination treatments? A simulation study.  International Journal of Environmental Research and Public Health, 18(1):1–19, 2021.
  * [7] National Cancer Instute.  Common Terminology Criteria for Adverse Events (CTCAE) v5.0, 2017.
  * [8] Beat Neuenschwander, Michael Branson, and Thomas Gsponer.  Critical aspects of the Bayesian approach to phase I cancer trials.  Statistics in Medicine, 27:2420–2439, 2008.
  * [9] John O’Quigley, Xavier Paoletti, and Jean Maccario.  Non-parametric optimal design in dose finding studies.  Biostatistics, 3(1):51–56, 2002.
  * [10] James Michael Pauff, Kyriakos P. Papadopoulos, Filip Janku, Anita Ahmed Turk, Lipika Goyal, Rachna T. Shroff, Toshio Shimizu, Masafumi Ikeda, Nilofer Saba Azad, James M. Cleary, Mary Linton Bounetheau Peters, Mitesh J. Borad, Kurt A. Jaeckle, Sani Haider Kizilbash, Rebecca Tupper, Carrie E. Furin, Matthew P. Hanley, Elizabeth Goodwin Hill, Xiaojian Xu, and James J. Harding.  A phase I study of LY3410738, a first-in-class covalent inhibitor of mutant IDH1 in cholangiocarcinoma and other advanced solid tumors.  Journal of Clinical Oncology, 39(3_suppl):TPS350–TPS350, 2021\.
  * [11] Barry E. Storer.  Phase I trials.  In C Redmond and T Colton, editors, Biostatistics in Clinical Trials, pages 337–342. John Wiley & Sons, Ltd, Chichester, 2001.
  * [12] Jessica J Tao, Nicholas A Cangemi, Vicky Makker, Karen A Cadoo, Joyce F Liu, Drew W Rasco, Willis H Navarro, Christopher M Haqq, and David M Hyman.  First-in-human phase i study of the activin a inhibitor, stm 434, in patients with granulosa cell ovarian cancer and other advanced solid tumors.  Clinical Cancer Research, 25(18):5458–5465, 2019.
  * [13] Graham M Wheeler, Adrian P Mander, Alun Bedding, Kristian Brock, Victoria Cornelius, Andrew P Grieve, Thomas Jaki, Sharon B Love, Christopher J Weir, Christina Yap, et al.  How to design a dose-finding study using the continual reassessment method.  BMC Medical Research Methodology, 19(1):1–15, 2019.
  * [14] Yinghui Zhou and John Whitehead.  Practical implementation of bayesian dose-escalation procedures.  Drug information journal: DIJ/Drug Information Association, 37(1):45–59, 2003.