S M Saberi Fathi
email: saberifathi@um.ac.ir

Jack A Tuszynski
email: jackt@ualberta.ca

Geometrical comparison of two protein structures using Wigner-D functions

Keywords: Protein structure, protein classification, computational methods, Wigner-D functions, root-mean-square deviation, tubulin, Fourier transform, reciprocal space, similarity value, correlation value

In this paper, we develop a quantitative comparison method for two arbitrary protein structures. This method uses a root-mean-square deviation (RMSD) characterization and employs a series expansion of the protein's shape function in terms of the Wigner-D functions to define a new criterion, which is called a "similarity value". We further demonstrate that the expansion coefficients for the shape function obtained with the help of the Wigner-D functions correspond to structure factors. Our method addresses the common problem of comparing two proteins with different numbers of atoms. We illustrate it with a worked example.

Introduction

A quantitative comparison of two protein tertiary structures to assess their similarity is a major challenge, but if properly investigated, it can offer answers to important questions in biochemistry and cell biology [START_REF] Li | Characterization of local geometry of protein surfaces with the visibility criterion[END_REF] . In particular, structural similarity between proteins is a very good predictor of their functional similarity. In order to classify proteins according to their structural characteristics, we first have to be able to determine the 3D structures of the proteins in question, which typically involves x-ray or electron crystallography, or in some cases other techniques such as nuclear magnetic resonance (NMR) or mass spectroscopy [START_REF] Rupp | Predictive models for protein crystallization[END_REF] . In the absence of crystallographic structures for a given protein, computational methods may still be used to predict a 3D structure based on sequence similarity with crystallographically resolved protein structures using a technique called homology modeling [START_REF] Arnold | The SWISS-MODEL workspace: a web-based environment for protein structure homology modelling[END_REF] . Assuming structural information is available, a number of methods have been developed to compare protein structures [START_REF] Kolodny | Protein structure comparison: implications for the nature of "fold space", and structure and function prediction[END_REF][START_REF] Carugo | Recent progress in measuring structural similarity between proteins[END_REF] . Some methods are based on numerical techniques such as geometric hashing [START_REF] Wolfson | Geometric hashing: an overview[END_REF] or spherical harmonic descriptors [START_REF] Funkhouser | A search engine for 3D models[END_REF] . A recently reported method uses so-called Zernike descriptors [START_REF] Kihara | Molecular surface representation using 3D Zernike descriptors for protein shape comparison and docking[END_REF] .

Traditionally, protein classifications have been performed manually with the aid of automated tools, and they take into account information available to biologists regarding both the functions and the phylogenetic origins of the proteins investigated. Examples of relevant databases include SCOP (Structural Classification of Proteins) [START_REF] Murzin | SCOP: a structural classification of proteins database for the investigation of sequences and structures[END_REF][START_REF] Andreeva | Data growth and its impact on the SCOP database: new developments[END_REF] , CATH (Class, Architecture, Topology, and Homologous superfamily) [START_REF] Orengo | CATH--a hierarchic classification of protein domain structures[END_REF] , and FSSP (Families of Structurally Similar Proteins) [START_REF] Holm | A database of protein structure families with common folding motifs[END_REF] to name but a few.

In order to match two distinct protein structures, there should exist a one-to-one map between their structural elements, which is called "correspondence". In addition, proper alignment of the structural elements of these proteins should be generated. A common measure that is used for this type of alignment is RMSD [START_REF] Maiorov | Significance of root-mean-square deviation in comparing threedimensional structures of globular proteins[END_REF][START_REF] Carugo | A normalized root-mean-square distance for comparing protein threedimensional structures[END_REF] . Until now, a complete geometrical comparison of two proteins has rarely been possible mainly because most proteins have different sizes and/or different numbers and types of atoms. Therefore, a complete match between an arbitrary pair of proteins is a difficult task to accomplish in general. This is why either partial or local similarity tests have frequently been used in the past [START_REF] Maiorov | Significance of root-mean-square deviation in comparing threedimensional structures of globular proteins[END_REF] . An example of using RMSD for partial similarity analysis is the STRUCTAL software [START_REF] Levitt | STRUCTAL. A structural alignment program[END_REF] . In the Results and Discussion section, we discuss in more detail different methods used for protein structure comparisons and compare and contrast them with our method In this paper, we introduce a fully automated method that enables one to compare protein structures and to perform identification of proteins. To this end we expand the protein shape function in terms of Wigner-D functions [START_REF] Wigner | Gruppentheorie und ihre Anwendungen auf die Quantenmechanik der Atomspektren[END_REF] and demonstrate mathematically that the expansion coefficients can be regarded as the structure factors of a protein. We then compare them to assess their similarity by introducing a new parameter referred to here as the "Similarity Value" (SV).

Our method obtains the similarity value in the reciprocal space (in relation to the spatial domain)

where two proteins have the same dimension (values of their structure factors). However, it is important to note that these proteins are allowed to have different numbers of atoms in the spatial domain. We demonstrated below that the SV is generally a good alternative parameter to the RMSD value. However, in comparing different-size structures, using the similarity value (SV) is strongly preferred as it permits a quantitative comparison between any protein structures independently of their sizes.

Basic mathematical idea

The Wigner D-functions describe the rotation on a sphere in 4-dimensional space (4-sphere), and they are analogous to the well-known spherical harmonic functions, which are commonly used to describe the rotation on a sphere in 3-dimensional space (3-sphere) [START_REF] Wigner | Gruppentheorie und ihre Anwendungen auf die Quantenmechanik der Atomspektren[END_REF] . A rigid body can be projected on a 4-sphere; thus, its shape function can be expanded using the Wigner-D functions.

Proteins are not typically thought of as rigid bodies due to their weak bonds, but instead they undergo sizeable thermal fluctuations at finite temperature and conformational changes due to ligand binding. However, the different conformations of a protein which are explored over time can be quantitatively characterized using shape functions in time series representations.

We start by expanding a hypothetical protein shape function, f , in terms of Wigner-D functions as

∑ ∑ ∑ (1)
where the C lmn factors are the coefficients of the series expansion, D mn l is a Wigner-D function and the parameters and satisfy: and . The Wigner-D function is defined by [START_REF] Potts | A fast algorithm for nonequispaced Fourier transforms on the rotation group[END_REF] (2)

where

√ ( 3
)
where and , and

