Sergei Grudinin
email: sergei.grudinin@inria.fr

Maria Kadukova

Andreas Eisenbarth

Simon Marillet

Frédéric Cazals
email: frederic.cazals@inria.fr

Predicting binding poses and affinities for protein -ligand complexes in the 2015 D3R Grand Challenge using a physical model with a statistical parameter estimation

Keywords: protein-ligand docking, machine learning, scoring function, ridge regression, parameter estimation

The 2015 D3R Grand Challenge provided an opportunity to test our new model for the binding free energy of small molecules, as well as to assess our protocol to predict binding poses for protein-ligand complexes. Our pose predictions were ranked 3-9 for the HSP90 dataset, depending on the assessment metric. For the MAP4K dataset the ranks are very dispersed and equal to 2-35, depending on the assessment metric, which does not provide any insight into the accuracy of the method. The main success of our pose prediction protocol was the re-scoring stage using the recently developed Convex-PL potential. We make a thorough analysis of our docking predictions made with AutoDock Vina and discuss the effect of the choice of rigid receptor templates, the number of flexible residues in the binding pocket, the binding pocket size, and the benefits of re-scoring.

However, the main challenge was to predict experimentally determined binding affinities for two blind test sets. Our affinity prediction model consisted of two terms, a pairwise-additive enthalpy, and a non pairwise-additive entropy. We trained the free parameters of the model with a regularized regression using affinity and structural data from the PDBBind database. Our model performed very well on the training set, however, failed on the two

Introduction

The Community Structure-Activity Resource (CSAR) was created in 2008 from the desire to provide high-quality experimental data to the community that is developing new computational tools. CSAR organized the first benchmark exercises in 2010 that highlighted the difficulties to predict the binding free energy of protein-ligand complexes, and discovered that some properties such as hydrogen bonding cause more trouble than the ligand size [1]. The goal of the 2011-2012 blind exercises was to compare improvements for pose prediction, enrichment, and relative ranking of congeneric series of compounds with a set of four protein targets. One of its conclusions was that correct ranking and correct pose prediction are not necessarily correlated [2]. This led to the CSAR Benchmark Exercise 2013 and 2014, which focused on the further improvement of the ability of computational methods to correctly predict binding poses and affinities. The successor of CSAR, the Drug Design Data Resource (D3R, www.drugdesigndata.org), is a new effort for docking and scoring data hosted at UCSD. It serves the same purpose as CSAR of providing high-quality experimental data for testing and improving ligand-protein docking algorithms and their scoring protocols.

The 2015 D3R Grand Challenge provided two unpublished datasets of protein-ligand co-crystal structures with the measured affinities to assess the current algorithms for pose and affinity predictions. The first dataset, referred to as HSP90, contained 8 crystal structures of the ATP site of the Hsp90 protein with resolution < 2.0 Å, as well as binding data for 180 compounds across five orders of magnitude and three chemical series. The challenge for this dataset was to predict the binding poses of 6 ligands spanning all three chemical series, and also to predict binding affinities for these 6 ligands as well as for the other 174 ligands. The second dataset, referred to as MAP4K4, comprised 18 compounds and crystal structures of a kinase protein, which is a member of the STE20 family, with resolution < 2.5 Å and inhibition/binding data over four orders of magnitude, in many cases confirmed by multiple assay methods. The challenge for this dataset was to predict the poses of 30 co-crystal structures containing diverse chemical series, and also to predict binding affinities of 18 of the ligand structures. Here, more co-crystal structures than affinity data were available.

We participated in the 2015 D3R Grand Challenge mainly to assess our new developments on the binding affinity predictions for small molecules, but also to validate our pose prediction protocol introduced earlier [3]. As we said, our main goal was to assess a new approximation to the binding free energy given by a combination of two terms, a pairwise-additive enthalpy, and a non pairwiseadditive entropy, which in principle consists of conformational entropy, vibrational entropy, and solvent entropy. Below we provide a detailed explanation of our model with the results obtained on the training set and the two test sets provided by the 2015 D3R Grand Challenge.

Test Datasets and Organization of the Challenge

The D3R Grand Challenge 2015 consists of two datasets. The MAP4K4 (protein kinase kinase kinase kinase) dataset was provided by Genentech [4]. It consists of 30 ligand compounds, for which co-crystal structures are available. For 18 of these complexes, there are also two experimental values provided in the results: IC 50 inactive and K i active.

The HSP90 (Heat Shock Protein 90) dataset was provided by AbbVie and Prof. Heather Carlson and CSAR colleagues at the University of Michigan and consists of 180 ligand compounds of three chemical classes (aminopyrimidines, benzimidazolones, benzophenone-like). Experimental IC 50 are available for all of these. For 8 of them, co-crystal structures are provided.

The 2015 D3R Grand Challenge was oraganized in two stages. For the HSP90 dataset, the oraganizers provided as input 4 protein-ligand co-crystal structures drawn from the PDB and 2 proteinligand co-crystal structures from the blinded dataset that were solved with compounds of each chemical class in three series. The oraganizers also provided SMILES strings of the 6 ligands to be docked, of the 2 ligands in the revealed co-crystal structures, and of the additional 172 compounds for affinity prediction or ranking. For the MAP4K4 dataset, the oraganizers provided as input 2 protein-ligand co-crystal structures drawn from the PDB with SMILES strings of the 30 ligands to be docked, among which 18 were selected for affinity prediction or ranking. In the first stage, the goal was to predict binding poses and binding affinities for the selected ligands. After the end of the first stage, the organizers made available the blinded co-crystal structures for both datasets. In the second stage, the goal was to predict binding affinities for the selected ligands using the released co-crystal structures of protein-ligand complexes. After the end of the second stage, the experimental binding affinities for both datasets were made available.

At the first stage of the Challenge, we submitted one docking pose prediction for each of the datasets. At this stage, we additionally submitted nine binding affinity predictions for the MAP4K4 dataset. All of these binding affinity predictions were computed using docking poses selected by nine different strategies.

At the second stage of the Challenge for each of the datasets, we submitted four binding affinity predictions using different prediction models. These were computed using co-crystal structures for the MAP4K4 dataset and top-one docking poses for the HSP90 dataset.

Binding Affinity Model

Free Energy Approximation

We approximated the logarithm of the binding affinity with the binding free energy, which is given as

∆G bind = ∆H -T ∆S, (1)
where ∆H is the enthalpic difference between the bound and the unbound states of the protein-ligand complex, T is the temperature, and ∆S is the entropic difference upon binding. The binding free energy is a much more sophisticated function compared to the potential energy, and involves not only interaction energy between the partners, but also changes in the internal energy of the compounds, interactions with solvent, rearrangement of solvent molecules and changes of conformational degrees of freedom corresponding to the entropic loss upon binding. Direct computation of the binding free energy of a molecular complex, or of the relative binding affinities associated with two ligands (i.e., the binding energy difference) can be carried out with thermodynamic integration techniques. However, while this method can be applied to complexes consisting of few partners of small size, it reaches a prohibitively high computational cost, unless specific properties of the system(s) studied can be exploited [5,6]. For example an improved relative binding affinity calculation was recently proposed, using free energy perturbation / replica exchange with solute tempering (FEP/REST, [7]). In this method, only a small region of interest is heated up, so that a small number of replicas is sufficient to increase sampling and reduce the effects of quasi-nonergodicity. However, even if the selection of this region can be automated to some extent [6], this task in non trivial in all generality. Therefore, in this work we adopted the following computational strategy. First, we approximate the enthalpic contribution ∆H using a pair-wise distance-dependent knowledge-base potential function that has been successfully applied to proteinprotein and protein-ligand scoring [3,8,9], and whose coefficients are obtained from knowledge-base using a machine-learning technique. Second, we approximate the entropic contribution T ∆S with a set of geometric descriptors. These represent the non pairwise-additive part of the binding free energy [10] and have been recently demonstrated to perform very well on the protein-protein structureaffinity benchmark [11]. Below, we provide a more detailed description of the two contributions into the binding free energy.

Enthalpic Contribution

We represent the enthalpic contribution to the binding free energy as a linear functional F defined for all possible configurations of a protein-ligand complex [3,9]. To simplify the shape of the functional, we make several assumptions. First, we assume that F depends only on the interface between the protein and the ligand. We defined the interface as a set of atom pairs at a distance smaller than a certain cutoff distance r max such that the first atom in each pair belongs to the protein and the second atom in each pair belongs to the ligand. We use the value of the cutoff distance of 10 Å, which was adapted from previous works [12,13,14,15,16,17,18,19]. It has been previously demonstrated that even smaller distances of 6 Å work somewhat well and that the quality of the scoring function does not improve after the cutoff distance of 12 Å. Second, we assume that both a protein and a ligand can be represented as a set of discrete interaction sites located at the centers of atomic nuclei. We divide all interaction sites into M types according to the properties of corresponding atomic nuclei, such as element symbol, aromaticity, hybridization, and polarity, resulting in a total of M ×(M +1)/2 pairs of interactions. An interaction between two sites can be also regarded as an interaction between two atoms of certain types. Third, we assume that F depends only on the distribution of the distances between the interaction sites, i.e. the number of site pairs at a certain distance, F (P ) ≡ F (n 11 (r), .., n kl (r), .., n M M (r)) ≡ F (n(r)), (2) where n kl (r) is the number density of site-site pairs at a distance r between two sites k and l, with site k on the protein, and site l on the ligand, where M is the total number of interaction types. Finally, we assume that F is a linear functional, such that we can write the approximation to the entalpic contribution as

∆H ≈ M k=1 M l=k rmax 0 n kl (r)U kl (r) dr, (3)
where U kl (r) are unknown functions that are deduced from the training set of binding affinities for protein-ligand complexes. From now on, we will call these functions scoring potentials. We compute the site-site number densities for each protein-ligand complex n kl (r) as a sum of Gaussian-distributed distances between interaction sites of types k and l as

n kl (r) = uv e -(r-r ij ) 2 2σ 2 , ( 4
)
where σ is the standard deviation, which is constant for all distributions. The sum is taken over all pairs of sites uv of types k and l separated by a distance r ij smaller than a certain threshold distance r max , with the first site located on the protein, and the second site located on the ligand. In order to determine unknown scoring potentials U kl (r) (Eq. 3), we decompose them along with the number densities n kl (r) in a polynomial basis,

U kl (r) = q w kl q ψ q (r), r ∈ [0; r max ] n kl (r) = q x kl q ψ q (r), r ∈ [0; r max ], (5)
where ψ q (r) are orthogonal basis functions on the interval [0; r max ], and w kl q with x kl q are the expansion coefficients of U kl (r) and n kl (r), respectively. Given this, the enthalpic contribution ∆H can be expanded up to the order Q as

∆H ≈ M k=1 M l=k Q q w kl q x kl q = (w • x), w, x ∈ R Q×M ×(M +1)/2 (6)
We will refer to the vector w as to the scoring vector, whose value is to be determined from the knowledge-base, and to the vector x as to the structure vector that is computed from the structural data. For the estimation of binding affinity, we used the maximum expansion order Q = 10 and the number of types M = 48, as it is explained below. Thus, the dimensionality of the scoring vector w and the structure vector x is Q × M × (M + 1)/2 = 11, 760.

Atom Types of the Enthalpic Contribution

Our knowledge-based enthalpic contribution to the binding free energy requires definition of the atom types. We assigned the types to the atoms of both proteins and small molecules according to their physical and chemical properties as well as their functional groups. To do so, we started with 164 internal atom types inspired by the fconv library [20] and grouped them by measuring the statistical similarity of distance distribution functions between different atom types in the training data set. Overall, our parameterization consists of 48 atom types.

More precisely, we have 17 types for nitrogen, 9 types for oxygen, 8 types for carbon, 4 types for sulphur, 2 types for phosphorus, and 8 types for halogens. Table S1 from Electronic Supplementary Material lists all the used atom types along with the chemical properties of the corresponding atoms.

We should specifically add that our model has no explicit hydrogen atoms. Generally, it is not a trivial task to correctly predict positions of all missing hydrogens. For example, their positions depend on experimental parameters of the system, such as pH. Neither our model has directional bonds. Introducing explicit hydrogen atoms or angular dependency would add additional degrees of freedom (or features) to our model, thus increasing the size of the feature space.

Entropic Contribution

Of particular interest when modeling macro-molecular complexes are terms meant to code entropic properties, in particular conformational entropy, vibrational entropy, and solvent entropy. To model the latter two terms, we use geometric features recently introduced in [10], which proved instrumental to obtain state-of-the-art affinity predictions for the protein -protein complexes of the structureaffinity benchmark [11]. In a nutshell, these descriptors dissect the surface and volume properties of the isolated subunits and of the complex whose affinity is to be estimated. Their rationale is to replace the buried surface area (BSA), namely the portion of the exposed surface of both partners which gets buried upon complex formation [21], with a more informative descriptor. While the BSA exhibits a remarkable correlation with various biophysical quantities [22], it does not account for interface shape and packing properties, these latter attributes being of interest to model entropic changes.

Interface definition. We define interface atoms using the parameter free Voronoi interface model. It is based on the α-complex derived from the Voronoi (power) diagram of the atoms, and has been proposed in [23,24]. In a nutshell, define the restriction of an atom as the intersection between its ball in the solvent accessible model and its cell in the Voronoi diagram. The Voronoi interface identifies pairs of neighboring restrictions, such that each pair involves either two different partners or a partner and the interfacial solvent. The atoms found in at least one such pair are denoted I and their complement I C .

Inverse volume-weighted internal path length (IVW-IPL) and variants. This parameter encodes the size and morphology of the interface and takes atomic packing into account. Its definition involves three ingredients.

The first one consists in the interface atoms I contributing to the Buried Surface Area, identified with the aforementioned Voronoi based model [24]. Once identified, the interface atoms of a partner define its binding patch (Fig. 1(A)). The second one is the notion of shelling order (SO), namely an integer value assigned to the interface atoms of a partner, and measuring the distance to the bounding patch boundary (Fig. 1(B)). The shelling order of atom a is denoted SO(a). The third one is the packing of an atom. Atomic packing is defined for an atom in a structure which may be the isolated subunit containing this atom or the complex defined by these subunits. Following classical works [START_REF] Gerstein | The international tables for crystallography[END_REF], we define the atomic packing as the volume of the solid defined by the intersection between the ball of the atom (in the solvent accessible model) and the Voronoi region of this atom in the structure (the isolated subunit or the complex) ( [START_REF] Cazals | [END_REF] and Fig. 1(CD)). This volume is denoted volume bound(a) for atom a.

Using these notions, the IVW-IPL is defined as follows:

IVW-IPL = a∈I SO(a) volume bound(a) (7)
On the one hand, the shelling order refines socalled core-rim models [21]. Borrowing to the notion of cooperative effects involving non-bonded weak interactions, an isotropic or disk-like interface is indeed expected to be more stable than an elongated one-even if their surface areas match.

On the other hand, the atomic packing encodes the local density of neighbors of a given atom, and thus provides a measure for local interactions (hydrogen bonds, van der Waals interactions). Note that packing is a subtle quantity related to the enthalpy -entropy compensation [27,28], as its properties strike a balance between enthalpy (a large number of neighbors favors interactions) and entropy (by restricting the moving freedom of atoms, tight packing is detrimental for dynamics yielding an entropic penalty assuming they were free to move in the first place).

We also use two variants of the previous quantity. First, the inverse volume sum which is the same as IVW-IPL except that SO(a) is replaced by 1 in Eq. (7). Second, the number of atoms at interface which can be written as in Eq. ( 7), but replacing SO(a) /Vol(a) with 1. These two quantities essentially isolate the contributions of the numerator and denominator in Eq. ( 7) so that the model parameters can be adjusted separately.

Remark 1