{ (4)
while is the associated Legendre polynomial is defined as .

(

The dimension of a Wigner-D function is given by ∑

We can express Eq. 1 in matrix notation simply as . Indeed, the discrete Fourier transform on SO(3) can be written in terms of the Wigner-D functions as [START_REF] Potts | A fast algorithm for nonequispaced Fourier transforms on the rotation group[END_REF][START_REF] Hielscher | The Radon transform on SO(3): a Fourier slice theorem and numerical inversion[END_REF] ∑ ∑ ∑ ̂

where ̂ is the Fourier transform of f . We can express the above relation in matrix form as ̂ . Thus, the coefficients can be viewed as Fourier transforms of a given function .

On the other hand, we know from crystallography that the Fourier transform of the shape function of an object is defined as the corresponding structure factor [START_REF] Lipson | Fourier transforms and X-ray diffraction[END_REF] . Thus, the coefficients describe the structure factors of a given protein with the shape function (which is obtained from the positions of the atoms of the protein).

Having generated the shape function , we can obtain the coefficients of the expansion by

∫ ∫ ∫ (8)
where we use the orthogonality of the Wigner-D function:

∫ ∫ ∫ (9)

Method and algorithm

In this section we discuss practical aspects of implementing our method for particular proteins.

First, we download each protein's atom positions from the Protein Data Bank (PDB) and convert these positions to obtain the corresponding Euler angles. Then, we define the protein shape function as follows: if a voxel contains a protein's atom then is equal to one, otherwise is taken to be zero. The next steps are to compute the resultant Wigner-D functions up to and obtain the matrix elements.

One simple way to measure the similarity between two arbitrarily selected proteins is equivalent to computing the correlation value between the structure factors of the two proteins: (10) where indicates the inner product and 'abs( ) indicates the absolute value of a variable.

However, the CV measure does not provide proper comparison results for proteins, as is explained below.

Representing a 3D shape by expansion in terms of Wigner-D functions effectively projects this shape on a 3-manifold as a part of the hyper-surface of a 4-sphere. The matrix elements are the points on the manifold constructed in this manner. The CV computed in Eq. 10 gives a fractional rate of the overlap between the two manifolds.

We illustrate this with a specific example. We have chosen a crystal structure for the tubulin heterodimer with PDB code 1JFF [START_REF] Löwe | Refined structure of alpha beta-tubulin at 3.5 A resolution[END_REF] . This PDB has two subunits: 1JFF-A for the α-tubulin monomer and 1JFF-B for the β-tubulin monomer. As shown in Table I, the CV for 1JFF and 1JFF-A is approximately 1. This is because the 1JFF-A manifold is a sub-manifold of 1JFF, and all the points of 1JFF-A are subsumed by 1JFF. A discussion about SVs between 1JFF, 1JFF-A and 1JFF-B which are obtained in Table I, is given in the Results and Discussion section.

Instead of using the CV measure, we define a solid measure by applying the RMSD concept to the structure factor distances using the following procedure. A structure factor is a complex number, so we can embed it as a vector in a 2D Euclidean space. Thus, for each protein, we can define a space with the dimension equal two times the number of computed structure factors. For example, for we will have 7770 structure factors, so our space's dimensionality is (we represent this space by a 2-column and a 7770-raw matrix). Subsequently, we compute the distances between each pair of elements in this matrix.

We obtain an matrix of the distances. In a similar way we obtain another matrix for the second protein. The next step is to compute the parameter

∑ ∑ ∑ ∑ (11)
where d ij and ¢ d ij are the elements of the distance matrix of each of the two proteins. This is an RMSD relation, except we have eliminated the average coefficient 1 . Equation 11in the vector form is (12) Where ∑ ∑ , is defined similarly (these are vector lengths, i.e., the sum of the squares of arrays), and is the scalar product of the two protein vectors (i.e., the sum of the corresponding array multiplications). This scalar product indicates the correlation between two proteins, because if there is no correlation, then , and if we have a maximum correlation (the two proteins are the same), then . To obtain a direct measure of the similarity between two proteins, we define SV by rewriting Eq. 12 as follows:

⟨ | ⟩ . ( 13
)
Based on the above discussion, SV will satisfy the following inequality:

.

In other words, when the two proteins are the same, then , and when they are completely different and there is no correlation between them, then .

The Fourier transformation is a linear transform [START_REF] Arfken | Mathematical methods for physicists[END_REF] , and it preserves the length and the inner product. Thus, the Fourier transform is an isometric mapping [START_REF] Lipson | Fourier transforms and X-ray diffraction[END_REF][START_REF] Curtis | Linear Algebra. Undergraduate texts in mathematics[END_REF][START_REF] Titchmarsh | Introduction to Theory of the Fourier Integrals[END_REF] . We have shown earlier in this paper that the coefficients are the Fourier transforms of . Therefore, SV is a good measure to compare two proteins. In the following box, we summarize the algorithm for computing the similarity between two proteins in several simple steps.

Algorithm

1. Obtain protein data from the PDB website ( position coordinates of all atoms).

2. Sort atoms by their distances to the center of mass. It is assumed that all masses are distributed equally for all atoms. Another possibility could be to consider real masses of atoms.

Convert Cartesian coordinates of all atoms to the corresponding Euler angles, relative

to the center of mass of the protein.

4. Define the shape function, : the number of atoms corresponding to the position coordinates of the ith-atom in Euler angles.

Compute the results of

Eqs. 2-4.

6. Compute the structure factor, from the discrete form of Eq. 8:

∑

We sum only over occupied atom positions because the shape function is zero when there is no atom in the voxel. Thus, the other terms are zero.

7. Repeat steps 1-7 for each protein analyzed.

8. Compute SV using Eq. 13 between two proteins selected for comparison.

Results and Discussion

Table I lists CV and SV measures for selected pairs of protein structures as examples.

We see that the correlation value, CV, does not give a good comparison between two proteins. This is because it is a criterion to compute the overlap between two manifolds in the reciprocal space. If the two proteins are similar, this criterion gives a good correlation between them because these two proteins have the same structure factors. However, for two different or partially different proteins the CV is not very accurate.

To check our SV criterion, we have calculated the atomic shape function for the 1JFF-A structure by using the structure factors, . Figure 1 shows the histograms and plots of and its reconstructions ∑ for and for 1JFF-A. We see that the reconstructed functions, , are in good agreement with , especially when increases.

Figure 2 shows that when the surface under a pocket of the structure factors is normalized to one, the structure factor for a given has the Poissonian distribution:

√∑ ∑ (15)
where and . The Poissonian distribution is usually considered to be a continuous distribution. However, here we make it discrete since we need to perform a numerical computation. The maximum probability value for the Poissonian distribution occurs when and the magnitude of the corresponding peak for probability is then equal to

√ ( 16
)
where we used Stirling's approximation relation, i.e. √ . The peaks in Fig. 2 are in good agreement with Eq. 11. This is another test to confirm the validity of our method, since it gives the same result as the one obtained in x-ray pattern intensity distributions and in Poisson's distribution for random interactions between radiation and matter [START_REF] Huldt | Diffraction imaging of single particles and biomolecules[END_REF][START_REF] Wilson | The probability distribution of X-ray intensities[END_REF] .

In the following discussion we wish to highlight the differences between our method and other methods. The methods introduced for comparing protein similarities are normally based on the following the proteome-scale protein structure modeling, score function comparison, obtaining moments or descriptors, comparing RMSD between residues or chains of two proteins [START_REF] Funkhouser | A search engine for 3D models[END_REF] .

Discussing all the methods is out of scope in this paper, but here, we review some methods, which may appear similar to our method. One of these methods involves spherical polar Fourier shape density functions (SPF) [START_REF] Ritchie | Protein docking using spherical polar Fourier correlations[END_REF] . This method uses the expansion of the 3D density function in terms of radial and spherical harmonic functions and computes the correlation coefficients between two density function expansion coefficients. The other method uses Zernike descriptors.

The Zernike functions are extensions of the spherical harmonic functions. The Zernike descriptors were first used by Novotni and Klein [START_REF] Novotni | Shape retrieval using 3D Zernike descriptors[END_REF] to compare two shapes in shape searching algorithms in computer science. Later they were adapted for protein comparison purposes [START_REF] Sael | Fast protein tertiary structure retrieval based on global surface shape similarity[END_REF][START_REF] Venkatraman | Protein-protein docking using region-based 3D Zernike descriptors[END_REF] .

The 3D Zernike method is a rotational invariant method and it finds a descriptor which represents a given shape. By comparing the descriptors, the similarities between any two shapes could be determined. Another method that should be mentioned here is the spherical harmonic method [START_REF] Funkhouser | A search engine for 3D models[END_REF][START_REF] Sael | Fast protein tertiary structure retrieval based on global surface shape similarity[END_REF] . This method expands a shape function in terms of spherical harmonic functions.

After some algebraic computations, the spherical harmonic method defines the descriptors and compares them. The above methods use moment or descriptor concepts to compare proteins.

Some methods have used RMSD values as a score to compare between two structures but because of the different protein sizes normally these methods use RMSD only partially. For example, some of these methods have used a difference between the intra-structural residueresidue distances, e.g. Dali [START_REF] Holm | Protein structure comparison by alignment of distance matrices[END_REF][START_REF] Holm | Dali: a network tool for protein structure comparison[END_REF] , CE [START_REF] Shindyalov | Protein structure alignment by incremental combinatorial extension (CE) of the optimal path[END_REF] , or between inter-structural residue-residue distances, such as STRUCTAL [START_REF] Levitt | STRUCTAL. A structural alignment program[END_REF] , SAL [START_REF] Kihara | The PDB is a covering set of small protein structures[END_REF] and TM-score [START_REF] Zhang | TM-align: a protein structure alignment algorithm based on the TMscore[END_REF] . Our method is different from these structural methods for the following reasons:

1. Using the Wigner-D function does not require a definition of a radial function, as is done in the 3D Zernike, spherical harmonic or SPE methods. The three angles in the Wigner-D functions are the Euler angles and it is well known from classical mechanics that moving through a 3D rigid body is possible by using three Euler angles.

2. We show that the expansion coefficients of the shape function defined by the Wigner-D functions are equivalent to the Fourier transform of the shape function (see more details in the Basic Mathematical Idea section). Thus, we introduce the expansion coefficients of a shape function in terms of Wigner-D functions as structure factors.

3. It is well known that the Fourier transform (consequently an expansion on the Wigner-D function) is a linear transform and hence it preserves isometry [START_REF] Potts | A fast algorithm for nonequispaced Fourier transforms on the rotation group[END_REF][START_REF] Hielscher | The Radon transform on SO(3): a Fourier slice theorem and numerical inversion[END_REF][START_REF] Curtis | Linear Algebra. Undergraduate texts in mathematics[END_REF][START_REF] Titchmarsh | Introduction to Theory of the Fourier Integrals[END_REF] . Thus, if we define a RMSD-type criterion, like Similarity Value between structure factors, we prove mathematically that the properties obtained in the reciprocal space reflect the same properties in the position space. That means that if in the reciprocal space two proteins are similar, the same result holds true in the position space (provided a method of comparison is defined).

4. The size of a protein analyzed does not affect our comparison. This is because we are able to compare two proteins with the same size in reciprocal space, even though they may be different in position space, and we can choose the dimension of reciprocal space according to a desired level of accuracy. Note that in reality the expansion terms have to go to infinity but similar to other computational calculations, we should choose a cutoff in order to terminate this divergence. The number of expansion coefficients used increases the level of accuracy but also the cost of computation.

5. Contrary to inductive methods, our method is deductive and it is proved by mathematics, thus it does not rely on a great deal of experience to be validated. This is why we did not need to compare a large number of proteins in our manuscript and only a few number of known [START_REF] Li | Characterization of local geometry of protein surfaces with the visibility criterion[END_REF] . When an RMSD value between two proteins structures increases, SV decreases which means that the similarity between these two proteins decreases. Some differences between two results that can be found in the tables are related to the partial computation of RMSD value, especially when the number of atoms in the two proteins compared is different. This good agreement between RMSD values and SVs indicates that our SV is a reliable parameter for comparing any two proteins. In some algorithms, using RMSD values in a part of the algorithm, RMSD can be replaced by SV as an alternative parameter. This is because SV is equivalent to RMSD and SV can be computed more precisely than RMSD for proteins with different sizes. However, it should be kept in mind that the cost of computation for SV is about 33 seconds for comparing two proteins using a laptop with an Intel Core i7 CPU.

Below we discuss a specific example by explaining the SVs results for 1JFF and their monomers (1JFF-A and 1JFF-B). As we know, 1JFF-A has , 1JFF-B has and 1JFF has the sum of both monomers atoms i.e.

. The SVs reported between these macromolecules in Table I are:

1JFF 1JFF-A SV = 0.2091 1JFF 1JFF-B SV = 0.4967 1JFF-A 1JFF-B SV = 0.2271
Here, we discuss these results in more detail. To simplify notation, we represent 1JFF-A with A, 1JFF-B with B, and 1JFF with AB. Their distance matrices are defined by which has arrays, similarly with s ´s matrix elements, and with arrays. Note that here we define distance matrix in position space. We can write as (17) where indicates direct sum between matrices (see Fig. 3). We note that and are transpose of each other. Let us assume an unknown direct way (in position space and not in reciprocal space), then we can find the RMSD between the above structures. We then compute the following terms ) where [ ] means that all arrays of the matrix in the bracket will be squared, ∑[ ] is defined as summation over all matrix arrays in the bracket, [ ],

∑[ ] ∑[ ] ∑[ ] (18
indicates "an imaginary minus sign" between two different size matrices (the difference between matrix sizes is a serious problem when trying to define a corresponding RMSD) and we do not know how it acts. Now, we expand between AB and A and B. First, we have

∑[ ] ∑[ ] (19)
Similarly, we find that

∑ [ ] ∑[ ] (20)
We readily see that . This shows that two monomers do not have the same RMSD values or SVs. Now, one can add leaks of arrays in distance matrix of A (or B) with respect to the distance matrix size of AB by adding zeros (see Fig. 4). Note that this assumption has not been proved and is only a heuristic. Thus, the following relation could be obtained (a derivation is given by Figs. 5 and6)

| | (21)
where has a similar definition to the one mentioned after Eq. 12.

However, this definition does not provide a normalized measure to compare with the SV. Now, we define SS as follows

| | (22)
where we have defined as follows (23) and is defined similarly. Here, to arrive at the right-hand side of Eq. 22 we approximated in the dominator. We will see later this is a very reasonable approximation.

The computation of for 1JFF and its two monomers 1JFF-A and 1JFF-B results in , and in units of square Angstroms. From these values, it is easy to obtain . To get a sense of the numerical values involved, in addition to our real case, we compute SS in Eq. 22 for a totally correlated case and for an uncorrelated ( ) case. Thus, we have

{ ( 24
)
The above equation shows that . We normalize SS to ½ so that we have {

Thus, for our case we have:

. Now, we come back to SVs obtained in Table I that yields: . We see that these two results are in good agreement. Note that SS Normalized is obtained by using an approximation.

Here in the structural form, we have shown why two monomers of 1JFF are different. Thus, if one of these monomers is similar to 1JFF the other could not be similar and vice versa.

Conclusions

This paper introduces a new method to compare protein structures; it can be generalized to compare arbitrary shapes defined as a set of 3D coordinates. The novelty of our method lies in expanding the shape function using Wigner-D functions, showing that the expansion coefficients correspond to the structure factors, and using the RMSD measure in the reciprocal space (for the structure factors) to define a similarity value, namely the SV parameter. We show that this measure gives a corresponding similarity in the spatial domain because of the isometric property of the Fourier transform. We have verified our method by obtaining the shape function by using the structure factors and Wigner-D functions (see Fig. 2). The absolute values of the structure factors are the same as the intensities measured by x-ray scattering. We also show that the structure factor distribution is a Poisson distribution; as is well known, the intensity distribution for x-ray scattering is also a Poisson distribution. This result demonstrates the reliability of our method. The numerical results shown in Tables I-III for SV also confirm the reliability and usefulness of our method.

An important problem for similarity comparison methods is that the number of the protein atoms in an arbitrary pair of proteins is generally not the same. To address this problem, some methods use partial similarity measures between two proteins. However, in our method, despite the fact that the number of atoms of the two proteins being compared is different, the number of structure factors is the same in reciprocal space. This is another important advantage of our method. " " means that all arrays of the matrix will be squared. The sum over all arrays of these matrices yields as defined by Eq. 11.

Bottom. Subtraction of two top matrices. Normally, we should subtract the sum of two top matrices. But, here to show the derivation of our formula in Eq. 21 before summation we subtract two top matrices and we see the result in the bottom. To solve Eq. 21, we have to sum over all arrays of bottom matrix. The minus sign causes the change of shading on in the bottom matrix.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Histograms (a) and plots (b) of the shape function (blue) and its reconstructions for(green) and (red) for 1JFF-A. We see that the reconstructed are in good agreement with when increases.

Fig. 3 .

 3 Fig. 3. A schematic of the distance matrix representation.

Fig. 4 .

 4 Fig. 4. Extended matrix. Here we add zeros to change the size of from to .

Fig. 5 .

 5 Fig. 5. A schematic representation of obtaining . The subtraction of the two top matrices yields the bottom matrix.

Fig. 6 .

 6 Fig. 6. Top. Squared arrays of and . " " means that all arrays of the matrix

Fig. 1 Fig. 2

 12 Fig.1

Fig. 3

 3 Fig.3

Fig. 4 Fig. 5 Fig. 6

 456 Fig. 4

Table I .

 I Comparison between pairs of proteins using correlation value (CV) and similarity value (SV). 1JFF, 1SA0, 1TUB, and 1FSZ are all structures of proteins in the tubulin-FtsZ superfamily. 1ATN is a structure for actin. Comparisons between unrelated protein pairs (tubulin-FtsZ superfamily with actin) are italicized.

	First	Second	Correlation	Similarity
	Protein's	Protein's	Value (CV)	Value (SV)
	PDB ID	PDB ID
	1JFF	1ATN	0.3962	0.0002
	1JFF	1FSZ	0.9661	0.3453
	1JFF	1SA0	0.9361	0.0981
	1JFF	1TUB	0.9933	0.4872
	1JFF	1JFF-A	0.9537	0.2091
	1JFF	1JFF-B	0.9987	0.4967
	1JFF-A	1JFF-B	0.9460	0.2271
	1JFF-A	1ATN	0.5634	0.0010
	1JFF-A	1SA0	0.9955	0.3771
	1JFF-A	1TUB	0.9300	0.2425
	1JFF-A	1FSZ	0.9795	0.4131
	1JFF-B	1SA0	0.9270	0.1073
	1JFF-B	1ATN	0.3782	0.0002
	1JFF-B	1TUB	0.9969	0.4962
	1JFF-B	1FSZ	0.9637	0.3689
	1ATN	1SA0	0.5973	0.0023
	1ATN	1FSZ	0.5016	0.0005
	1ATN	1TUB	0.3639	0.0002
	1SA0	1FSZ	0.9806	0.2322
	1SA0	1TUB	0.9085	0.1154
	1FSZ	1TUB	0.9489	0.3852

Table II .

 II Set of 48 protein structures with SV and RMSD from Li et al.[START_REF] Li | Characterization of local geometry of protein surfaces with the visibility criterion[END_REF] for comparison. The SVs are computed from structure factors for l max = 7.

	1ACJ	1QIF	0.497	0.31
	1SNC	1STN	0.495	0.70
	1STP	1SWB	0.145	0.33
	First Protein's 1ULB	Second Protein's 1ULA	Similarity 0.474	0.79
	PDB ID 2YPI	PDB ID 1YPI	Value (SV) 0.165	RMSD [1] 1.27
	1A6W 2H4N	1A6U 2CBA	0.198 0.498	0.34 0.20
	1MRG 2CTC	1AHC 2CTB	0.401 0.499	0.43 0.15
	1RNE 5CNA	1BBS 2CTV	0.316 0.034	0.61 0.40
	1RBP 1FBP	1BRQ 2FBP	0.108 0.494	0.62 1.06
	1BYB 2SIM	1BYA 2SIL	0.499 0.499	0.43 0.14
	1HFC 1MTW	1CGE 2TGA	0.399 0.159	0.37 0.42
	3GCH 1APU	1CHG 3APP	0.070 0.498	1.10 0.40
	1BLH 1QPE	1DJB 3LCK	0.497 0.465	0.23 0.28
	1INC 5P2P	1ESA 3P2P	0.397 0.480	0.21 0.42
	1GCA 4PHV	1GCG 3PHV	0.499 0.045	0.32 1.23
	1HEW 3PTB	1HEL 3PTN	0.498 0.122	0.21 0.26
	1IDA 1BID	1HSI 3TMS	0.083 0.499	1.07 0.24
	1DWD 1OKM	1HXF 4CA2	0.150 0.472	0.27 0.22
	2IFB 4DFR	1IFB 5DFR	0.382 0.496	0.37 0.82
	1IMB 3MTH	1IME 6INS	0.498 0.381	0.22 1.09
	2PK4 6RSA	1KRN 7RAT	0.445 0.440	0.39 0.18
	2TMN 1CDO	1L3F 8ADH	0.266 0.403	0.62 1.34
	1IVD 7CPA	1NNA 5CPA	0.426 0.132	1.23 0.40
	1HYT 1ROB	1NPC 8RAT	0.332 0.469	0.88 0.28
	1PDZ 1IGJ	1PDY 1A4J	0.499 0.411	0.66 0.80
	1PHD	1PHC	0.499	0.17
	1PSO	1PSN	0.499	0.33
	1SRF	1PTS	0.498	0.26

Table III .

 III Set of 86 protein structures with SV and RMSD from Li et al.[START_REF] Li | Characterization of local geometry of protein surfaces with the visibility criterion[END_REF] for comparison. The SVs are computed from structure factors for l max = 7.

	6ALD 1K70 1WBL	1EWD 1K6W 1WBF	0.477 0.497 0.371	0.44 1.08 0.39
	1NLM 1M6P 1YDB	1F0K 1KEO 1YDC	0.163 0.136 0.491	1.66 1.05 0.12
	1F4X 3KIV 1H0S	1F4W 1KIV 2DHQ	0.488 0.467 0.499	0.25 0.30 0.26
	First Protein's 1JBW 1KMP 1LLO	Second Protein's 1FGS 1KMO 2HVM	Similarity 0.430 0.498 0.498	1.48 0.64 0.12
	PDB ID 1FR8 2NGR 43CA	PDB ID 1FGX 1KZ7 43C9	Value (SV) 0.498 0.467 0.491	RMSD [1] 0.54 1.61 0.23
	1AD4 1LD8 2MIN 5BIR	1AD1 1FT1 1L5H 4BIR	0.499 0.416 0.084 0.487	0.50 0.92 0.55 0.61
	1AHX 1HVC 1LL2 5EUG	1AHG 1G6L 1LL3 4EUG	0.499 0.345 0.496 0.498	0.24 0.46 0.37 0.21
	1AUR 1LSP 1LMC 5EAU	1AUO 1GBS 1LMN 5EAS	0.499 0.360 0.499 0.064	0.20 0.26 0.10 0.40
	1AXZ 1LC3 1EYN 7TAA	1AXY 1GCU 1NAW 6TAA	0.498 0.458 0.208 0.499	0.12 0.77 1.02 0.24
	1GN8 1GJW 1BHT	1B6T 1GJU 1NK1	0.491 0.499 0.033	0.51 0.29 0.58
	1B9Z 1N75 1PBO	1B90 1GLN 1OBP	0.494 0.422 0.143	0.54 1.47 0.38
	1LRI 1GOY 1OPB	1BEO 1GOU 1OPA	0.498 0.476 0.295	1.05 0.47 0.68
	1BUL 1H46 1I75	1BUE 1GPI 1PAM	0.499 0.193 0.499	0.18 0.15 0.13
	1BYD 1GUZ 1NME	1BYA 1GV1 1PAU	0.499 0.383 0.499	0.43 0.62 0.29
	1C3R 1YDD 1KEV	1C3P 1HEA 1PED	0.202 0.491 0.281	0.39 0.18 0.81
	1C5I 1YDA 1PIG	1C5H 1HEB 1PIF	0.494 0.498 0.495	0.13 0.20 0.32
	1QJW 1KIC 1PJC	1CB2 1HOZ 1PJB	0.498 0.420 0.498	0.63 0.35 0.61
	1CTE 1A80 1KLT	1CPJ 1HW6 1PJP	0.499 0.466 0.168	0.29 0.93 0.97
	1SZJ 1I3A 1QHG	1CRW 1I39 1PJR	0.499 0.498 0.499	0.33 0.40 0.23
	1ESW 4AIG 1CEB	1CWY 1IAG 1PKR	0.498 0.494 0.041	0.38 0.26 0.58
	1CY7 1JZS 2PK4	1CY0 1ILE 1PMK	0.155 0.497 0.417	1.12 0.69 0.71
	1DED 1JQ3 1BK9	1D7F 1INL 1PSJ	0.481 0.493 0.494	0.26 0.35 0.24
	1P7T 1JAY 1QBB	1D8C 1JAX 1QBA	0.406 0.435 0.497	0.66 0.60 0.11
	1DMY 1UEH 1PYY	1DMX 1JP3 1QME	0.499 0.499 0.157	0.19 0.67 0.59
	1DQY 1JSO 1OSS	1DQZ 1JSM 1SGT	0.052 0.499 0.367	0.75 0.10 0.27
	1LP6 1JYL 1SWN	1DV7 1JYK 1SWL	0.471 0.208 0.497	0.56 0.94 0.31
	1E2S 1JVS 1LBT	1E1Z 1K5H 1TCA	0.499 0.351 0.440	0.13 1.16 0.24
	1ESE	1ESC	0.499	0.19

Acknowledgements S.M.S.F. acknowledge grant number 2/22306 from Ferdowsi University of Mashhad. J.A.T. gratefully acknowledges research support received from the National Science and Engineering Research Council (Canada), the Canadian Breast Cancer Foundation, and the Allard Foundation.


Bruno Lecoutre
email: bruno.lecoutre@univ-rouen

Gérard Derzko

Comparing performances of several response-adaptive designs in dose finding studies

Keywords: Bayesian inference, dose selection, power, response-adaptive design, treatment failures

The main objectives of dose-finding studies are (i) to select the best dose(s) for future phase 3 and/or (ii) to evaluate the dose-response relationship. At that phase 2B stage of the drug development, little is known about the drug's efficacy, and only a small-to moderate-size sample can be experimented. Response-Adaptive (RA) designs, where the dose allocation to new patients depends on the information got from previous allocations and responses, do match dose-finding study conditions well. Therefore one can wonder whether they are valuable alternative designs to Balanced Randomization (BR) designs to meet the above-mentionned dose-finding objectives. Here the operating characteristics and the capability of selecting the best dose(s) and of detecting pre-specified dose contrasts are compared in four multiple dose RA designs namely two Play-the-Winner type, Generalized Drop-the-Loser and Doubly-adaptive Biased Coin designs. The restriction of these designs to two treatments have been shown to be the best choices indeed, as well as the Bayesian procedure used here for inference. The critical concept of target allocation rule is elaborated, and made explicit to match the specific dose-finding objectives. It is shown that the studied RA designs are efficient competitors to BR designs, but none of them outperforms the others on all purposes. Directions for RA design selection are given to experimenters, depending on the primary objective of their planned dose-finding study.

INTRODUCTION

Dose-finding clinical studies are phase 2 trials performed in patients with specific medical indication. They are a critical step in drug development, as they determine the design of the next large confirmatory phase 3 studies. Recently the design of dose finding studies have gained a renewed interest, as it was recognized that incorrect dose selection and poorly estimated dose-response relationship obtained from these studies may explain the high failure rate of phase 3 clinical studies. An overview and a comparative study of different designs can be found in [START_REF] Bornkamp | Innovative approaches for designing and analyzing adaptive dose-ranging trials[END_REF] and [START_REF] Dragalin | A simulation study to compare new adaptive dose-ranging designs[END_REF]. In particular practical aspects optimal designs were considered in [START_REF] Bretz | Practical considerations for optimal designs in clinical dose finding studies[END_REF], and specifically studied in [START_REF] Chambaz | Targeting the optimal design in randomized clinical trials with binary outcomes and no covariate[END_REF], or adaptive allocation rule for grouped sequential trials in [START_REF] Morgan | A comparison of adaptive allocation rules for group-sequential binary response clinical trials[END_REF]. The techniques developed here may be classified (see [START_REF] Bornkamp | Innovative approaches for designing and analyzing adaptive dose-ranging trials[END_REF], as another design-focused adaptive approach, based on empirical or optimal dose allocation rules.

Main objectives of dose finding studies are of (i) selecting the most effective dose(s) among a set of candidate doses; (ii) demonstrating dose-effect relationship. Objective (i) generally consists in selecting one or two doses for further investigations in the phase 3 trial. Objective (ii) may enhance in the demonstration of efficacy, if the dose-effect relationship is found increasing, or may detect some non monotonous dose-response profile (Gamma-shaped, bell-shaped, etc). These objectives require that a set of t doses indexed by t = 1, 2, . . . , t be studied, four doses including possibly a placebo being a minimum. This article concentrates on binary responses (success or failure), with success probability ϕ t , and number of patients n t at dose t. Note that, when the term "treatment" is used in this article, it designates some dose of one drug; doses are naturally ordered and expected to induce responses in some order as well. The spacing of the doses is another peculiar feature compared to completely distinct treatments.

At phase 2B development stage, however, very little is known about the drug's efficacy, and samples can only have small (a few tens patients) to moderate (a few hundreds patients) sizes. Therefore response-adaptive (RA) designs, where the dose allocated to the next patient depends on the already performed allocations and observed responses according to some prespecified rule, seem especially wellfitted. In these designs indeed the treatment allocation rule is devised so that the more effective doses are allocated more often in the trial, and the distribution of the total sample size N to the different dose groups is expected to account as soon as possible as soon as possible for what is learnt from the recorded effectiveness of the doses. How RA designs do perform in dose-finding studies is thus a naturally arising question.

As RA designs have been early advocated for ethical reasons, but study power is always a concern, several authors have developed optimal RA designs, which minimize the expected number of treatment failures t t=1 n t (1 -ϕ t ) while controlling the study power. A feature of paramount importance in RA designs is their so-called target allocation rule,

ρ : t → ρ t (ϕ 1 , . . . , ϕ t ) [in short ρ t ] t t=1 ρ t = 1
which depends on the probabilities of success {ϕ t }, and whose values are the asymptotic allocation proportions in each dose group t. Since the {ϕ t } are unknown and must be estimated again each time a treatment is to be allocated to a new patient, {ρ t } are only "ideal" (asymptotic) allocation proportions.

In some RA designs, the target allocation rule is univocally determined by the allocation rule of this design, and it is said to be standard, while in other RA designs, the experimenter can chose the optimal target allocation rule, so that it matches the trial objectives. In any case, given any finite sample, the closeness of the expected allocation proportions to their asymptotic values (accuracy) and the sample variance of allocation proportions (precision) are critical features, which characterize the convergence of the RA design.

The first RA designs were proposed for two-group trials, hence they are not suitable for dose-finding studies. These designs have been generalized for more than two treatments, and new designs have been developed. The RA designs can be classified according to the type of their target allocation rule, either standard or optimal.

Standard target

• Play-the-Winner [PW], introduced for 2 treatments by [START_REF] Zelen | Play the winner rule and the controlled clinical trial[END_REF] and extended to t treatments in [START_REF] Lecoutre | Adaptive designs for multi-arm clinical trials: The play-the-winner rule revisited[END_REF];

• Randomized Play-the-Winner [RPW] [START_REF] Wei | The randomized play-the-winner rule in medical trial[END_REF][START_REF] Wei | The generalized Pólya urn design for sequential medical trials[END_REF] and its Friedman's urn extensions (Bai, Hu and[START_REF] Bai | An adaptive design for multi-arm clinical trials[END_REF][START_REF] Bai | Asymptotic properties of adaptive designs for clinical trials with delayed response[END_REF];

• Drop-the-Loser [DL] [START_REF] Ivanova | A play-the-winner-type urn design with reduced variability[END_REF].

Optimal target

• Doubly-adaptive Biaised-Coin [DBC] [START_REF] Eisele | The doubly adaptive biased coin design for sequential clinical trials[END_REF];

• Generalized Drop-the-Loser [GDL] [START_REF] Sun | A generalized drop-the-loser rule for multi-treatment clinical trials[END_REF].

For two treatments, [START_REF] Rosenberger | Optimal adaptive designs for binary response trials[END_REF] identified the optimal target allocation rule, which, for a fixed asymptotic variance (reflecting the power of the test), minimizes the expected number of failures. When the rate difference ϕ 1 -ϕ 2 is chosen as the effect measure, the solution is

ρ t (ϕ 1 , ϕ 2 ) = √ ϕ t √ ϕ 1 + √ ϕ 2 (1)
to be contrasted with the Neyman allocation

ρ N eyman t (ϕ 1 , ϕ 2 ) = ϕ t (1 -ϕ t ) ϕ 1 (1 -ϕ 1 ) + ϕ 2 (1 -ϕ 2 ) (2)
which would be the target allocation rule for minimizing the variance, if no condition on treatment failures were imposed. They also derived the appropriate target allocation rules for other effect measures, which can be used as alternatives to the rate difference; yet they neither study these functions further, nor extend them to more than two treatments. [START_REF] Sun | A generalized drop-the-loser rule for multi-treatment clinical trials[END_REF] considered the "natural extension" of (??) for more than two treatments but did not justified that formula. [START_REF] Tymofyeyev | Implementing Optimal Allocation in Sequential Binary Response Experiments[END_REF] derived two equivalent optimal allocation proportions strategies based on the noncentrality parameter from a multivariate test, and described a targeted RA procedure, which provides some benefit in terms of power over balanced randomization (BR) designs. However, their proposals are difficult to implement in practice. Furthermore, they are not well adapted to the objectives of dose-finding studies. [START_REF] Basak | An urn model for odds-ratio-based response-adaptive phase III clinical trials for two or more treatments[END_REF] also considered a natural (but not justified) extension of the standard target allocation rule for odds-ratio-based urn models to more than two treatments, and compared its performance with DL and RPW with respect to a list of criteria that they defined.

Frequentist performances, namely (i) coverage probabilities of interval estimation procedures, (ii) power, and (iii) overall failure count, of PW, RPW, DL, GDL and DBC designs were thoroughly compared in the case of two treatments in [START_REF] Lecoutre | Frequentist performance of Bayesian inference with response-adaptive designs[END_REF]; they used a Bayesian approach to inference. They showed or confirmed that (i) the Bayesian approach has excellent frequentist properties, (ii) PW and DBC designs are particularly efficient, and in contrast RPW performs very poorly, (iii) the different target allocation rules have their specific interest, and (iv) the overall failure count is only modestly smaller in RA designs than in BR designs.

In the setting of dose-finding studies, where more than two doses are involved, it appears inadequate to optimize the power of detecting an overall effect, as in [START_REF] Tymofyeyev | Implementing Optimal Allocation in Sequential Binary Response Experiments[END_REF], and it is more in the line of phase 2 dose-finding studies' objectives to optimize: (i) the selection of the best dose, or the two best doses; (ii) the power of detecting a linear contrast of doses' success rates; (iii) the power of detecting a quadratic contrast of doses' success rates. In (i) only standard targets will be studied here. (ii) and (iii) can be obtained on constraint that the overall failure rate is fixed. This should be obtained by selecting the appropriate optimal target allocation rule. Yet it is still interesting and useful to include designs associated with standard targets for comparisons.

In this paper, we project to compare how well several RA designs perform and reach the abovementioned objectives, and we concentrate on two kind of designs. (i) We will consider two standard target RA designs, which extend the PW design to more than two treatments, one developed in [START_REF] Lecoutre | Adaptive designs for multi-arm clinical trials: The play-the-winner rule revisited[END_REF], and the other a new one; we exclude RPW and DL, as they performed more poorly. (ii) We will consider the GDL and the DBC designs, which can be used both with standard and optimal targets and showed excellent performances when restricted to two treatments [START_REF] Lecoutre | Frequentist performance of Bayesian inference with response-adaptive designs[END_REF][START_REF] Rosenberger | Maximizing power and minimizing treatment failures[END_REF]. As in [START_REF] Lecoutre | Frequentist performance of Bayesian inference with response-adaptive designs[END_REF] we will use non-informative Bayesian inference procedures for comparing those designs, as their (frequentist) operating characteristics were showed to compete advantageously with the frequentist inference procedures themselves in the case of two treatments, not to mention their specific asset in terms of coherence and interpretability.

The paper is organized as follows: in Section 2, the selected RA designs for more than two treatments are either defined or briefly reminded. The expressions of the asymptotic treatment allocation proportions targeted in the studied standard and optimal RA designs are obtained in Section 3; dose schemes for illustrations are also presented in that section. The Bayesian inference technique is reminded and coverage properties obtained in section 4. In section 5, capabilities of RA designs to select the best dose(s) are illustrated. Section 6 focuses on failure counts and power of detecting linear and quadratic contrasts. We conclude that none of the studied RA design is optimal for all the main objectives of a dose-finding study, but their performance depends on the particular study objective considered. A technical appendix is devoted on building the target allocation rules.

MULTI-TREATMENTS RESPONSE-ADAPTIVE DESIGNS

2.1 Play-the-Winner type designs

Basic rule for two treatments

In case of success to treatment t the next patient is assigned t ("repeat after success") and in case of failure the next patient is assigned the other treatment ("change after failure"). This ensures two characteristic properties. (i) The allocations are made according to repeated blocks of treatments. A typical sequence for two doses 1 and 2 is block 1

1 + 1 + 1 -| 2 + 2 - block 2 1 -| 2 + 2 + 2 + 2 - block 3 1 + 1 + 1 -| 2 . . .
where + indicates a treatment success and -a failure.

(ii) Sample numbers of failures in the two treatments are equal or differ by only one unit.

Direct extension of Play-the-Winner type designs: Repeated blocks [PWext]

This direct extension, which had never been studied yet, preserves the above characteristic properties.

In case of success to treatment t the next patient is assigned t. In case of failure the next patient is assigned one of the other t -1 treatments in such a way that the allocations are made according to repeated blocks of the t treatments. Consequently, the sampling numbers of failures for each treatment differ by at most one unit. The simplest possibility (which will not be considered here) is to always use the same order within blocks. A typical sequence for three doses 1, 2, 3 is block 1

1 + 1 + 1 -| 2 + 2 -| 3 - block 2 1 -| 2 + 2 + 2 -| 3 + 3 - block 3 1 + 1 + 1 -| 2 -| 3 .
. . We will prefer the following alternative: in case of failure the next patient is assigned the one of the other t -1 treatments that has the higher estimated rate of success. It is indeed a better rule with respect to the objective of allocating more often the more effective doses. A typical sequence is block 1 This alternative, studied by [START_REF] Lecoutre | Adaptive designs for multi-arm clinical trials: The play-the-winner rule revisited[END_REF], preserves the basic PW rule, but not its characteristic properties. In case of success to treatment t the next patient is assigned t. In case of failure the next patient is randomly assigned one of the other t -1 treatments, either equiprobably or proportionally to the observed rates of success (other possibilities could be considered), which implies different standard target allocation rules. It must be noted that PWraf with equiprobable assignment has the same standard target as PWext. As demonstrated by [START_REF] Lecoutre | Adaptive designs for multi-arm clinical trials: The play-the-winner rule revisited[END_REF], the convergence speed of PWraf towards the target allocation rule is excellent, boh with equiprobable or proportional assignment.

1 + 1 + 1 -| 2 + 2 -| 3 -[ordered success rates after block 1: φ1 = 2 3 > φ2 = 1 2 > φ3 = 0 1 ] block 2 1 -| 2 + 2 + 2 -| 3 + 3 -[ordered success rates after block 2: φ2 = 3 5 > φ1 = 2 4 > φ3 = 1 3 ] block 3 2 + 2 -| 1 -| 3 . . .

Randomized Play-the-Winner [RPW] and Friedman's urn extensions

The (completely) randomized Play-the-Winner rule, introduced by [START_REF] Wei | The randomized play-the-winner rule in medical trial[END_REF], is generally presented as a Generalized Friedman's Urn (GFU) model (also called as generalized Pólya urn model [START_REF] Freedman | Bernard Friedman's urn[END_REF]. A typical GFU model can be described as follows. Before the j-th patient comes in, the urn contains balls (or "particles", since the numbers of balls can be non integers) associated with each treatment t. A ball is drawn at random and replaced. The subject is assigned to the corresponding treatment (say t). When the response is known, balls are added to the urn. In the original RPW rule, one ball is added after success and no ball is added after failure. This does not preserve the basic PW rule, repeat after success and change after failure. [START_REF] Bai | An adaptive design for multi-arm clinical trials[END_REF] developed GFU extensions for t treatments that generalize in a straigthforward way the initial RPW rule. However, by contrast with PWraf, the corresponding GFU extensions of RPW perform poorly. Their convergence speed is always slower and can even be very bad [START_REF] Lecoutre | Adaptive designs for multi-arm clinical trials: The play-the-winner rule revisited[END_REF].

Drop-the-Loser type designs

The Drop-the-Loser [DL] rule, proposed by [START_REF] Ivanova | A play-the-winner-type urn design with reduced variability[END_REF], is also presented as a GFU model. DL has also relatively poor performance compared to PWraf [START_REF] Lecoutre | Adaptive designs for multi-arm clinical trials: The play-the-winner rule revisited[END_REF]. [START_REF] Zhang | Generalized drop-the-loser urn clinical trials with delayed response[END_REF] and [START_REF] Sun | A generalized drop-the-loser rule for multi-treatment clinical trials[END_REF] developed a very general model they called the Generalized Drop-the-Loser [GDL] rule. This rule involves an urn model with additional "immigration" balls.

For practical purposes, these authors recommended the following particular rule (which includes DL). Before the j-th patient comes in, if a treatment ball is drawn, that treatment is assigned, but the ball is never replaced, regardless of the response. If an immigration ball is drawn, no treatment is assigned and the ball is returned into the urn, together with u j t balls for each treatment. The total of added balls is fixed C.

For the simple DL rule, u j t = 1, which implies a standard target. In GDL, the number of added balls u j t are chosen for adjusting the desired target allocation rule. A fundamental property is that if u j t = Cρ t the sample allocation proportions converge to the target allocation rule. In practice, for each j {ρ t } are estimated from the estimated probabilities { φt }, so that u j t depends on j. The immigration balls play two roles. They allow to add treatment balls to the urn according to the current estimates of the success rates, and in return they eliminate the possibility that certain type of treatment balls will become extinct. The initial proportion of immigration balls do not yield much difference in terms of the allocation proportions [START_REF] Sun | A generalized drop-the-loser rule for multi-treatment clinical trials[END_REF].

Doubly-adaptive Biased coin Design

Eisele (1994) introduced and studied the asymptotics of the Doubly-adapted Biased Coin Design (DBCD), originally indicated when responses are independent random variables from standard exponential family. This design was considered in the case of two treatments with immediate responses by [START_REF] Rosenberger | Maximizing power and minimizing treatment failures[END_REF] (see also [START_REF] Hu | Asymptotic properties of Doubly-adaptive Biased Coin Designs for multitreatment clinical trials[END_REF].

The probability z j t of assigning treatment t to the j-th patient is a function of both the desired target allocation rule and the proportion P j t of patients who have received dose t before the j-th patient. The most frequently proposed allocation function is

j, t → z j t = ρ t ρ t /P j t γ t ℓ=1 ρ ℓ ρ ℓ /P j ℓ γ t t=1 z j t = 1
where γ is a non-negative integer, which is shown to have a limited influence on the performance of the allocation process [START_REF] Rosenberger | Maximizing power and minimizing treatment failures[END_REF].

As GDL, under general conditions for the allocation function, the sample allocation proportions converge to the target allocation rule. For each j, {ρ t } are estimated from the estimated probabilities { φt }.

TARGET ALLOCATION RULES

Standard targets

PWext target allocation rule is

ρ t = 1 1-ϕt t s=1 1 1-ϕs t = 1, . . . , t (S1)
Indeed within each block the rank of the first failure in t has geometric distribution G(1-ϕ t ), with expectation 1/(1 -ϕ t ). This design can be compared with GDL and DBC, as those two may accommodate any target allocation rule. PWraf (with equiprobable assignment after failure), a version of RPW and DL designs have that same target allocation rule, but, in our experience, RPW and DL performed more poorly than PWraf, and preliminary results indicate that PWraf was outperformed by PWext, therefore these three designs will not be included in comparisons. PWraf, with proportional assignment after failure, target allocation rule is

ρ t = ϕt 1-ϕt t s=1 ϕs 1-ϕs t = 1, . . . , t (S2)
Indeed the expectation 1/(1 -ϕ t ) must be weighted by ϕ t , because the allocation rule in this design specified proportionality to success rates. This design can be compared with GDL and DBC, as those two may accommodate any target allocation rule. A version of RPW has also the same target allocation rule, but performs more poorly and will not be included in comparisons.

Optimal targets

Optimal targets aim to control the power of detecting some specified effect, while minimizing the expected number of failures. If controlling the power of the test of homogeneity of the success rates is the objective, it can be achieved by minimizing the expected number of failures under the condition of a fixed asymptotic variance for the number of treatment failures, and the target allocation rule can be shown to be

ρ t = √ ϕ t t s=1 √ ϕ s t = 1, . . . , t (O1)
If controlling the power of the test of a specific contrast of the success rates is the objective, it can be achieved by minimizing the expected number of failures under the condition of a fixed asymptotic variance for the contrast t t=1 c t ϕ t ( t t=1 c t = 0), and the target allocation rule can be shown to be

ρ t = |c t | √ ϕ t t s=1 |c s | √ ϕ s t = 1, . . . , t (O2)
Proofs are based on Lagrange multiplier technique (see Appendix). Of note, the optimal target also depends on the effect measure of interest; optimal target for contrasts of success rates' odds, logarithms or odds logarithms are demonstrated in the appendix. As usual, optimal targets are optimal only if the underlying condition of interest is present; if not, the power may decrease considerably. GDL and DBC designs can be compared when O1 or O2 are the selected target allocation rule. But there is also some interest in assessing PWext and PWraf respectively associated with their standard target allocation rules S1 and S2; indeed this allows appreciating the gain obtained in using optimal target, compared to standard ones.

Dose schemes for illustration

In the next sections, comparisons between PWext, PWraf, GDL and DBC designs will be illustrated with 4-dose schemes. As the proposed approach to evaluate the dose-response relationship is based on contrasts, all the results will be illustrated with 3 typical scenarios involving linear and quadratic contrasts of equi-spaced doses equal to 1, 2, 3 and 4 times some unit dose u. In this situation, the coefficients {c L t } of the linear contrast are (-.3, -.1, +.1, +.3), which are such that δ L = t c L t ϕ t is the slope of the regression line. The coefficients {c Q t } of the quadratic contrast are (+.25, -.25, -.25, +.25), which are such that δ Q = t c Q t ϕ t is the second-degree coefficient of the adjusted quadratic polynomial. Note that these coefficients would be the same for a situation involving a placebo and three doses equal to 1, 2 and 3 times the unit dose. Of course, situations involving unequally spaced doses could easily be dealt with by using the appropriate coefficients. The three scenarios are illustrated in Figure 1. Scheme 1 [success probabilities .2 .4 .5 .6] shows a notable linear but no quadratic effect, scheme 2 [.25 .45 .5 .2] shows a notable quadratic but no linear effect, and scheme 3 [.1 .4 .5 .35] shows both linear and quadratic notable effects. Seven other monotonic dose schemes will be used as illustrations where contrasts are not the concern.

BAYESIAN INFERENCE

Let n t1 and n t0 be the respective numbers of successes and failures to treatment t. For all considered adaptive designs, the likelihood function is proportional to

t=T t=1 ϕ n t1 t (1 -ϕ t ) n t0 (3)
i.e. proportional to the likelihood function associated with the comparison of T independent Binomial (or Negative Binomial) proportions. A simple and usual Bayesian solution assumes T independent Beta priors for ϕ t , Beta(ν t1 , ν t0 ). This is a conjugate prior and the marginal posterior distributions are again independent Beta distributions

Beta(ν t1 + n t1 , ν t0 + n t0 ) (4)
The Bayesian approach allows obtaining the distributions of any derived parameter of interest from the joint posterior distribution [START_REF] Lecoutre | Bayesian predictive approach for inference about proportions[END_REF]: in particular, in the particular case of two treatments, efficient inferences can be made for the ratio, the difference or the oddsratio [START_REF] Lecoutre | Adaptive designs for multi-arm clinical trials: The play-the-winner rule revisited[END_REF]. The posterior distribution for any parameter of interest can be PW: play-the-winner; GDL: generalized drop-the-loser DBC: doubly-adaptive biased coin; BR: balanced randomization approximated by simulating a large sample from two independent Beta distributions, i.e. repeatedly drawing random ϕ t from their marginal posterior Beta distributions and computing the associated value for the parameter of interest. The α/2 and 1 -α/2 quantiles of the generated values give the (approximate) desired confidence limits.

In all rigor, the Jeffreys rule gives different priors for the different designs, since it is based on the Fisher information. So, it could be argued that an objective prior should depend on the design (e.g., [START_REF] De Cristofaro | On the foundations of likelihood principle[END_REF][START_REF] Bunouf | Bayesian priors in sequential binomial design[END_REF][START_REF] Sun | Objective Bayesian analysis under sequential experimentation[END_REF]. Most possible refinements of the inference would however need more complicated numerical tools for a small expected improvement in terms of efficiency. This is in agreement with the solution proposed by Berger (2004, pages 5-6) for a problem of medical diagnosis involving three proportions. While expressing his theoretical preference for the reference prior approach [START_REF] Berger | On the development of the reference prior method (with discussion)[END_REF], involving a different prior for each particular parameter of interest, he uses the simple Jeffreys prior in practice.

Coverage properties

95 percent simultaneous Bayesian credible region for φ = (ϕ 1 , . . . , ϕ t ) were simulated from 10 5 replications for the sample size N = 100. The probability of error -the region does not contain φ -is reported in Table ?? for the following designs: PWext and PWraf, two versions of DBC with the corresponding standard targets (S1 and S2), BR.

The Bayesian credibility regions have always reasonable coverage probabilities, in most cases close to the nominal level. Similar results have been obtained for all designs considered in this paper, including the case of optimal targets (with GDL and DBC designs). Interval estimates for specific contrasts between the {ϕ t } were also investigated and have also good coverage properties. This confirms the conclusions previously obtained for the Jeffreys prior in other situations (see [START_REF] Lecoutre | Frequentist performance of Bayesian inference with response-adaptive designs[END_REF][START_REF] Berger | The case for objective Bayesian analysis[END_REF][START_REF] Lecoutre | Bayesian procedures for prediction analysis of implication hypotheses in 2 × 2 contingency tables[END_REF][START_REF] Agresti | Frequentist performance of Bayesian confidence intervals for comparing proportions in 2×2 contingency tables[END_REF][START_REF] Lecoutre | Bayesian methods for experimental data analysis[END_REF].

COMPARISON OF DESIGNS AND ILLUSTRATIONS

Simulation study

The initial conditions and parameter values recommended by the authors were used (we found no clear differences with other settings). In all models, the first treatment was randomly selected. For PWext and DBC, the first four patients were randomly assigned to each of the four treatments. GDL was used with no initial treatment balls and one immigration ball. The parameter values were C = 2 for GDL and γ = 2 for DBC (PW-type designs have no parameter).

When needed, the {ϕ t } were estimated by the Bayesian posterior mean associated with the Jeffreys prior for a Binomial proportion.

Allocation proportions for standard targets

Table ?? gives, for the 3+7 doses schemes, the expectation and standard deviation of the proportion of subjects assigned to each treatment, as well as the asymptotics (target values) for the standard targets S1 and S2. Numbers in Table ?? were obtained from simulations with 10 5 replications. The sample size was fixed to N = 100.

With target S1 all compared designs have an excellent accuracy, i.e. expected sample allocation proportions nicely close to their target allocation rule values, although PWext is always better performer than GDL and DBC. The convergence speed of PWext, and in a less degree of DBC, is especially good.

For illustration, simulations for sample size N = 25 and success probabilities .5, .6, .7, .8 give the respective expected sample allocations: [.161 .199 .262 .378] (PWext) [.211 .230 .258 .301] (GDL) [.178 .210 .261 .351] (DBC) which are to be compared to the asymptotic values [.156 0.195 0.260 0.390]. The superiority of PWext also appears in terms of precision (see the standard deviations of sample allocation proportions).

Target S2 is more significant of an already noticed trade-off between precision and accuracy: GDL has better precision (standard deviation) than PWraf and DBC designs, but worsened accuracy (expected allocation proportion). Indeed, as commented in Lecoutre, Derzko and ElQasyr (2010) "the gain in variance reduction obtained with GDL is highly correlated with the loss in convergence rate". In fact, this gain appears to be an artefact, due to the poor performance of the GDL rule in terms of allocation proportion. For illustration, simulations obtained by increasing the sample size to N = 500 for the probabilities of success .5, .6, .7, .8 give the respective expected sample allocations: [.123 .178 .266 .434

Selecting doses

The performance of RA designs for dose selection is evaluated by the chance of finding a correct ordering of the allocation and success sample proportions and by the chance of observing more success for the best dose and for the two best doses.

Success rates ϕ t are not necessarily increasing with t, but it can be easily verified that, for standard targets S1 and S2, the application ϕ t → ρ t is monotonic, so that dose ordering can be consistently considered through either φ or ρ orderings. Of course, the sample estimates φt do not always have the Table 2: Allocation proportions of treatments for standard target RA designs (N = 100). For each dose scheme, target allocation rule values and expected sample allocation proportions obtained with selected RA designs (simulated from 10 5 replications).

Target S1 PW: play-the-winner; GDL: generalized drop-the-loser DBC: doubly-adaptive biased coin; BR: balanced randomization true ordering, and the application φt → ρt is not always monotonic. Hence looked-for performance measures are possibly the capability of RA design to generate as large as possible proportions of correctly ordered sample values φt and ρt . As ties may occur in sample φt and/or ρt , a correct order may be considered either in the strict sense (counted as correct orderings are only those with no tied doses), or in the extended sense (all correct orderings up to possible tied doses are considered as correct). Numbers are always larger with the latter approach. Similar issue occurs for best dose(s) selection: for example if the estimated success rates of dose 2 and 3 are maximal and identical, we can consider this selection of the best dose as incorrect, because of the tie (strict sense), or correct (extended sense). Since all considered dose schemes include only distinct doses, we will only report results for correct order in the strict sense. In consequence, special care should be taken when comparing RA designs versus BR, as the balance in the latter increases the probability of generating tied estimates.

1 1-ϕ t / t s=1 1 1-ϕs Target S2
For the dose schemes, target allocation rules and RA designs considered in Table ??, Table ?? gives the following statistics (simulated from 10 5 replications): the expected proportion of treatment failures (standard error); the proportions of sample with correct (strict) ordering of the allocation proportions and of success proportions; the proportion of samples with more success for the best dose and for the two best doses.

The total numbers of failures are quite similar, and minimal, in the three designs compared under target S2; in contrast with what has been observed with only two treatments(see [START_REF] Lecoutre | Frequentist performance of Bayesian inference with response-adaptive designs[END_REF][START_REF] Rosenberger | Maximizing power and minimizing treatment failures[END_REF], the gain in the failure counts obtained with RA designs is notable compared with the one obtained with an equally sized balanced randomization design (BR). A similar observation is also valid, yet to a lesser degree, for the three designs compared under target S1. Of note in these two cases, GDL design has slightly poorer performance than PW-types or DBC designs. Accuracy in failure counts is quite similar in all six RA designs. The picture is quite different in terms of ranking correctly the allocated proportions or the success proportions: PWext outperforms all the other designs, including BR. Finally, the PWext shows a very significant advantage over all others for selecting the best dose, and a significant advantage over all the other RA designs for selecting the two best doses.

In order to refine the comparisons of correct rankings of allocations and successes, the sample sizes N required to get 80 per cent chance of selecting the best dose (i.e. such that 80 per cent of simulated samples have a strictly higher estimated success rate for the best dose) were estimated for each compared design.

For this sample size, Table ?? gives the following statistics (simulated from 10 5 replications): the expected number of failures; the expected number of sample with correct ordering of allocation proportions and of success proportions; the proportion of samples with more success for the two best doses. Results are only reported for the dose schemes involving the linear and quadratic contrasts.

The sample size required is generally minimal for DBC design (target S2), which also outperforms others in terms of failure counts (very notably so compared BR). Second performer for sample size and failure count is the PWext design, which outperforms other RA designs for obtaining the correct ordering of allocated doses and successes. BR remains the best performer for selecting the two best doses, and the six RA designs perform as well as each other under these conditions. Table 3: Performance of RA designs for selecting doses (N = 100). For each dose scheme, line 1: expected proportion of treatment failures (standard error), line 2: proportion of samples with correct ordering of the allocation proportions / of the success proportions, line 3: proportion of samples with more success for the best dose / for the two best doses (simulated from 10 5 replications).

Table 4: Comparisons of the ordering of allocations and successes. N is selected in order to get 80 per cent chance of selecting the best dose. For each dose scheme, line 1: expected number of treatment failures, line 2: proportion of samples with correct ordering of the allocation proportions / of the success proportions, line 3: proportion of samples with more success for the two best doses (simulated from 10 5 replications).

Target S1

NUMBER OF TREATMENT FAILURES AND POWER FOR CON-TRASTS

For each target, we consider the power for the two contrasts of interest, the linear contrast δ L with coefficients (-.3, -.1, +.1, +.3) and the quadratic contrast δ Q with coefficients (+.25, -.25, -.25, +.25).

Suppose we aim to demonstrate that δ L > 0. With the Bayesian approach, it can be concluded that δ L > 0 when the posterior probability Pr(δ L > 0 | data) is larger than a given 1 -α (the lower limit of the 100(1 -2α) per cent credible interval for δ L is larger than zero). The sampling probability of this event, assuming true values ϕ t * , evaluates the power of the procedure. The same computation applies to the quadratic contrast, the aim being in this case to demonstrate that δ Q < 0. Power arguments should be balanced against the performance in terms of allocation proportions and consequently in terms of expected total numbers of treatment failures.

Both the standard and optimal targets are considered here. Since the coefficients c Q t are equal in absolute value, the optimal target for the quadratic contrast O2 Q happens in this particular case to be identical to the optimal target O1, which controls the power of the test of homogeneity of the success rates. For each target, results for the best designs are given in Table ??. The expected sample allocation proportions and the expected power (probability of the event above, averaged over samples) are reported. The expected standard deviation of the posterior distribution of each contrast is also reported; it evaluates the precision of inference.

For target O2 L , optimal for the linear contrast, the allocation proportions of the smallest and the largest doses are the largest, while those of intermediate doses are smaller: this is a reminder that O2 L is optimal if the linear contrast is truly notable; but if this condition happens to be wrong, the target may be far from optimal. In contrast allocation proportions for standard target S1 or S2, and incidentally here for O2 Q (quadratic contrast) are monotonous in ϕ t , so that the information collected Table 5: Power study for linear and quadratic contrasts (N = 100, 25 subjects for each dose in BR). For each dose scheme, lines 1-4: expected sample allocation proportions, line 5: expected number of treatment failures (standard error), lines 6 and 7: estimated power (standard error of the posterior distribution of δ), respectively for the linear and quadratic contrasts (simulated from 10 4 replications). O1: optimal target for the test of homogeneity O2 L : optimal target for the linear contrast; O2 Q : optimal target for the quadratic contrast is all the more important, since the dose is truly more effective; in consequence S1 and S2 can be expected to be robust against discrepancy to any particular dose-response relationship assumption. In addition the failure counts are notably smaller with the standard targets S1 (PWext) and even more so with S2 (DBC), versus the two optimal targets, or versus BR. As expected O2 L gives the best power when the linear contrast is truly notable, and so does O2 Q (or O1) when the quadratic contrast is truly notable; the power may even be larger than BR's one sometimes.

O1 √ ϕ t / t s=1 √ ϕ s Target S1 S2 O2 Q |c Q t | √ ϕ t / |c Q s | √ ϕ s O2 L |c L t | √ ϕ t /

CONCLUDING REMARKS

In this article, we focused on RA design's operating characteristics, which match the particular objectives of phase 2B dose-finding studies well, namely dose selection and dose-response relationship. Restricting attention to success-failure responses, four promising multi-treatment RA designs were selected on the basis of the excellent operating characteristics of their restriction to two treatments; noticeably, based on extended experience, RPW and DL have been discarded as less performing. Then appropriate standard or optimal target allocation rules have been elicited; finally the general (coverage, convergence), and specific to dose-finding (ability for dose selection, power of detecting contrasts and associated failure count) properties of the selected RA designs have been compared. All the considered designs, PWext, PWraf, GDL and DBC designs had close-to-nominal coverage probabilities, using the well-tried Bayesian approach to confidence. It was found that in PW-type and DBC (with the same target allocation rule) designs, the sample allocation proportions converge excellently to their targets, especially the PWext design.

It was demonstrated that none of the studied RA design is optimal for all the main objectives of a dose-finding study, but their performance depends on the particular study objective considered: if the main objective is to select the best dose(s) for further phase 3 trial, the PWext outperforms the competitors, yet optimal designs minimize the numbers of failures in these conditions. If the objective is to infer linear or quadratic contrasts, optimal designs perform the best as expected. Yet, the required allocation proportions may not reflect the true success rates; this exposes to the risk of decreased power if the dose-effect relationship had been misspecified. In addition they perform worse than standard designs in terms of expected number of failures, and the relevant target allocation rule depends on the effect measure used for inference. Therefore a more robust option may consist in using the slightly less performing standard designs, whose allocations proportions always reflect the success rates correctly. Of note the conclusion that RA designs generally provide a modest advantage to BR in the case of two treatments (see [START_REF] Lecoutre | Frequentist performance of Bayesian inference with response-adaptive designs[END_REF][START_REF] Rosenberger | Maximizing power and minimizing treatment failures[END_REF]) is not confirmed with more than two treatments: sizeable savings of failures can even be obtained in some settings.

Comparison with the most usual BR should not be skipped: the main argument in favour of RA designs is that doses are to be allocated in conformity with some specific phase 2B objective, but one generally expects BR to have better power. In this article the discrepancy from balanced samples obtained with the diverse RA designs, and any possible worsening of power in comparison to BR can be evaluated. It was found that the best RA designs are at least equally powerful as BRs, while failure counts are smaller. When the main objective is to select the best dose(s), RA designs seem to be interesting competitors to BR group sequential designs, where involved techniques may be used to drop the worst dose(s) early.

The two extensions of PW rule to more than two treatments, PWext and PWraf, are of course subjected to potentially biasing the dose allocation process: if the trial takes place in a single experimental center indeed, any success with one dose is followed by the same allocated dose. In multi-center with centralized dose allocation, or when the evaluations are independent of the inclusion process, any selection bias becomes very unlikely. But those two designs are useful mainly as references for design comparisons. It was shown that DBC designs perform excellently, but cannot be suspected to selection bias: they can be used when complete protection of selection bias is required. The implementation of RA designs is more involved than that of BR, as in particular a centralized interactive system is required for dose allocation and response recording, and the drug supply management is more difficult. It is not the place in this article to develop these important practical aspects, but a detailed analysis of them can be found in [START_REF] Quinlan | Barriers and opportunities for implementation of adaptive designs in pharmaceutical product developmentBarriers and opportunities for implementation of adaptive designs in pharmaceutical product development[END_REF].

GDL and DBC designs cope with delayed responses. PW designs can be adapted for this circumstance. Either some default random allocation, or preferably an allocation based on available responses can be used when the adaptive rule cannot be applied.

APPENDIX: DERIVATION OF THE TARGET ALLOCATION RULE

  2.1.3 Alternative generalization of Play-the-Winner type designs: Randomized after failure [PWraf ]

Figure 1 :

 1 Figure 1: Three doses schemes involving equi-spaced doses equal to 1, 2, 3 and 4 times the unit dose u. The coefficients of the linear and quadratic contrasts are respectively (-.3, -.1, +.1, +.3) and +.25, -.25, -.25, +.25).

Table 1 :

 1 Simultaneous 95 percent Bayesian credibility region for ϕ t : Probability of error for N = 100 (simulated from 10 5 replications).

	Target	S1 1 1-ϕt / t s=1	1 1-ϕs	S2 ϕt 1-ϕt / t s=1	ϕs 1-ϕs
	ϕ 1 ϕ 2 ϕ 3 ϕ 4 PWext GDL DBC PWraf GDL DBC	BR
	.2 .4 .5 .6	.045	.052 .048	.048	.052 .048	.049
	.25 .45 .5 .2	.054	.057 .053	.047	.042 .053	.054
	.1 .4 .5 .35	.045	.045 .047	.051	.053 .047	.056
	.1 .15 .2 .3	.043	.047 .046	.030	.036 .035	.058
	.1 .2 .3 .4	.051	.050 .051	.043	.044 .041	.046
	.2 .3 .4 .5	.057	.062 .055	.045	.049 .046	.051
	.5 .6 .7 .8	.052	.051 .052	.061	.060 .066	.050
	.4 .6 .8 .9	.049	.054 .050	.046	.054 .045	.045
	.6 .7 .8 .9	.060	.054 .056	.057	.070 .060	.046
	.7 .8 .85 .9	.061	.055 .058	.067	.065 .074	.056

  Let {n t } t t=1 the sizes of t treatment groups and N = t t=1 n t the total sample size. Assume the overall failure count t t=1 n t (1 -ϕ t ) is to be minimized, while keeping constant (= K) the variance t of a contrast {c t }, with t t=1 c t = 0. Function g(.) depends only on the effect measure f (.) and on ϕ t . The table below shows a few examples.Let λ be a Lagrange multiplier, the function L(n t , λ) =

	• Contrasts of rates		n t N	=	|c t | √ ϕ t t s=1 |c s | √ ϕ s
	t=1 c 2 t • Contrasts of odds g(ϕt)		n t N	=	√ ϕt (1-ϕt) 2 √ ϕs s=1 |c s | |c t | t (1-ϕs) 2
	Rate n t N Odds • Contrasts of odds logarithms (log(odds)) • Contrasts of rate logarithms Arithmetic	=	f (ϕ t ) ϕ t |ct| √ ϕ i t s=1 √ ϕs |cs| ϕt 1-ϕt	g(ϕ t ) ϕ t (1 -ϕ t ) ϕt (1-ϕt) 3
			Rate n t N =		log(ϕ t ) |ct| (1-ϕt) √ ϕt t s=1 |cs| (1-ϕs) √ ϕs	1-ϕt ϕt
	Geometric
			Odds log( ϕt 1-ϕt )	1 ϕt(1-ϕt)
		t t=1	n t (1 -ϕ t ) + λ	t t=1	c 2 t	g(ϕ t ) n t	-K
	is to be minimized in n t 's. One gets
	∂L(n t , λ) ∂n t ∂L(n t , λ) ∂λ	= (1 -ϕ t ) 1 -= t i=1 c 2 t g(ϕ t ) n t -K = 0 λg(ϕ t )c 2 t (1 -ϕ t )n 2 t	= 0
	and solutions
			n t = |c t |	λg(ϕ t ) 1 -ϕ t
	The asymptotic treatment allocation rates are
	N =	t t=1	n t =	t t=1	|c t |	λg(ϕ t ) 1 -ϕ t
	n t N	=		|c t | g(ϕt) 1-ϕt t s=1 |c s | g(ϕs) 1-ϕs
	In particular

nt

Target S1