The previous quantities were computed using the following programs from the Structual Bioinformatic Library (http://sbl.inria.fr > Applications): sbl-vorshell-bp-ABW-atomic.exe [24] for binding patches and atomic shelling orders 1 , and sbl-vorlume-pdb.exe [START_REF] Cazals | [END_REF] for atomic volumes and surface areas2 .

Solvent Interactions and Electrostatics. In order to account for the interaction between the partners and water molecules, we use the fractions of charged and polar residues on the non interacting surface (NIS) [29], NIS polar and NIS charged . The NIS consists of the exposed surface of the partners not involved in the interface (Fig. 1(A), I C and SASA > 0).

We also include the atomic solvation energy, σ, from Eisenberg et al [30], describing the free energies of transfer from 1-octanol to water per surface unit ( Å2 ). This value is positive for non-polar atoms, and negative for polar atoms. The corresponding variable, ATOM SOLV, is a weighted sum of atomic solvent accessible surface areas (Eq. ( 8)), and may be seen as the atomic-scale counterparts of NIS charged and NIS polar :

ATOM SOLV = a∈I C SASA(a) • σ(a) (8)
Finally, we include an intermediate-grained description of the non-interacting surface which consists in the total area of polar atoms of the complex (according to Eisenberg's solveation parameters). The corresponding term, POLAR SASA (Eq. ( 9)), is also a weighed sum of exposed areas.

POLAR SASA = a∈I C and σ(a)<0 SASA(a) (9)
Using the introduced interactions terms, we represent the entropic contribution to the binding free energy as a linear combination of the following form,

T ∆S ≈ w 1 IVW-IPL + w 2 NIS polar + w 3 NIS charged +w 4 ATOM SOLV + w 5 POLAR SASA (10) Similar to the vector of free parameters w of the enthalpic term, we will deduce free parameters w 1 ...w 5 using a machine-learning procedure, as explained below.

Remark: Note that all the atomic descriptors except NIS polar and NIS charged are computed for both partners. Since the ligand is not made of aminoacids, these two residue-level descriptors are computed for the protein receptors only.

Training Set

To train our model of binding free energy given by Eqs. 6 and 10, we collected structural and binding affinity information from the PDBBind [32,33] database, which provides experimentally determined protein-ligand complexes deposited in the Protein Data Bank supplied with the measured binding affinity data. To construct the PDBBind database, its authors manually examined the primary references for each protein-complex and collected experimentally determined binding affinity data (K d , K i and IC 50 ). These constitute the "general set" of the database consisting of 14,620 complexes. Then the authors of PDBBind additionally compiled the "refined set" to select protein-ligand complexes with better quality out of the general set. To do so, they applied a number of filters to the selection regarding binding data, crystal structures, as well as the nature of the complexes. For the training, we used 3,706 complexes from the latter "refined set" of PDBBind release 2015, which contains three-dimensional structures of resolution equal to or better than 2.5 Å with the corresponding binding data, which includes K d , K i , and IC 50 values.

Training Models

Entropy only model. We trained a simple model consisting only of the five variables encoding the entropic contribution. Because of its small size, no regularization was needed and we resorted to ordinary least-squares regression to compute the parameters. Namely, the following function is minimized:

C(w) = i y i -w T x i 2 , (11)
where w i , 1 ≤ i ≤ 5 are the parameters corresponding to variables x i , 1 ≤ i ≤ 5, and y i are experimentally measured log K d /K i values.

Enthalpy only and combined model. Our model for the binding affinity contains 11,765 free parameters (11,760 and 5, respectively introduced in sections 3.2 and 3.3) that were to be deduced from 3,706 experimentally measured data points. Because the dimensionality of our problem is several times larger compared to the size of the training set, to deduce the free parameters we chose a regression model with regularization. After several initial experiments with different regression models, we stopped at the linear ridge (LRR) and kernel ridge regressions (KRR), mainly because of their speed for our data [START_REF] Barber | Bayesian reasoning and machine learning[END_REF]. The linear and kernel ridge regressions minimize the following cost function, which consists of the regression term and the regularization term,

C(w) = i y i -w T x i 2 + λ||w|| 2 2 , (12)
where y i are experimentally measured log K d /K i values, x i are features computed from 3D structures of protein-ligand complexes as described above, and λ is the adjustable regularization parameter that determines the importance of the regularization term with respect to the regression term. Its optimal value is determined using the cross validation procedure. For the kernel regression model, we used a Gaussian kernel to compute the inner product between feature vectors. This allowed us to introduce a non-linear relation between our features and the binding affinity data.

Computational Details

To train our prediction model with LRR and KRR, we used the SHOGUN C++ library [START_REF] Sonnenburg | [END_REF]. We found two adjustable parameters, the regularization parameter λ, and the width of the Gaussian kernel, using a 5-fold cross validation. The confidence interval for the cross-validation was set to 95% and

I c , SASA = 0 I c , SASA > 0 I (A) (B) (C) (D)
Fig. 1 Structural parameters used in this work. (A) A fictitious 2D complex and its interface atoms label highlighted in color. The binding patch on each partner consists of one layer of atoms (I, colored solid balls), as identified by a Voronoi interface model [23,24]. The non interface atoms (I c ) are ascribed to two categories: those which retain solvent accessibility (SASA > 0, dashed balls), and those which do not (SASA = 0, dash-dotted balls) (B) Each interface atom is assigned its shelling order, namely an integer equal to the smallest number of atoms traveled to reach an exposed non interface atom, i.e. an atom belonging to I c and with SASA > 0 (in grey) [31]. (C,D) Atomic packing measured by atomic volumes. The volume of an atom is defined as the volume of the intersection between its ball in the solvent accessible model, and its Voronoi cell [START_REF] Cazals | [END_REF]. Note that this quantity is defined even if the atom retains solvent accessibility. Practically, interfaces and binding patches are computed with Vorshell [24], while atomic surface areas and volumes are computed with sbl-vorlume-pdb.exe [START_REF] Cazals | [END_REF]. Both programs are available from the Structural Bioinformatics Library (SBL), see http://sbl.inria.fr.

the number of runs was set to 3. Atomic typization for the enthalpic features was computed using the Knodle library [36] available at https:// team.inria.fr/nano-d/software/Knodle. Geometric descriptors for the entropic features were computed using the Structural Bioinformatics Library available at http://sbl.inria.fr.

Pose Prediction Protocol

To generate putative binding poses of a ligand with a protein, we used the AutoDock Vina software package with the default scoring function [37]. Using AutoDock Vina, we generated 1,000 docking solutions for the subsequent re-scoring. To start the AutoDock Vina docking, we detected the binding pockets of each protein-ligand complex by analogy with the known homologous complexes. Using PyMOL and the AutoDock Vina plugin available at http://wwwuser.gwdg.de/ ~dseelig/adplugin. html [38], bounding boxes were fitted to the binding pockets. For the MAP4K4 dataset, we progressively tested three bounding boxes of increasing volumes, ranging from 6,125 to 8,685 Å3 for the 4OBO receptor, and ranging from 5,610 to 12,472 Å3 for the 4U44 receptor. In the obtained AutoDock Vina configuration files, the parameter num modes was set to 1,000 and exhaustiveness to 100. We chose all the ligands to be flexible during the docking procedure, and we kept all the protein residues inside the binding pockets rigid. To test the effect of flexible receptor residues, we ran two additional docking simulations for the MAP4K4 dataset. First, we allowed residues Tyr 36, Lys 54, Met 105, Phe 107, and Asp 171 to be flexible. Second, we allowed residues Tyr 36, Met 105, and Phe 107 to be flexible. However, we were not satisfied by the initial sampling results and did not re-score the predictions. We used the default assignment of ligand bond orders predicted by AutoDock Tools. However, one of the conjugated double bonds in the HSP90 set was falsely marked as rotatable. Also, we did not have flexibility in non-aromatic cycles, i.e. in cyclohexyl. To create the initial three-dimensional models of ligands given in the SDF format, we used the Open Babel conversion library [39]. We used the provided templates as the models of protein receptors, 4OBO 4U44 for all the ligands in the MAP4K4 dataset and 2XDX (HSP90 73 and HSP90 179 ligands), 3OW6 (HSP90 40 and HSP90 44 ligands), 4YKR (HSP90 40 and HSP90 44 ligands), and 4YKY (HSP90 164 and HSP90 175 ligands) for the Hsp90 dataset. Then, we re-scored the docking solutions with the Convex-PL scoring function [3], without any further structure optimization, and ranked the poses according to the computed values. If multiple rigid receptor templates or multiple bounding boxes around the binding pockets were used, we re-scored the docking solutions together from all the predictions. For the first stage of the challenge, we submitted five top predictions for each protein-ligand complex for the two datasets. At the second stage of the challenge, we used top-one predictions for the HSP90 dataset for the subsequent re-scoring with multiple scoring functions to predict the binding affinities. At this stage, we did not predict binding poses for the MAP4K4 dataset, since the co-crystal structures were made available by the organizers.

For the assessment of docking poses we used the LigAlign analysis program [40] written for the PyMOL package [41]. We visually inspected the structural alignment of all the ligands in MAP4K and HSP90 datasets and corrected the alignment weight function of LigAlign when needed for the correct alignment.

Results and Discussion

Predicted Docking Poses

Table 1 lists results for our docking predictions for the two sets, MAP4K and HSP90. This table shows the averaged assessment results for the top-1 binding poses, best binding poses, and all the five poses for each protein-ligand complex. As it can be seen from the table, we obtained good predictions for the HSP90 set with the mean RMSD for the best pose of 0.65 Å. These poses along with the crystallographic solutions are shown in Fig. 2. For the MAP4K set, however, AutoDock Vina experienced sampling problems for some of the large ligands. More precisely, in some regions of space no samples were generated while other were sampled densely. As a result, long ligands did not find their place into the corresponding binding pockets. For example, when using 4U44 template to model the receptor, the minimum RMSD in top 1,000 docking poses for MAP12 ligand was 11.2 Å, for MAP13 ligand was 6.5 Å, and for MAP17 ligand was 7.9 Å. Similarly, when using the 4OBO template to model the receptor, the minimum RMSD in top 1,000 docking poses for MAP12 ligand was 5.5 Å, for MAP13 ligand was 5.0 Å, and for MAP15 ligand was 5.0 Å. Conversely, for small ligands such as MAP03, MAP04, MAP21, MAP22, MAP30, and MAP32, we could obtain the minimum RMSD in top AutoDock Vina 1,000 docking solutions within 1 Å at least for one out of two rigid receptor templates. Also, since we did not relax the initial protein structures, some of our binding poses had steric clashes with bulky protein residues, such as phenylalanines. Overall, for the MAP4K set, we achieved the mean RMSD for the best pose of 3.52 Å. Fig. 3 shows some of the successful predictions for this set.

In the MAP4K set predictions there is a big difference of 3.58 Å between RMSD values of mean pose 1 and mean best pose. This is mostly due to the choice of different templates, namely, both 4OBO and 4U44 for all of the predicted ligands. If we discard 4U44-based predictions from the submission, the difference between RMSD values of mean pose 1 and mean best pose reduces to 1.08 Å. Similarly, if we discard 4OBO-based predictions, this difference reduces to 1.71 Å. The difference between RMSD values of mean pose 1 and mean best pose for the HSP90 set is explained by a flipped ring in the HSP90 164 ligand ranked first by the Convex-PL potential. An interesting question is to which extent rescoring with an additional potential is useful for successful pose predictions. To answer this question, we analyzed 1,000 top poses predicted by AutoDock Vina. Table 1 Assessment of the docking pose predictions for the MAP4K and HSP90 datasets. Here, RMSD states for root mean square deviation, RSR states for real space R-factor, and RSCC states for real space correlation coefficient. RMSD is a similarity measure of an atomistic model with the one obtained experimentally. RSR and RSCC are measures of the similarity between an electron-density map calculated directly from a structural model and one calculated from experimental data.

computed with respect to co-crystal structures. We can clearly see that re-scoring with the Convex-PL potential significantly improves the quality of the best pose among top 5 poses. This is true for all AutoDock Vina simulations of this set. However, re-scoring did not improve the mean best pose, neither it improved the mean among top 5 poses. Table 3 shows the analysis of 100 poses predicted by AutoDock Vina for the HSP90 dataset. Here we see that re-scoring with the Convex-PL potential significantly improves all the RMSD values, namely, the mean top pose, the mean best pose, and the mean among all top 5 poses. To conclude, we strongly believe that the default scoring function of AutoDock Vina is too simplistic, since it is mainly designed to describe steric interactions, thus an additional re-scoring step with a more sophisticated scoring function is required. We also think that a subsequent gradient-based optimization step can further improve the quality of predictions.

Another interesting question is how much the quality of pose predictions depends on the choice of the rigid receptor template, the number of flexible residues in the binding pocket, and the binding pocket size. We cannot answer this question generally, since the conclusions should depend on the sampling method. However, we can draw certain conclusions from our docking experiments with AutoDock Vina, which is one of the most popular docking methods. First, if we analyze the effect of the size of the search space in the binding pocket, we can clearly see that the search smaller the volume, the higher the probability to find a better docking pose (see column 4 in Table 2). This emphasizes the insufficient sampling of the binding pocket by AutoDock Vina even with the sampling parameters chosen much higher than the default values (num modes set to 1,000 and exhaustiveness set to 100). Second, flexibility of some of the receptor residues does not, generally, improve the quality of the predictions. This is probably connected to the previous problem, i.e. the insufficiently sampled search space, which is much larger in this case. Finally, the success of the pose prediction critically depends on the chosen receptor template. This dependency cannot be seen clearly in results presented in Table 2 for the MAP4K dataset. However, we analyzed AutoDock Vina predictions for the HSP90 set, where several co-crystal structures are available for the same family of ligands. Table 3 shows RMSD values to co-crystal structures for the HSP90 73 and HSP90 179 ligand predictions if started with 2XDX and 2JJC templates, and also for HSP90 40 and HSP90 44 ligand predictions if started with 4YKR, 3OWD, and 3OW6 templates. We can see very large deviations in the best pose RMSD value depending on the chosen receptor template. Thus, it is very critical to preselect the best receptor model for each of the ligand in a set prior to starting the docking simulations.

Affinity Model Training

As a first affinity prediction computational experiment, we trained our affinity models using the linear and the kernel ridge regressions. Table 4 lists the training results after a five-fold cross validation. We should mention that all the features were scaled to an interval [0; 1] prior to the training. As further computational experiments demonstrated, in order to have an affective regularization, it is sufficient to properly scale the five entropic features, since they have very different amplitudes compared to all the enthalpic features. We also noted that the regularization yielded a dense model, more specifically, the linear model had 5134 / 11765 i.e. 43% of non null terms. This is mainly because some of our 11760 enthalpic features are never met in protein-ligand complexes from the training set, and also because the l-2 norm used in our regression models produces a dense solution. This said, it appears that enthalpy is more important for the prediction of the binding affinities, compared to entropy. We can also see that when a linear model is used it is important to have both b Five flexible residues on the receptor protein were chosen as Tyr 36, Lys 54, Met 105, Phe 107, and Asp 171. c Three flexible residues on the receptor protein were chosen as Tyr 36, Met 105, and Phe 107.

Table 2 Assessment of the docking pose predictions for the MAP4K set when only using the AutoDock Vina scoring function. For each of the two datasets five docking simulations were performed, three with the rigid receptor and bounding boxes around the binding site of a different size, and two with some residues of the receptor inside the binding pocket allowed to be flexible. Only the rigid simulations were used for the subsequent re-scoring. RMSD to the ligands of experimentally obtained structures are computed for the top 5 poses, top 100 poses, and top terms in the approximation of the binding free energy. However, if we use a non-linear model, then the additional entropic term does not contribute to the overall accuracy of the prediction. This can be explained by the fact that the entropic term represents non-additive contributions to the binding free energy and it is very important to have it in the linear model. On the other hand, when we use the kernel ridge regression, the non-additivity of the entropic term is efficiently replaced by the non-linear kernel. Figure 4 shows predicted binding affinities versus experimentally measured values for both linear and non-linear models. We can see both from Table 4 and Fig. 4

Affinity Model Testing

The next computational experiment was to test the trained binding affinity model on the blind data of the MAP4K and HSP90 datasets. Particularly interesting for us was the MAP4K dataset, since the D3R challenge organizers provided 18 co-crystal structures for the prediction of binding affinities. Figure 5 shows the performance of our binding model on the MAP4K dataset. Notably, the accuracy of the model is much lower compared to the training set. We should add that at the first stage of the Challenge for this dataset, we also submitted nine different binding affinity predictions. These were computed by re-scoring with the Convex-PL potential the top binding poses predicted using nine different docking strategies. More precisely, we predicted top-one binding poses with three bounding boxes of different size around the binding site. These were permutated with receptors chosen solely from the 4OBO template, solely from the 4U44 template, and a mix of two. The obtained Spearman ρ coefficients varied from 0.25 to 0.4. Respectively, the Kendal τ coefficients varied from 0.15 to 0.29.

At the second stage of the Challenge, we submitted four different binding affinity predictions. These were computed by re-scoring the co-crystal structures with the following scoring functions : a) the Convex-PL potential (ρ = 0.25, τ = 0.24); b) the entropic term given by Eq. 10 trained with the LRR model without feature rescaling (ρ = 0.44, τ = 0.33, Eq. 11), which provided the best results; c) the free energy approximation given by Eqs. 1,6,10 trained with the LRR model without feature rescaling (ρ = 0.3, τ = 0.19, Eq. 12); d) the free energy approximation given by Eqs. 1,6,10 trained with the KRR model without feature rescaling (ρ = 0.12, τ = 0.07, Eq. 12).

In case of the HSP90 dataset, we obtained very low correlation coefficients. More precisely, we submitted four different binding affinity predictions following the same strategy as in the previous dataset. Re-scoring with the Convex-PL potential resulted in ρ = 0.07 and τ = 0.05, re-scoring with the en-tropic term trained with the LRR model (Eq. 11) without feature rescaling gave ρ = 0.03 and τ = 0.02, re-scoring with the free energy approximation trained with the LRR model (Eq. 12) without feature rescaling gave ρ = 0.07 and τ = 0.05, rescoring with the free energy approximation trained with the KRR model (Eq. 12) without feature rescaling gave ρ = 0.15 and τ = 0.12. The quality of the predictions is very low presumably because most of the co-crystal structures for this dataset were not available. Indeed, in our experiments we used the top-1 putative docking predictions for 180 protein-ligand complexes. Even though we obtained very high quality predictions for six complexes with experimentally solved structures (with a mean RMSD to crystallographic solution of about 1 Å), there is absolutely no guarantee that the same quality of predictions holds for the rest of 174 structures in this set.

To conclude the discussion, we would like to speculate about interpretability of the experimental affinity data. Indeed, IC 50 is a measure of inhibition of a specific process and K i can be used to relate it to binding affinity, but solely for competitive agonists and antagonists. On the other hand, K d is related to the binding energy. Since binding does not necessarily result in inhibition and conversely, the conversion between IC 50 /K i and K d is non-trivial. Therefore, since both the model selection and fitting were performed onlog(K d ), it is not unexpected that predictions of IC 50 /K i are not accurate including the very different scales.

Conclusion

In this challenge we assessed our docking pipeline and tested a new model to predict binding affinity. Our docking predictions appeared to be successful, especially for the HSP90 set. Nonetheless, we experienced certain difficulties in sampling large conformational space for long ligands of the MAP4K set, and this was the main source of inaccuracies in our pose predictions. We critically assessed the quality of pose predictions depending on the choice of rigid receptor templates, the number of flexible residues in the binding pocket, the binding pocket size, and the subsequent re-scoring with the Convex-PL potential. We demonstrated that re-scoring with the Convex-PL potential significantly improves at least the quality of the best pose among the top 5 poses. We also showed that it is very important to choose the right receptor template and to reduce the size of the binding pocket (consequently, the size of the search space) as much as possible. Finally, we showed that flexibility of protein side chains inside the binding pocket does not improve the quality of docking poses. However, our main goal in this challenge was to test the free energy model that consists of the enthalpic term and the entropic term. Free parameters of the model were deduced from the PDBBind database using regression procedures. Even though our model performed very well on the training set, its success rates on the MAP4K and, especially, HSP90 test sets were rather low. The mild performances of our model owe to several facts. First, the very large number of free parameters in the enthalpic term possibly causes overfitting of training data. Second, as in any regression problem, the coverage of the test set by the training set appears critical to ensure performance guarantees-as shown by our post hoc analysis. Given this observation, an interesting research direction to be explored is the application of other regression methods, possibly more tolerant to mild coverages of the parameter space. Finally, it should be kept in mind that binding affinity is a subtle thermodynamic property, calling for an exhaustive analysis of dynamics, including conformational entropy, vibrational entropy, and solvent entropy. While recent work has shown that reliable proxies existed for the latter two terms, accurately accounting for changes in conformational entropy remains an open problem. Filling this gap calls for efficient sampling procedures able to reveal low regions of the energy landscape of the system, together with their occupancy probabilities. As noticed above, such procedures would also be highly valuable to sample large conformational space.

To conclude, this work calls for two complementary extensions aiming at fostering our understanding of the emergence of function (i.e., binding) from a combination of structure and dynamics. The first one will address cases where compact models (such as our entropic model) can be used reliably to estimate affinity from static pauses. The second one will aim at going beyond such cases, so as to obtain equally reliable estimates exploiting conformational ensembles accurately coding more general entropic changes (in particular conformational entropy).

Fig. 2

 2 Fig. 2 Best binding poses for the HSP90 set. Crystallographic structures are shown in blue. Predictions are shown in red. Molecular surface is shown for the crystallographic structure. A) HSP90 40 ligand with the 4YKQ co-crystal structure. B) HSP90 44 ligand with the 4YKT co-crystal structure. C) HSP90 73 ligand with the 4YKW co-crystal structure. D) HSP90 164 ligand with the 4YKX co-crystal structure. E) HSP90 175 ligand with the 4YKZ co-crystal structure. F) HSP90 179 ligand with the 4YKU co-crystal structure.

Fig. 3

 3 Fig. 3 Examples of some of the relatively successful predictions for the MAP4K set. Crystallographic structures are shown in blue. Predictions are shown in red. Molecular surface is shown for the crystallographic structure. A) MAP07 ligand prediction. B) MAP11 ligand prediction. C) MAP14 ligand prediction. D) MAP18 ligand prediction. E) MAP19 ligand prediction. F) MAP23 ligand prediction.

Fig. 4

 4 Fig. 4 Accuracy of the ∆H -T ∆S model on the training set after a five-fold cross validation. Left : training of the linear model. Right : training of the non-linear model.

Fig. 5 Fig. 6

 56 Fig. 5 Accuracy of the ∆H -T ∆S model on the test set MAP4K. Left : training of the linear model. Right : training of the non-linear model.

  Table 2 shows RMSD values of ten docking simulations for the MAP4K dataset

						Obtained results ( overall rank )
	dataset	# submissions	mean	RMSD, Å mean	mean	mean	RSR mean	mean	mean	RSCC mean	mean
			pose 1	best pose	all poses	pose 1	best pose	all poses	pose 1	best pose	all poses
	MAP4K	35	7.1 (28)	3.52 (12)	6.54 (25)	3.39 (7)	2.72 (2)	3.43 (8)	0.37 (20-29)	0.15 (27-35)	0.35 (21-31)
	HSP90	52	1.07 (8)	0.65 (5)	1.31 (5)	3.57 (5)	3.12 (3)	3.8 (3)	0.71 (8-9)	0.57 (7-9)	0.7 (5-6)

Table 3

 3 Assessment of the docking pose predictions for the HSP90 set when only using the AutoDock Vina scoring function. RMSD to the ligands of experimentally obtained structures are computed for the top 5 poses, top 10 poses, and top 100 poses. See main text for details.

	1,000

Table 4

 4 Affinity model training accuracy measured as a root mean squared error.

User manual: http://sbl.inria.fr/doc/Space_filling_ model_shelling_diagram_surface_encoding-user-manual.html

User manual: http://sbl.inria.fr/doc/Space_filling_ model_surface_volume-user-manual.html

Acknowledgments

The authors thank Dr. Petr Popov from MIPT Moscow for the initial analysis of the HSP90 targets.


Brice Arnaud

Sébastien Lebègue

Guillaume Raffy

Final draft post-refeering (pre-proofs

Keywords: numbers: 63, 20, dk, 65, 40, De

come

Anisotropic thermal expansion of bismuth from first principles

I. INTRODUCTION

The semimetal bismuth is of interest both scientifically and technologically. Indeed, it exhibits many fascinating properties, like giant magnetoresistance 1,2 , thermoelectricity 3,4 , large diamagnetism 5 that can be ascribed to the peculiar electronic structure of bismuth, namely the small overlap between the valence and the conduction bands giving rise to a Fermi surface made of tiny electron and hole pockets 6 . Real-world applications of bismuth related to the aforementioned properties range from hall magnetometry 7 to diamagnetic levitation on the microscale 8 .

The possibility to drive bismuth strongly out of equilibrium by an ultrashort laser pulse is also behind a huge amount of experimental [9][10][11][12][13][14] and theoretical work [15][16][17][18][19] . From a theoretical point of view, the ultrafast dynamics of coherent optical phonons has been tackled by means of first-principles calculations where the lattice parameters are kept constant. However, the development of strain from coherent acoustic phonons on a picosecond timescale is still poorly understood 20 and has never been addressed by ab-initio calculations. A prerequisite for achieving such a goal is to demonstrate that the thermal expansion of bismuth at equilibrium can be understood and predicted by performing ab-initio calculations. A good strategy is to resort to the quasiharmonic approximation, where the atoms of the crystal are considered to undergo harmonic oscillations, but with frequencies that depend on strain. This approximation, when combined with density functional perturbation theory 21 , has been found to produce thermal expansion coefficients in good agreement with experimental results well below the melting temperature of isotropic [22][23][24] and anisotropic solids [25][26][27][28] .

The paper is organized as follows. In section II, we give an account of the technicalities used to perform our first-principles calculations. In section III, we describe the crystallographic structure of bismuth and compare our calculated lattice constants with and without spin-orbit interaction (SOI) to the experimental lattice constants at 4 K obtained from X-ray measurements. In section IV, we explain how the elastic constants at 0 K can be computed using a finite strain method and make a comparison with available experimental results indirectly obtained by measuring the sound wave velocities for different directions and polarizations. We also discuss the impact of SOI on the calculated elastic constants. In section V, we introduce the theory allowing the calculation of the thermal expansion tensor of bismuth within the quasiharmonic approximation and compare the thermal expansion co-Final draft post-refeering (pre-proofs) efficients parallel and perpendicular to the trigonal axis to some measurements made using an optical lever dilatometer. We also unravel the role respectively played by the elasticity and the anharmonicity in the anisotropy of the thermal expansion coefficients of bismuth.

Finally, the specific heat at constant pressure is reported and compared with calorimetry measurements.

II. COMPUTATIONAL DETAILS

All the calculations are performed using the ABINIT code 29 . We use a plane-wave basis set, the Hartwigsen-Goedecker-Hutter (HGH) pseudopotentials 30 and the local density approximation (LDA) for the exchange-correlation functional. We carefully check the convergence of our results with respect to the wave function cut-off and the k-point sampling of the Brillouin zone. A 40 Ry cut-off and a 16 × 16 × 16 mesh for Brillouin zone sampling ensure that our results (lattice parameters, elastic constants, phonon frequencies) are well converged. The dynamical matrix is explicitly calculated on a 8 × 8 × 8 q-point mesh using density functional perturbation theory 31 and the phonon frequencies are Fourier interpolated on a 32 × 32 × 32 q-point mesh in order to compute the thermal expansion tensor. We include SOI in all calculations but also present the lattice parameters and elastic constants obtained without SOI in order to highlight the crucial role played by SOI.

Final draft post-refeering (pre-proofs)

III. LATTICE PARAMETERS AT ZERO TEMPERATURE

Bismuth crystallizes in a rhomboedral structure, also called A7 structure, with two atoms per unit cell. The vectors spanning the unit cell are given by

a 1 = aξ, - aξ √ 3 , h ; a 2 = 0, 2aξ √ 3 , h ; a 3 = -aξ, - aξ √ 3 , h , (1)
where ξ = sin[ α 2 ] and h = a 1 -4 3 ξ 2 . The length of the three lattice vectors is equal to a and the angle between any pair of vector is α. The two atoms belonging to the unit cell are located at ±u (a 1 + a 2 + a 3 ) where u is a dimensionless parameter and a 1 +a 2 +a 3 is parallel to the ternary axis (C 3 axis). Alternatively, the structure can be viewed as an hexagonal structure spanned by the following three lattice vectors ã1 = a 1a

2 ; ã2 = a 2 -a 3 ; ã3 = a 1 + a 2 + a 3 , (2)
where ã1 = ã2 ≡ a ⊥ and ã3 ≡ a . The lattice cell parameters of the two structures are related to each other by the following relations

a ⊥ = 2a sin α 2 a = 1 3 3a 2 ⊥ + a 2 a = a 3 + 6 cos(α) sin α 2 = 3 2 a ⊥ / 3a 2 ⊥ + a 2 (3)
All the calculations have been performed using the rhombohedral structure because it contains three times atoms less than the hexagonal structure. However, the hexagonal structure, as will be seen later, is more convenient to define thermal expansion coefficients. Our calculated LDA lattice parameters with and without SOI are given in Table I along with the experimental results at 4.2 K 32 . The agreement between theory and experiment is significantly improved when SOI is included 33 . Indeed, a ,0 and a ⊥,0 are respectively underestimated from 0.69 % (1.58 %) and 0.26 % (1.27 %) with respect to experiments, leading to an underestimation of the equilibrium volume V 0 of 1.2 % (4.1 %) when SOI is included (neglected). Thus, the inclusion of SOI is mandatory to achieve a better description of the equilibrium lattice parameters of bismuth.

Final draft post-refeering (pre-proofs)

IV. ELASTIC CONSTANTS AT ZERO TEMPERATURE

The elastic properties of a bismuth crystal can be inferred from the theory of elasticity.

The Lagrangian strain tensor η is defined as

η ab = ǫ ab + 1 2 k ǫ ak ǫ kb ( 4
)
where ǫ is the linear strain tensor which transforms a vector a into (1 + ǫ)a. The energy of the crystal per unit cell E can be expanded in power series with respect to the strain η as

E[η] = E 0 + V 0 2 ijkl C ijkl η ij η kl + • • • (5)
where E 0 and V 0 are the energy and the volume of the unstrained unit cell and C ijkl are the elastic stiffness constants of the crystal. Using Voigt's notation, Eq. 5 can be written as

E[η] = E 0 + V 0 2 αβ C αβ η α η β + • • • (6)
where the fourth-rank stiffness tensor has been replaced by a 6 × 6 matrix C. By virtue of the rhombohedral structure A7 of bismuth (space group R3m), the matrix C can be cast in the form

C =                  C 11 C 12 C 13 C 14 0 0 C 12 C 11 C 13 -C 14 0 0 C 13 C 13 C 33 0 0 0 C 14 -C 14 0 C 44 0 0 0 0 0 0 C 44 C 14 0 0 0 0 C 14 C 66                  (7)
provided that the z axis is taken along the trigonal axis. Here, only 6 elements are independent and

C 66 = 1 2 [C 11 -C 12 ]
. In order to compute these elements, we consider six sets of deformations parametrized by η η 1 = (η, η, η, 0, 0, 0) ; η 2 = (η, 0, 0, 0, 0, 0) ; η 3 = (0, 0, η, 0, 0, 0) η 4 = (0, 0, 0, 2η, 0, 0) ; η 5 = (η, η, 0, 0, 0, 0) ; η 6 = (η, 0, 0, 2η, 0, 0), (8) where η is varied between -0.01 and 0.01 with a step of 0.001. For each deformation labelled by i and each value of η, we build the Lagrangian strain matrix η and solve Eq. 4 in an iterative way to obtain the matrix ǫ. Thus, we generate a distorted cell from the undistorted

Final draft post-refeering (pre-proofs)

one by using the matrix ǫ and we compute the total energy E i (η) ≡ E[η i ] of the distorted structure where the atomic positions are fully relaxed. Then, the energy per unit volume is fitted by a polynomial of order 4

E i (η) V 0 = 4 j=0 A i j η j , (9)
where A i 0 = E 0 /V 0 , A i 1 = 0 and A i 2 can be expressed as a function of the unknown second order elastic constants as

A 1 2 = C 11 + C 12 + 2C 13 + 1 2 C 33 ; A 2 2 = 1 2 C 11 ; A 3 2 = 1 2 C 33 A 4 2 = 2C 44 ; A 5 2 = C 11 + C 12 ; A 6 2 = 1 2 C 11 + 2C 14 + 2C 44 (10)
by using Eq. 6 and 7. The computed total energy per unit cell at 0 K including SOI (circles)

together with the polynomial fit (full lines) are displayed in Fig. 1 for the two deformations respectively labelled η 3 (a varies while keeping a ⊥ constant) and η 5 (a ⊥ varies while keeping a constant). The elastic constant C 33 /2 (see Fig. the trigonal axis than against a strain perpendicular to it. We do not show the energy per unit cell for the four remaining deformations but compare in Table II our calculated elastic constants with and without SOI to the experimental elastic constants indirectly obtained from the measurements of ultrasonic wave velocities by the pulse echo technique 34 . The calculated elastic constants without SOI are all overestimated with respect to experiments.

Taking into account SOI leads to a decrease of all elastic constants that can be explained in a very qualitative way as follows: SOI mixes bonding and antibonding states and not only increases the equilibrium volume by 2.9 % (see Table I) but also softens the elastic constants by as much as 45 % (the larger effect being for C 44 ). As shown in Table II, the overall agreement between theory and experiment is significantly improved when SOI is included. The only exceptions are the C 14 and C 44 elastic constants which are underestimated by about 30% with respect to experiment. Such a discrepancy might be ascribed to experimental uncertainties associated with the sample and its bonded transducer and/or to long range effects, like van der Waals interactions, not captured by the LDA exchange-correlation functional 35 . It is also worth mentioning that all the calculated elastic constants with or without SOI satisfy Born's mechanical stability for a rhomboehedral structure 36,37 ensuring that bismuth is stable at 0 K.

Inverting the matrix C defined in Eq. 7 leads to the following expression for the compliance matrix :

S =                 
S 11 S 12 S 13 S 14 0 0 S 12 S 11 S 13 -S 14 0 0 S 13 S 13 S 33 0 0 0 S 14 -S 14 0 S 44 0 0

0 0 0 0 S 44 2S 14 0 0 0 0 2S 14 S 66                  (11)
Final draft post-refeering (pre-proofs) where S 66 = 2(S 11 -S 12 ). The anisotropy in the elastic properties of Bi can be ascertained by introducing a parallel and a perpendicular compressibility respectively defined as

χ = - ∂ǫ 33 ∂P = 2S 13 + S 33 = (-2C 13 + C 11 + C 12 )/M (12)
and

χ ⊥ = - ∂ǫ 11 ∂P = - ∂ǫ 22 ∂P = S 11 + S 12 + S 13 = (C 33 -C 13 )/M, (13)
where

M = C 33 (C 11 + C 12 ) -2C 2 13
. The computed values of χ and χ ⊥ with and without SOI are reported in Table III. The compressibilities are underestimated with respect to experiments when SOI is neglected. Taking into account SOI leads to a very good agreement between the theoretical and experimental compressibilities in accordance with the fact that the elastic constants that play a role in the compressibilities are much better described when SOI is included (see Table II). The ratio of linear compressibilities χ /χ ⊥ calculated when SOI is included indicates that the contraction along the trigonal axis is about 2.7 times larger than the contraction perpendicular to it upon applying an hydrostatic pressure. As stated before and illustrated in Fig. 1, Bi is stiffer perpendicularly to the trigonal axis than parallel to it. The bulk modulus B that measures material's resistance to uniform compression is defined as B = -V 0 ∂P ∂V = 1/χ where the hydrostatic compressibility χ is defined as χ = χ + 2χ ⊥ . Using the values of χ and χ ⊥ reported in Table III, we obtain a theoretical value of 33.65 (39.15) GPa for the bulk modulus B when the SOI is included (neglected). Thus, the theoretical value including the SOI agrees well with the experimental value of 34.23 GPa at 4 K 34 . All the forthcoming calculations include SOI.

Final draft post-refeering (pre-proofs)

V. THERMAL EXPANSION

We now present an analysis of the thermal expansion of bismuth using Grüneisen's theory.

We follow the approach of Schelling and Keblinski 26 and emphasize the differences in the formalism used to treat the thermal expansion in anisotropic and cubic solids. Since bismuth is an anisotropic solid, the thermal expansion can be described in terms of the 3 × 3 thermal expansion tensor α whose components are written in terms of the strain tensor ǫ as

α ab = ∂ǫ ab ∂T σ ( 14
)
where T is the temperature and where the subscript σ means that the temperature derivative is taken at constant stress. According to the generalized form of Hookes' law, we have:

ǫ ab = de S abde σ de , (15)
where σ de and S abde are the stress and the compliance elastic tensor respectively. Eq. 15 can be inverted and leads to

σ de = ij C deij ǫ ij , (16)
where C deij is the stiffness elastic tensor. Thus, we have de

S abde C deij = δ ai δ bj (17)
By differentiating Eq. 17 with respect to T and by using Eq. 15 and 16, it is straightforward to show that the components of α can be written as

α ab = - de ∂ǫ ab ∂σ de T × ∂σ de ∂T ǫ , (18)
where the temperature derivative of the stress tensor is taken at constant strain. This equation generalizes the expression

α = 1 3B ∂P ∂T V (19)
defining the thermal expansion coefficient α of a cubic solid as a function of the bulk modulus B and the partial derivative of the pressure with respect to temperature at constant volume.

The stress tensor σ de appearing in Eq. 18 is defined as

σ de = 1 V 0 ∂F ∂ǫ de ( 20
)
where F is the Helmholtz free energy per unit cell of the crystal defined as

F [ǫ] = E[ǫ] + F vib [ǫ, T ] = E[ǫ] + 1 N q,λ hω λ (q) 2 + k B T 1 N q,λ ln 1 -exp - hω λ (q) k B T (21)
where the electron entropy is discarded and the vibrational contribution F vib [ǫ, T ] is computed within the harmonic approximation. Here, ω λ (q) is the frequency of the phonon mode (q, λ) corresponding to wavevector q and polarization λ and N is the number of qpoints included in the summation. Hence, combining Eq. 20 and Eq. 21 leads to

σ de = 1 V 0   ∂E ∂ǫ de - 1 N q,λ γ de q,λ hω λ (q) n q,λ + 1 2   , (22)
where n q,λ is the Bose occupation factor at temperature T for a phonon with frequency ω λ (q) and γ de q,λ is a generalized mode Grüneisen parameter given by

γ de q,λ = - ∂ ln ω λ (q) ∂ǫ de (23)
Note that σ de is temperature dependent because of the second term in Eq. 22 and that σ de (T → 0) is renormalized by zero point atomic motions. By derivating Eq. 22 with respect to temperature T at constant strain ǫ, we get

∂σ de ∂T ǫ = - q,λ γ de q,λ C q,λ (24)
where

C q,λ = 1 N hω λ (q) V 0 ∂n q,λ ∂T ǫ ( 25
)
is the contribution of mode (q, λ) to the lattice specific heat per unit volume at constant volume C V . Thus C V is given by

C V (T ) = q,λ C q,λ = k B V 0 1 N q,λ hω λ (q) 2k B T 2 1 sinh 2 hω λ (q) 2k B T . ( 26
)
Finally, by using Eq. 15 and inserting Eq. 24 in Eq. 18, we obtain

α ab = q,λ C q,λ de S abde γ de q,λ . (27)
In a completely harmonic lattice, the frequencies would be independent of the strain and the γ de q,λ would be zero, leading to a zero thermal expansion. Eq. 27 generalizes the expression

α = 1 3B q,λ C q,λ γ q,λ , (28)
Final draft post-refeering (pre-proofs)

giving the thermal expansion coefficient α for a cubic solid. Here C q,λ is defined in Eq. 25 and γ q,λ is the mode Grüneisen parameter defined as

γ q,λ = -∂ ln[ω λ (q)]/∂ ln V. (29)
Usually, the phonon frequencies ω λ (q) decrease as the volume V increases giving rise to positive Grüneisen parameters and thus to a positive thermal expansion coefficient α at any temperature as can be inferred from Eq. 28. Inserting the non-zero matrix elements of S allowed by symmetry (see Eq. 11) in Eq. 27 and taking into account the fact that γ de q,λ = 0 when d = e leads to the following expressions

α ⊥ ≡ α 11 ≡ α 22 = q,λ C q,λ (S 11 + S 12 ) γ ⊥ q,λ + S 13 γ q,λ (30)
and

α ≡ α 33 = q,λ C q,λ 2S 13 γ ⊥ q,λ + S 33 γ q,λ (31)
where α ⊥ and α are the thermal expansion coefficient respectively within the basal plane and along the ternary axis. Here,

γ ⊥ q,λ ≡ 1 2 γ 11 q,λ + γ 22 q,λ = - a ⊥,0 2ω 0 λ (q) ∂ω λ (q) ∂a ⊥ (32)
and γ q,λ ≡ γ 33 q,λ = -

a ,0 ω 0 λ (q) ∂ω λ (q) ∂a ( 33
)
where a ⊥,0 and a ,0 are the LDA lattice parameters reported in Table I and ω 0 λ (q) are the phonon frequencies calculated for these lattice parameters. We do not compare our computed ω 0 λ (q) with available experimental data because the agreement between theory and experiments has already been highlighted in a thorough study based on calculations performed using the ABINIT code 33 . A finite difference scheme based on a relative variation of ±0.2 % of a ⊥ and a around a ⊥,0 and a ,0 is used to compute the partial derivatives of the phonon frequencies ω λ (q) with respect to a ⊥ and a . Hence, the mode Grüneisen parameters respectively defined in Eq. 32 and Eq. 33 are computed for a 32 × 32 × 32 qpoints grid and for all polarizations λ. The Figure . 2 shows that the mode Grüneisen parameters γ q,λ and γ ⊥ q,λ are rather scattered for the acoustic modes with quite large positive values but also negative values. About 40 % and 18 % of the mode Grüneisen parameters are negative for the first and second Final draft post-refeering (pre-proofs) transverse acoustic branch (TA) while less than 5 % of the mode Grüneisen parameters are negative for the longitunal acoustic branch (LA).

By comparing Figs 2(a) and 2(b), we note that the mode Grüneisen parameters γ q,λ are slightly larger than the mode Grüneisen parameters γ ⊥ q,λ for acoustic modes with energy ranging from 0 to 7.5 meV (TA+LA). On the contrary, the mode Grüneisen parameters γ q,λ are almost 3 times smaller in average than the mode Grüneisen parameters γ ⊥ q,λ for optical modes with energy ranging from 8 to 13.5 meV (TO+LO). In other words, the optical phonon frequencies are more sensitive to a variation of a ⊥ than to a variation of a while the opposite is true for the acoustic phonon frequencies.

We can also introduce macroscopic Grüneisen functions

γ ⊥, =   q,λ γ ⊥, q,λ C q,λ   /C V (34)
where C q,λ and C V are respectively defined in Eq. 25 and 26. The calculated lattice specific heat at constant volume C V displayed in Fig. 3(a) is in good agreement with the experimental lattice specific heat at constant pressure up to the Debye temperature θ D = 119 K 41 .

As shown in Fig. 3(b), the behaviour of the temperature dependent Grüneisen functions γ ⊥ and γ is quite complex with a crossover around 40 K. However, γ ⊥ and γ saturate towards γ ⊥ ∞ = 1.32 and γ ∞ = 1.11 when T ≫ θ D in accordance with the fact that lim

T →∞ γ ⊥, (T ) = 1 6N q,λ γ ⊥, q,λ , (35)
since C q,λ (T ) → k B /V 0 N and C V (T ) → 6k B /V 0 when T → ∞. Interestingly, the high temperature limiting values of γ ⊥ ∞ and γ ∞ extracted from experimental results 38 are estimated to be 1.32 and 1.10 and are in excellent agreement with our calculated values.

Following the approach of Munn 40 , the principal coefficients of thermal expansion defined Final draft post-refeering (pre-proofs) in Eq. 30 and 31 can also be expressed as

α ⊥ = C V χ ⊥ γ ⊥ + S 13 γ -γ ⊥ (36) α = C V χ γ -2S 13 γ -γ ⊥ (37)
where the directional Grüneisen functions γ ⊥, are defined in Eq. 34 and the compressibilities χ and χ ⊥ are respectively defined in Eq. 12 and 13. It is a good first approximation to treat both the compressibilities and the compliance matrix element S 13 as constant, and regard the temperature dependence of the coefficients of thermal expansion as due solely to that of the heat capacity C V and the Grüneisen functions γ ⊥, . Such expressions for the coefficients of thermal expansion allow to disentangle the role of the anisotropy in either the Grüneisen functions or the elastic constants. Fig. 3(b) shows that γ = γ ⊥ for T=40 K. Thus, the anisotropy in the thermal expansion coefficients measured by α /α ⊥ is given by χ /χ ⊥ = 2.71 (see Table III) and is only due to the anisotropy in the elastic properties.

When moving away from the crossover temperature (T=40 K), the term proportionnal to γγ ⊥ starts to play a role as the cross-compliance S 13 = -C 13 /M (∼ -0.94 MBar -1 ) has the same order of magnitude as χ and χ ⊥ (see Table III). For α , the correction arising from the anisotropy in the Grüneisen functions is given by -2C V S 13 γγ ⊥ and remains very small at low temperature since C V → 0 when T → 0. Thus, the correction is small and positive for T < 40 K since S 13 < 0 and γγ ⊥ > 0 in this low temperature regime.

However, the correction becomes non negligible at higher temperatures and negative as the sign of γγ ⊥ changes when T > 40 K. The high temperature limit of this correction is given by -12k B S 13 (γ ∞ -γ ⊥ ∞ )/V 0 and is sketched as a vertical downward arrow in Fig. 4. For α ⊥ , the same type of conclusion holds but the sign and the magnitude of the correction is changed since it is given by C V S 13 γγ ⊥ . The high temperature limit of this correction is given by 6k B S 13 (γ ∞ -γ ⊥ ∞ )/V 0 and is depicted as a vertical upward arrow in Fig. 4 since this quantity is positive. The thermal expansion coefficients parallel (α ) and perpendicular (α ⊥ ) to the ternary axis calculated neglecting (dashed curves) and including (full curves) the anisotropy in the Grüneisen functions are displayed in Fig. 4. At low temperature (T < 40 K), the anisotropy in the Grüneisen functions plays a minor role since the full curves almost coincide with the dashed curves. However, it starts to contribute at higher temperature since it reduces α and increases α ⊥ significantly, thereby reducing the anisotropy in the thermal expansion coefficients given by α /α ⊥ and bringing our calculated Final draft post-refeering (pre-proofs) 39 . The full curves are calculated using Eq. 36 and 37 while the dashed curves are calculated by neglecting the anisotropy in the Grüneisen functions, which is tantamount to put S 13 = 0 in Eq. 36 and 37.

thermal expansion coefficients in closer agreement with the measurements made using an optical lever dilatometer 38,39 . Thus, the quasi-harmonic approximation based on quantities calculated at 0 K provides quite accurate results compared to experimental data 38,39 , even near the melting temperature of bismuth where renormalization effects due to the temperature dependence of elastic constants might be important 24 and where the applicability of the quasiharmonic approximation is also questionable because large anharmonic effects are expected. It's worth coming back to the calculated lattice specific heat at constant volume C V shown in Fig. 3(a). Above θ D , the experimental lattice specific heat at constant pressure deviates from C V and should in principle be compared to C p . A very simple reasoning based on thermodynamics [START_REF] Callen | Thermodynamic[END_REF] shows that

C p -C V = α 2 BT, (38)
where α ≡ 1

VI. CONCLUSION

We performed first-principles calculations in order to understand the anisotropic thermal expansion of bismuth. First, we computed the elastic constants of bismuth at 0 K using a finite strain approach. All the elastic constants, with the exception of the C 14 and C 44 elastic constants, are found to be in good agreement with experimental results 34 when the SOI is included. We also calculated the hydrostatic compressibilities along the ternary axis (χ ) and perpendicular to it (χ ⊥ ) and found that the anisotropy in the directional compressibilities is large since χ /χ ⊥ ∼ 2.7. Then, we computed the thermal expansion coefficients parallel (α ) and perpendicular (α ⊥ ) to the ternary axis using the quasiharmonic approximation. These quantities are found to be in close agreement with experiments 38,39 . Another outcome of our calculations is that the anisotropy in the thermal expansion coefficients is essentially governed by the anisotropy in the mechanical properties below the Debye temperature θ D while both the anisotropy in the directional compressibilities and in the directional Grüneisen functions play a role at higher temperatures. Finally, this work is a first step towards a first-principles description of the thermal/non thermal expansion in laser-excited bismuth 20 , where the electron system is not equilibrated with the phonon system 17,18 .

FIG. 1 .

 1 FIG.1. Energy per unit cell (in meV) as a function of η (dimensionless quantity) for the deformations respectively labelled η 3 (Left panel) and η 5 (Right panel). The circles denote the results of the LDA calculations including SOI whereas the full lines are the result of the polynomial fit according to Eq. 9. The zero of energy corresponds to the energy of the fully relaxed structure (unstrained reference structure of volume V 0 ).

FIG. 2 .

 2 FIG.2. Calculated Grüneisen parameters γ q,λ and γ ⊥ q,λ as a function of phonon energy ω in meV (Left vertical scale) together with the phonon density of states (Right vertical scale).

FIG. 3 .

 3 FIG. 3. (a) Calculated lattice specific heat per unit volume at constant volume C V (solide curve) or at constant pressure C p (dashed curve) compared to experimental data (open circles) from Ref. 41 for temperatures T up to the melting temperature of 545 K. (b) Grüneisen functions γ ⊥, as a function of temperature T up to the Debye temperature θ D of 119 K 41 .

FIG. 4 .

 4 FIG.4. Coefficient of linear thermal expansion of bismuth within the basal plane (α ⊥ ) and along the ternary axis (α ) as a function of temperature T up to the melting temperature. Experimental data are denoted by circles38 and losanges39 . The full curves are calculated using Eq. 36 and 37

=

  α + 2α ⊥ is the volumetric expansion coefficient and B is the temperature dependent bulk modulus, which can be approximated by it's 0 K value. The Final draft post-refeering (pre-proofs) calculated C p displayed in Fig.3(a) coincides with the calculated C V below the Debye temperature θ D because α(T ) → 0 when T ≪ θ D , making the correction α 2 BT very small in this temperature range. When T > θ D , α(T ) becomes constant and the correction increases linearly with T . The agreement between theory and experiment is improved since all the experimental data (open circles) collapse on the calculated C p (dashed curve) up to room temperature.

TABLE I .

 I Calculated LDA lattice parameters with and without SOI compared to the experimental results of Ref 32 .

		Rhombohedral structure Hexagonal structure
		a 0 ( Å) α 0 ( o ) V 0 ( Å3 ) a ,0 ( Å) a ⊥,0 ( Å)
	Experiment	4.724 57.35	69.97	11.796	4.533
	Theory (without SOI) 4.653 57.48	67.12	11.610	4.475
	Theory (with SOI)	4.697 57.53	69.10	11.714	4.521

TABLE II .

 II Calculated LDA elastic constants with and without SOI of Bi compared to the experimental results at 4.2 K of Ref.34 .C 11 (GPa) C 12 (GPa) C 13 (GPa) C 14 (GPa) C 33 (GPa) C 44 (GPa)

	Experiment 34	69.3	24.5	25.4	8.40	40.4	13.5
	Present work (without SOI)	84.6	30.2	27.9	9.8	46.1	16.0
	Present work (with SOI)	67.7	25.0	24.3	5.9	40.6	8.7

TABLE III .

 III Hydrostatic compressibilities χ, χ and χ ⊥ of Bi (in MBar -1 ) obtained by using the elastic constants computed with and without SOI (see TableII) compared to the values inferred from the experimental results at 4.2 K of Ref.34 . The ratio χ /χ ⊥ is also given to quantify the anisotropy of the elastic properties of Bi.

		χ χ χ ⊥ χ /χ ⊥
	Experiment 34	2.92 1.72 0.60 2.86
	Present work (without SOI) 2.55 1.58 0.49 3.24
	Present work (with SOI) 2.97 1.71 0.63 2.71

P. B. Alers and R. T. Webber, Phys. Rev. 91, 1060 (1953).

ACKNOWLEDGMENTS

Calculations were performed using HPC resources from GENCI-CINES (project 095096).


Paul Feliot

Julien Bect

Emmanuel Vazquez

Bmoo

BMOO: a Bayesian Multi-Objective Optimization algorithm

come

We address the problem of derivative-free multi-objective optimization of real-valued functions subject to multiple inequality constraints. The problem consists in finding an approximation of the set

Γ = {x ∈ X : c(x) ≤ 0 and x ∈ X such that f (x ) ≺ f (x)} (1)
where

X ⊂ R d is the search domain, c = (c i ) 1≤i≤q is a vector of constraint functions (c i : X → R), c (x)
≤ 0 means that c i (x) ≤ 0 for all 1 ≤ i ≤ q, f = (f j ) 1≤j≤p is a vector of objective functions to be minimized (f j : X → R), and ≺ denotes the Pareto domination rule. Both the objective functions f j and the constraint functions c i are assumed to be continuous. The search domain X is assumed to be compact-typically, X is a hyper-rectangle defined by bound constraints. Moreover, the objective and constraint functions are regarded as black boxes and, in particular, we assume that no gradient information is available. Finally, the objective and the constraint functions are assumed to be expensive to evaluate, which arises for instance when the values f (x) and c(x), for a given x ∈ X, correspond to the outputs of a computationally expensive computer program. In this setting, the emphasis is on building optimization algorithms that perform well under a very limited budget of evaluations (e.g., a few hundred evaluations).

We adopt a Bayesian approach to this optimization problem. The essence of Bayesian optimization is to choose a prior model for the expensive-to-evaluate function(s) involved in the optimization problem-usually a Gaussian process model [START_REF] Santner | The design and analysis of computer experiments[END_REF][START_REF] Williams | Gaussian processes for machine learning[END_REF] for tractability-and then to select the evaluation points sequentially in order to obtain a small average error between the approximation obtained by the optimization algorithm and the optimal solution, under the selected prior.

Our contribution is twofold. First, we propose a new sampling criterion, which can handle both multiple objectives and multiple non-linear constraints. This criterion corresponds to a one-step look-ahead Bayesian strategy, using the dominated hyper-volume as a utility function (following in this respect [START_REF] Emmerich | Single-and multi-objective evolutionary optimization assisted by Gaussian random field metamodels[END_REF]). The dominated hyper-volume is defined using an extended domination rule, which handles both the objectives and constraints jointly, in the spirit of [START_REF] Fonseca | Multiobjective optimization and multiple constraint handling with evolutionary algorithms. I. A unified formulation[END_REF]. This makes it possible to deal with problems where no feasible solution is known at the beginning of the optimization procedure. The sampling criterion may be written under the form of an integral:

ρ n (x) = Gn P n (ξ(x) y) dy , x ∈ X (2)
where G n is the set of all non-dominated points, P n stands for the posterior probability given n evaluations, ξ is the vector-valued random process modeling the objectives and the constraints, and denotes our extended domination rule. Several criteria from the literature are recovered as special cases.

The second part of our contribution resides in the numerical methods employed to compute and optimize the sampling criterion. Indeed, this criterion takes the form of an integral over the

BMOO: a

 a Bayesian Multi-Objective Optimization algorithm P. Feliot IRT SystemX / Laboratoire des Signaux et Systèmes Supervisor(s): J. Bect, E. Vazquez Ph.D. expected duration: 2014-2017 Adress: IRT SystemX, 8 avenue de la Vauve, 91120 Palaiseau Email: paul.feliot@irt-systemx.fr Abstract:

Acknowledgements. This research work has been carried out in the frame of the Technological Research Institute SystemX, and therefore granted with public funds within the scope of the French Program Investissements d'Avenir. Short biography -I am titular of an engineering diploma from the french engineering school ISAE-Supaero (Toulouse) and of a Master of Research from the Paul Sabatier University (Toulouse). I started my PhD thesis with Supelec at the Technological Research Institute SystemX in January 2014 under the ROM project (model Reduction and Multi-physic Optimization). Industrial partners of the project are SAFRAN (Snecma), Cenaero, Airbus group, Renault and ESI group.

space of constraints and objectives, for which no analytical expression is available in the general case. Besides, it must be optimized at each iteration of the algorithm to determine the next evaluation point. In order to compute the integral, we use an algorithm similar to the subset simulation method [START_REF] Au | Estimation of small failure probabilities in high dimensions by subset simulation[END_REF][START_REF] Cérou | Sequential Monte Carlo for rare event estimation[END_REF], which is well known in the field of structural reliability and rare event estimation. For the optimization of the criterion, we resort to Sequential Monte Carlo (SMC) techniques, following earlier work by [START_REF] Benassi | Bayesian optimization using sequential Monte Carlo[END_REF] for single-objective bound-constrained problems. The resulting algorithm is called BMOO (for Bayesian multi-objective optimization).

The new algorithm is compared with reference methods on both single-and multi-objective optimization problems with very promising results.