Dzmitry Padhorny

Andrey Kazennov

Brandon S Zerbe

Kathryn A Porter

Bing Xia

Scott E Mottarella

Yaroslav Kholodov

David W Ritchie

Sandor Vajda

Dima Kozakov

Michael Levitt

Protein-protein docking by fast generalized Fourier transforms on 5D rotational manifolds

Energy evaluation using fast Fourier transforms (FFTs) enables sampling billions of putative complex structures and hence revolutionized rigid protein-protein docking. However, in current methods, efficient acceleration is achieved only in either the translational or the rotational subspace. Developing an efficient and accurate docking method that expands FFT-based sampling to five rotational coordinates is an extensively studied but still unsolved problem. The algorithm presented here retains the accuracy of earlier methods but yields at least 10-fold speedup. The improvement is due to two innovations. First, the search space is treated as the product manifold SO(3) × (SO(3)∖S 1 ), where SO(3) is the rotation group representing the space of the rotating ligand, and (SO(3)∖S 1 ) is the space spanned by the two Euler angles that define the orientation of the vector from the center of the fixed receptor toward the center of the ligand. This representation enables the use of efficient FFT methods developed for SO [START_REF] Vajda | Convergence and combination of methods in proteinprotein docking[END_REF]. Second, we select the centers of highly populated clusters of docked structures, rather than the lowest energy conformations, as predictions of the complex, and hence there is no need for very high accuracy in energy evaluation. Therefore, it is sufficient to use a limited number of spherical basis functions in the Fourier space, which increases the efficiency of sampling while retaining the accuracy of docking results. A major advantage of the method is that, in contrast to classical approaches, increasing the number of correlation function terms is computationally inexpensive, which enables using complex energy functions for scoring.

protein docking | manifold | FFT D etermining putative protein-protein interactions using genome-wide proteomics studies is a major step toward elucidating the molecular basis of cellular functions. Understanding the atomic details of these interactions, however, requires further biochemical and structural information. Although the most complete structural characterization is provided by X-ray crystallography, solving the structures of protein-protein complexes is frequently very difficult. Thus, it is desirable to develop computational docking methods that, starting from the coordinates of two unbound component molecules defined as receptor and ligand, respectively, are capable of providing a model of acceptable accuracy for the bound receptor-ligand complex [START_REF] Ritchie | Recent progress and future directions in protein-protein docking[END_REF][START_REF] Andrusier | Principles of flexible proteinprotein docking[END_REF][START_REF] Vajda | Convergence and combination of methods in proteinprotein docking[END_REF][START_REF] Smith | Prediction of protein-protein interactions by docking methods[END_REF]. In view of the large number of putative proteinprotein interactions, the computational efficiency of docking is also a concern.

Most global docking methods start with rigid body search that assumes only moderate conformational change upon the association, accounted for by using a smooth scoring function that allows for some level of steric overlaps [START_REF] Vajda | Convergence and combination of methods in proteinprotein docking[END_REF]. Rigid docking was revolutionized by the fast Fourier transform (FFT) correlation approach, introduced in 1992 by Katchalski-Katzir et al. [START_REF] Katchalski-Katzir | Molecular surface recognition: Determination of geometric fit between proteins and their ligands by correlation techniques[END_REF]. The major requirement of the method is to express the interaction energy in each receptor-ligand orientation as a sum of P correlation functions, i.e., in the form Eðα, β, γ, λ, μ, νÞ = X P p=1 Z R p ðx, y, zÞ Tðλ, μ, νÞ Dðα, β, γÞL p ðx, y, zÞdV , [START_REF] Ritchie | Recent progress and future directions in protein-protein docking[END_REF] where R p and L p are defined on the receptor and ligand, respectively, T and D denote translational and rotational operators, and α, β, γ and λ, μ, ν are the rotational and translational coordinates.

To illustrate how such functions can be used for docking, consider the very simple case with P = 1, R p = -1 on a surface layer and R p = 1 on the core of the receptor, L p = 1 on the entire ligand, and R p = L p = 0 everywhere else. It is clear that this scoring function, which is essentially the one used by Katchalski-Katzir et al. [START_REF] Katchalski-Katzir | Molecular surface recognition: Determination of geometric fit between proteins and their ligands by correlation techniques[END_REF], reaches its minimum on a conformation in which the ligand maximally overlaps with the surface layer of the receptor, thus providing optimal shape complementarity. In later FFT-based methods, the scoring function has been expanded to include electrostatic and solvation terms (6, 7) and, more recently, structure-based interaction potentials [START_REF] Kozakov | PIPER: An FFT-based protein docking program with pairwise potentials[END_REF][START_REF] Mintseris | Integrating statistical pair potentials into protein complex prediction[END_REF], substantially improving the accuracy of docked structures. As mentioned, in all scoring functions, the shape complementarity term allows for some overlaps, thereby accounting for the differences between bound and unbound (separately crystallized) structures.

Most FFT-based methods (6-8, 10-12) define R p and L p on grids, and use a 3D Cartesian FFT approach to accelerate the sampling of the translational space. The method is based on the idea that the energy function, given by Eq. 1, can be expressed in

Significance

Expressing the interaction energy as sum of correlation functions, fast Fourier transform (FFT) based methods speed the calculation, enabling the sampling of billions of putative protein-protein complex conformations. However, such acceleration is currently achieved only on a 3D subspace of the full 6D rotational/translational space, and the remaining dimensions must be sampled using conventional slow calculations. Here we present an algorithm that employs FFT-based sampling on the 5D rotational space, and only the 1D translations are sampled conventionally. The accuracy of the results is the same as those of earlier methods, but the calculation is an order of magnitude faster. Also, it is inexpensive computationally to add more correlation function terms to the scoring function compared with classical approaches.

terms of the Fourier transforms r p of R p and l p of L p . Because the translational operator applied to l p in the Fourier space is given by Tðλ, μ, νÞl p ðn, m, lÞ = e -2πi=Nðnλ+lμ+mνÞ l p ðx, y, zÞ, [START_REF] Andrusier | Principles of flexible proteinprotein docking[END_REF] where i = ffiffiffiffiffiffi -1 p , accounting for the orthonormality of Fourier basis functions and interchanging the order of integration and summation yield Eðα, β, γ, λ, μ, νÞ = X P p=1 X nlm r p ðn, l, mÞl p ðα, β, γ, n, l, mÞe -2πi N ðnλ+lμ+mνÞ , [START_REF] Vajda | Convergence and combination of methods in proteinprotein docking[END_REF] which is the expression for the inverse Fourier transform of the product of the Fourier images r p ðn, m, lÞ and l p ðα, β, γ, n, l, mÞ as stated by the convolution theorem. Thus, for a given rotation, E can be calculated over the entire translational space using P forward and one inverse FFT. If N denotes the size of the grid in each direction, then the efficiency of this approach is OðN 3 logN 3 Þ compared with OðN 6 Þ when energy evaluations are performed directly. Owing to the high numerical efficiency of the FFT-based algorithm, it became computationally feasible, for the first time, to systematically explore the conformational space of protein-protein complexes evaluating the energies for billions of conformations, and thus to dock proteins without any a priori information on the expected structure of their complex. Despite the usefulness of the above algorithm, using FFTs only in translational space has three major limitations. First, FFTs on a new grid must be computed for each rotational increment of the rotating molecule; thus acceleration applies only to half of the degrees of freedom (Fig. 1). Second, each term in the scoring function requires a separate FFT calculation. Thus, accounting for electrostatics, desolvation, and, particularly, pairwise interactions substantially increases the required computational efforts. Third, experimental techniques such as NMR Nuclear Overhauser effect measurements and chemical cross-linking yield information on approximate distances between interacting residues across the interface, and this information can be used to perform the docking subject to pairwise distance restraints. Unfortunately, each pairwise distance restraint requires a new correlation function term. Because the required computational effort is proportional to P, the number of correlation functions in the energy expression, the increasing complexity reduces the numerical advantage of the FFT approach.

In principle, the above problems can be avoided by applying the transforms first, and then moving the proteins in the Fourier space without the need for recomputing the transforms. However, it is difficult to carry out rotations in the translational Fourier space, and, thus, to perform rotations efficiently, it is natural to use spherical coordinates. This approach was applied to crystallography in the early 1970s by Tony Crowther, who realized that the rotation function can be computed more quickly using the FFT, expressing the Patterson maps as spherical harmonics [START_REF] Crowther | The Molecular Replacement Method[END_REF]. A few groups also used this idea for the development of docking algorithms [START_REF] Ritchie | Protein docking using spherical polar Fourier correlations[END_REF][START_REF] Ritchie | Accelerating and focusing protein-protein docking correlations using multi-dimensional rotational FFT generating functions[END_REF]. Most notable is the Hex method of Ritchie and Kemp [START_REF] Ritchie | Protein docking using spherical polar Fourier correlations[END_REF], which represents protein shapes using Fourier series expansions of spherical harmonic and Gauss-Laguerre polynomials. This representation allows rotational searches to be accelerated by angular FFTs, and it enables translations to be calculated analytically in the Fourier basis [START_REF] Ritchie | Accelerating and focusing protein-protein docking correlations using multi-dimensional rotational FFT generating functions[END_REF]. A similar approach has been developed by Chacon's group [START_REF] Kovacs | Fast rotational matching of rigid bodies by fast Fourier transform acceleration of five degrees of freedom[END_REF][START_REF] Garzon | FRODOCK: A new approach for fast rotational protein-protein docking[END_REF], in which translations are calculated numerically. However, both approaches were found to have lower accuracy than traditional Cartesian FFT sampling [START_REF] Ritchie | Accelerating and focusing protein-protein docking correlations using multi-dimensional rotational FFT generating functions[END_REF]. This may be attributed to three main factors. Firstly, the energy functions used were less detailed than in some of the Cartesian approaches. In particular, we used only van der Waals and electrostatic terms [START_REF] Ritchie | Accelerating and focusing protein-protein docking correlations using multi-dimensional rotational FFT generating functions[END_REF]. Secondly, because the computational cost of the polar Fourier translation matrices grows as OðN 5 Þ, the polar Fourier representation is limited to using relatively low order expansions, which limits the achievable accuracy. Finally, the manifold structure of the 5D rotational space was not fully considered, and this resulted in a memory-intensive algorithm that mapped less efficiently onto modern multiprocessor computer architectures than simple 1D FFTs [START_REF] Ritchie | Ultra-fast FFT protein docking on graphics processors[END_REF]. Although we showed previously that the polar representation allows an elegant 5D factorization of multiterm potentials [START_REF] Ritchie | Accelerating and focusing protein-protein docking correlations using multi-dimensional rotational FFT generating functions[END_REF], previous efforts to exploit this property have, until now, had limited success.

In this paper, we describe a fast manifold Fourier transform (FMFT) algorithm that eliminates the above shortcomings and, on Fig. 1. Schematic representation of FFT-based docking methods. In Cartesian FFT sampling (upper path), the ligand protein is translated along three Cartesian coordinates in Fourier space using the translational operator T. The translation must be repeated for each rotation of the ligand. In 5D FMFT docking (lower path), the direction of the vector from the center of the receptor to the center of the ligand is defined by two Euler angles, and the ligand is rotated around its center, resulting in the search space ðSOð3Þ∖SÞ × SOð3Þ. All rotations are performed in generalized Fourier space, where D denotes the rotational operator. The only traditional search is the 1D translation along the vector between the centers of the two proteins.

the average, results in a 10-fold decrease in computing time while retaining the accuracy of the traditional Cartesian FFT-based docking. As will be further emphasized, even more important is that, using FMFT, the computational efforts required are essentially independent of the number of correlation function terms in the scoring function, thus enabling the efficient use of more accurate but also more complex energy expressions, as well as accounting for any number of pairwise distance restraints. Developing the method, we took advantage of the generalization of the Cartesian FFT approach to the rotational group manifold SOð3Þ by Kostelec and Rockmore [START_REF] Kostelec | FFTs on the rotation group[END_REF]. The basis for using this algorithm was recognizing that the 5D rotational search space can be regarded as the product manifold SOð3Þ × ðSOð3Þ∖S 1 Þ, where the rotation group SOð3Þ represents the space of the rotating ligand and ðSOð3Þ∖S 1 Þ is the space spanned by the two Euler angles that define the orientation of the vector from the center of the fixed receptor to the center of the ligand (Fig. 1 and Fig. S1). This is important, because the algorithm by Kostelec and Rockmore [START_REF] Kostelec | FFTs on the rotation group[END_REF] can be easily extended to the SOð3Þ × ðSOð3Þ∖S 1 Þ manifold.

As already mentioned, a general shortcoming of using Fourier decomposition in spherical spaces is the relatively slow convergence of the series of spherical basis functions. Thus, using a large number of terms reduces computational efficiency, whereas truncating the series limits the accuracy of the energy values calculated by the method. Therefore, a key factor explaining the success of our manifold FFT docking method is that we select the centers of highly populated clusters of low-energy docked structures rather than simply low-energy conformations as predictions of the complex. Such clusters occur in low-energy regions around the local minima in the conformational space. The size of each cluster represents the width of the corresponding energy well, and hence provides some information on entropic contributions to the free energy. Model selection based on cluster size has been used in our very successful docking server ClusPro and, in a substantial fraction of docking problems, enabled the identification of the docked structure closest to the native complex [START_REF] Kozakov | How good is automated protein docking?[END_REF]. We note that a similar clustering step is implemented in the protein structure prediction program Rosetta [START_REF] Gray | Protein-protein docking with simultaneous optimization of rigidbody displacement and side-chain conformations[END_REF]. For a somewhat more formal justification of the cluster-based approach to model selection, we argue that, using FFT, we globally and systematically sample the energy landscape of the interacting protein pair on a grid, and hence we can calculate an approximate partition function of the form Z = P j expð-E j =RTÞ, where E j is the energy of the jth pose, and we sum over all poses. For the kth low-energy cluster, the partition function is given by Z k = P j expð-E j =RTÞ, where the sum is restricted to poses within the cluster. Based on these values, the probability of the kth cluster is given by P k = Z k =Z. However, because the low-energy structures are selected from a relatively narrow energy range, and the energy values are calculated with considerable error, it is reasonable to assume that these energies do not differ from each other, i.e., E j = E for all j in the low-energy clusters. This simplification implies that P k = expð-E=RTÞ × ðN k =ZÞ, and thus the probability P k is proportional to N k , where N k is the number of structures in the kth cluster. Therefore, we select the centers of highly populated clusters of docked structures, rather than low-energy conformations, as predictions of the complex. Although neglecting the energy differences within the low-energy clusters seems to be arbitrary, the success of the ClusPro server demonstrates that the approximation is valid in a large fraction of cases. The significance of model selection based on cluster size rather than energy values is that it does not require very accurate energy evaluation, and hence, in FMFT, it is sufficient to use a limited number of spherical basis functions in the Fourier space, increasing numerical efficiency without noticeable loss of docking accuracy.

The high efficiency of the FMFT algorithm enables solving very demanding docking problems, way beyond what was considered feasible in the past. After demonstrating that the accuracy of FMFT is comparable to that of the traditional Cartesian FFTbased docking, we present here a few applications that require a large number of docking calculations. Such problems include docking ensembles of models obtained by NMR or homology modeling, and exploring a large number of putative peptide conformations in peptide-protein docking. As will be described, an additional and very favorable property of the FMFT algorithm is that the required computational efforts are almost completely independent of the number P of the correlation function terms in the energy expression given by Eq. 1, and hence the method can be efficiently used with scoring functions of arbitrary complexity. In contrast, in the traditional FFT approach, the efforts are proportional to P, and hence it is difficult to perform docking subject to pairwise distance restraints, as each restraint gives rise to an additional term in the scoring function. Using FMFT, we demonstrate that this problem can be solved effectively without significant increase in running times (Fig. S2).

Results and Discussion

FFT-Based Docking on 5D Rotational Manifolds. Here we demonstrate that, by taking advantage of the special geometry of the space characterizing molecular movement upon protein-protein association, it is possible to construct an extremely efficient FFT-based docking algorithm. We present the basic idea of this algorithm as the generalization of the translational FFT method described in the Introduction. Because we plan to work in the rotational space, we change the Cartesian coordinates to polar coordinates ðx, y, zÞ → ðr, θ, ϕÞ, and consider the generalization of the Fourier transform on the sphere Rðr, θ, ϕÞ = X N nlm rðn, l, mÞR nl ðrÞd lm ðcosθÞe -imϕ [START_REF] Smith | Prediction of protein-protein interactions by docking methods[END_REF] where R nl ðrÞ are radial basis functions, rðn, l, mÞ are generalized Fourier coefficients, d lm ðcosθÞ are Legendre polynomials [START_REF] Zare | Angular Momentum: Understanding Spatial Aspects in Chemistry and Physics[END_REF], and N is the order of expansion used. Eq. 4 looks like a Fourier transform, but e -imϕ is replaced by d lm ðcosθÞe -imϕ , which shows the non-Cartesian properties of the sphere [START_REF] Driscoll | Computing Fourier transforms and convolutions on the 2-sphere[END_REF].

Consider again the derivation of the convolution theorem (Eq. 1) but, this time, on the manifold ðSOð3Þ∖SÞ × SOð3Þ shown in the lower path of Fig. 1. The translation of the ligand can be represented as the rotation of the receptor, followed by the translation of the ligand along the z axis, Eðz, β, γ, α′, β′, γ′Þ =

X P p=1 Z
Tð-zÞ Dð0, β, γÞR p ðρ, θ, ϕÞ × Dðα′, β′, γ′ÞL p ðρ, θ, ϕÞdV .

[5]

Rotations of the receptor can be expressed as follows:

Dðα, β, γÞRðr, θ, ϕÞ = X nlm R nl ðrÞY lm ðθ, ϕÞ X m1 D l mm 1 ðα, β, γÞrðn, l, m 1 Þ, [6]
where Y lm ðθ, ϕÞ denotes spherical harmonics, and

D l mm′ ðα, β, γÞ = e -imα d l mm′ ðβÞe -im′γ [7]
are Wigner rotation matrices with d l mm1 ðβÞ denoting Wigner d functions, related to Jacobi polynomials [START_REF] Kostelec | FFTs on the rotation group[END_REF]. Eqs. 6 and 7 show that the rotational operator in the rotational group SOð3Þ acts on generalized Fourier coefficients the same way as the translation operator acts on Fourier coefficients in the Cartesian space (Eq. 2), apart from the asymmetry of the middle angle β, which requires special treatment. Describing the translation of the ligand along the z axis in the Fourier space is far from simple, and requires updating a set of coefficients. However, it is only 1 degree of freedom (as opposed to 3 degrees in the Cartesian space), and hence it can be accomplished relatively efficiently [START_REF] Ritchie | High-order analytic translation matrix elements for real-space sixdimensional polar Fourier correlations[END_REF]. Now we apply the translation operator and the rotation operator (Eq. 7) to the integral in Eq. 5. Based on the orthonormality of the generalized Fourier basis functions, interchanging the order of integration and summation yields Eðz, β, -γ, α′, β′, γ′Þ

= X mm 1 m2 ll 1 X nn 1 X p r p ðn 1 , l 1 , m 1 Þl p ðn, l, m 2 ÞT jmj nln 1 l 1 ðzÞ ! × d l 1 mm 1 ðβÞd l mm2 ðβ′Þe -iðmα′+m 1 γ+m2γ′Þ
.

[8]

Note that Eq. 8 is similar to Eq. 3 in Cartesian coordinates, with the difference that, instead of a 3D inverse Fourier transform, we have a generalized FMFT, which involves the Wigner d functions d l mm1 ðβÞ. However, the really important difference is in the order of the transforms and the summation of correlation functions. In Eq. 3, for each rotation of the ligand, we have to calculate the Fourier transforms l p ðα, β, γ, n, l, mÞ for each of the P components of the ligand energy function separately, form the product with the transform r p ðn, m, lÞ of the pth component of the receptor energy function, sum all terms, and take the inverse transform. In contrast, according to Eq. 8, we calculate the sum of initial precalculated generalized Fourier coefficients in the internal loop only once, and perform all rotations in Fourier space rather than calculating an FFT for each rotation. This allows us to calculate multiple energy terms using a single FMFT for each translation. Thus, as already emphasized, the computational efforts are essentially independent of the number P of the correlation function terms in the energy expression. Because inverse manifold Fourier transforms can be efficiently calculated by methods due to ref. [START_REF] Kostelec | FFTs on the rotation group[END_REF], this approach provides substantial computational advantage, particularly if P is high.

Execution Times. Execution times of the FMFT sampling algorithm were measured by docking unbound structures of component proteins in 51 enzyme-inhibitor pairs from the established Protein Docking Benchmark [START_REF] Hwang | Protein-protein docking benchmark version 4.0[END_REF] (Table S1). The times were compared with those required for docking the same proteins using PIPER, a protein docking program based on the Cartesian FFT approach [START_REF] Kozakov | PIPER: An FFT-based protein docking program with pairwise potentials[END_REF]. The FFTW (Fastest Fourier Transform in the West) library (26) was used for FFT calculations. All runs were performed using the standard PIPER scoring function, consisting of eight correlation function terms. Execution times were measured on one or several Intel Xeon E5-2680 processors. Using the FMFT algorithm, the average execution time was 15.39 min. In comparison, the average execution time for the same set of proteins using PIPER was 232.15 min, indicating that FMFT speeds up the calculations ∼15-fold. Using parallel versions of the algorithms on 16 CPU cores, the average execution times measured were 2.67 min and 20.19 min for FMFT and PIPER, respectively, which shows about a 7.5-fold speedup.

Application 1: Constructing Enzyme-Inhibitor Complexes. The quality of FMFT and PIPER results was determined by docking the same 51 enzyme-inhibitor pairs that we have used for comparing execution times (Table S2). In both cases, the scoring function was the same one normally used in PIPER for docking enzyme-inhibitor pairs, and it consisted of attractive and repulsive van der Waals, Coulombic electrostatics, generalized Born, and knowledge-based Decoys As the Reference State terms, the latter representing nonpolar solvation [START_REF] Chuang | DARS (Decoys As the Reference State) potentials for protein-protein docking[END_REF]. The docking procedure for these cases was the one normally used by PIPER [START_REF] Kozakov | How good is automated protein docking?[END_REF]. First, the conformational space was sampled using either the FMFT or the PIPER protocol. After docking, the 1,000 lowest energy poses were retained and clustered using interface C α rms deviation (RMSD) as the distance metrics with a fixed 9 Å clustering radius. The clusters were ranked according to cluster populations (i.e., number of poses in the cluster), and the centers of up to 30 largest clusters were reported as putative models of the complex (Table S2).

Fig. 2 A-C shows the results of docking. The number of hits shown in Fig. 2A is the number of near-native poses, defined as having less than 10 Å C α interface RMSD (IRMSD) from the native complex, generated by each of the two algorithms. Note that IRMSD is calculated for the backbone atoms of the ligand that are within 10 Å of any receptor atom after superimposing the receptors in the X-ray and docked complex structures. We found that the number of poses with less than 10 Å IRMSD is a good measure of the quality of sampling of the energy landscape in the vicinity of the native structure. Fig. 2 B and C shows the properties of models obtained by clustering low-energy poses using pairwise IRMSD as a distance metric. A large number of lowenergy poses typically yields a well-populated and thus highly ranked near-native cluster, reported as one of the final models. Based on all these results, FMFT and PIPER show comparable docking performance, both in terms of the number of near-native structures (Fig. 2A), the ranks of the clusters that define the final near-native models (Fig. 2B), and the IRMSD (Fig. 2C) of these models.

Application 2: Docking Interacting Protein Domains. We further compared FMFT and PIPER by docking interacting domains extracted from proteins that are defined as "Other" type in the Protein Docking Benchmark (25) (Tables S3 andS4). This problem is generally more challenging than docking inhibitors to enzymes because the Other category includes complexes with highly variable properties. Restricting consideration to individual domains eliminates the additional problem that the domains in multidomain proteins may shift relative to each other, affecting the docking results. Thirty cases representing domain-domain binding were selected from the Others section of the Protein Docking Benchmark (Table S4). Nineteen cases from this set represent binding of singledomain proteins (or single domains taken from larger proteins), and thus full protein structures were used for docking. In another 11 cases, receptor and/or ligand are composed of several domains, so reduced representations of protein structures were prepared: Only the binding domains were retained for docking, and the rest of the structure was cleaved. Residue ranges for binding domains were assigned according to structural classification of proteins (SCOP) domain classification [START_REF] Murzin | SCOP: A structural classification of proteins database for the investigation of sequences and structures[END_REF]. To prevent possible association at intraprotein domain-domain binding interfaces exposed by the cleavage, additional repulsion grids were used in the docking procedure. These were constructed by taking the backbone atoms of the original structure lying within 10 Å (but not closer than 5 Å) of the binding domain and placing repulsive spheres with 0.5-Å radius at the positions of those atoms. The 5-Å lower bound to the distance range specifying the thickness of this "repulsive padding" was introduced to ensure that additional repulsion doesn't affect binding to the relevant portion of protein surface. During the docking process, such repulsive padding grid was correlated with the standard repulsive van der Waals grid of the binding partner. The docking procedure overall was the same as that used for enzyme-inhibitor targets, except that 1,500 low-energy poses were used for clustering, generated from three docking runs (500 poses from each) performed with differently weighted components of the scoring function [START_REF] Kozakov | How good is automated protein docking?[END_REF]. Similarly to the results obtained for enzyme-inhibitor complexes, FMFT and PIPER show comparable performance (Fig. 2 D-F and Table S3). Although PIPER generates large numbers (>200) of nearnative structures for more complexes than FMFT, the number of complexes with very few (<10) such near-native structures is substantially smaller using FMFT than using PIPER. Thus, FMFT shows better performance for the more difficult-to-dock complexes (Fig. 2D). In addition, using PIPER, the number of models that are not ranked in the top 10 is much higher than using FMFT (Fig. 2E). Based on these results, FMFT performs as well as PIPER.

Application 3: Accounting for Pairwise Distance Restraints. An important consideration for selecting a docking method is the maximum complexity of the scoring function that still allows for solving problems with reasonable execution times. As mentioned, all FFTbased approaches require the use of scoring functions that can be written as sums of correlation functions. This is not a major limitation, because such functions may include many commonly used physics-based energy terms, such as steric repulsion, van der Waals interaction, and Coulombic electrostatics. It has also been shown that some energy terms that are not inherently correlation-based, such as the widely used pairwise interaction potentials, can be efficiently approximated by a sum of several correlation functions [START_REF] Chuang | DARS (Decoys As the Reference State) potentials for protein-protein docking[END_REF]. Altogether, this makes the number of correlations a crucial parameter, because this number effectively defines the complexity of the scoring function in the particular sampling run.

One important task, especially demanding in terms of scoring function complexity, is incorporating pairwise distance restraints, based on known interactions between residue pairs, into the docking procedure. Such restraints can be derived in a variety of experiments, including NMR, cross-linking, and mutagenesis assays [START_REF] Dominguez | HADDOCK: A protein-protein docking approach based on biochemical or biophysical information[END_REF]. The restraints can be implemented as short-distance attractive terms in the scoring function, but each will add a correlation function term. As emphasized, in Cartesian FFT, the number of transforms required is proportional to the number P of correlation functions (Eq. 3), whereas, in FMFT, the number of transforms is independent of P. To demonstrate this difference, we determined the structure of the glucose-specific enzyme IIA (E2A)-histidinecontaining phosphocarrier protein (HPr) complex [Protein Data Bank (PDB) entry 1GGR] (30) from the structures of its constituents in their unbound form (PDB entries 1F3G and 1POH) and 20 ambiguous interaction restraints (AIRs) based on NMR titration data [START_REF] Dominguez | HADDOCK: A protein-protein docking approach based on biochemical or biophysical information[END_REF]. The docking procedure was the one used for docking enzyme-inhibitor pairs, but with 20 additional correlations terms in the scoring function due to the restraints [START_REF] Dominguez | HADDOCK: A protein-protein docking approach based on biochemical or biophysical information[END_REF]. Each restraint is specified as a residue in one of the proteins, and a set of residues on the partner protein that are in contact with the first residue, where "contact" means ≤3 Å distance between any two atoms of the residue pair. To represent these restraints, receptor and ligand correlation components were constructed by placing 3-Å radius attractive spheres on the atoms of the particular residue on the first protein and the attracted point "charges" on the atoms of the interacting residues on the partner protein. Docking was performed using both FMFT and PIPER. Incorporation of restraints increased the population of the near-native cluster from 201 to 410, which became the most populated cluster and thus provided the putative model of the complex (Fig. 3 and Table S5) without any significant change in the IRMSD of the cluster center (5.25 Å for the unrestrained case versus 5.15 Å for the restrained). Adding the restraints increased the number of correlation function terms in the scoring function from 8 to 28. For PIPER, this resulted in a proportional increase in execution time (from 96.15 min to 373.80 min). In contrast, running FMFT, the execution time barely changed, from 12.32 min to 15.30 min. This result demonstrates that FMFT can be used with very complex scoring functions (Fig. S2).

Application 4: Docking Ensembles of NMR Models. Multiple docking runs may be required when one or both component proteins are given as ensembles of structures, obtained by NMR experiments or by extracting snapshots from molecular dynamics simulations. Because accounting for multiple structures may substantially improve docking results, the high efficiency of the FMFT method is particularly useful. As an example, we considered calculating the complex formed by the Escherichia coli Colicin E9 DNase domain and its cognate immunity protein IM9. Four different X-ray structures of the unbound E9 DNase domain (chains B, C, D, and E of PDB entry 1FSJ) were docked in a pairwise manner to 20 NMR models of the IM9 protein (PDB entry 1E0H), thus performing 80 docking calculations. Unstructured termini of the receptor were masked and didn't contribute to the calculated energy scores. The 50 lowest energy poses were extracted from each of the 80 docking runs and merged, yielding a total of 4,000 poses that were then clustered as usual in PIPER. Fig. 4 A-C shows the docking results. In short, merging the 50 lowest energy poses from each docking run, followed by clustering, provided a 2.94-Å IRMSD model of the complex ranked fifth, where the 1IBX structure of the native complex was used for comparison when evaluating the accuracy of results. To emphasize the advantage of ensemble docking, we also docked a single pair of structures, chain B of 1FSJ and the first NMR model of the ligand from the PDB file 1E0H. The standard docking protocol was used, and we retained the 1,000 lowest energy poses for clustering. Docking the single pair, the best near-native model obtained was ranked 13, and had the IRMSD value of 3.45 Å. Thus, in the case of structural uncertainty of the component proteins, ensemble docking can substantially improve the results, and, in this type of application, the higher speed of FMFT is a major advantage. Computational efficiency will be particularly important for genome-wide analyses of protein-protein interactions, but we think that, for such applications, it will be necessary to better understand the docking of homology models (see Application 5: Identification of Binding Sites by Docking Homology Models), because, generally, structures are not available for a substantial fraction of proteins.

Application 5: Identification of Binding Sites by Docking Homology Models. It has been shown that protein-protein interaction sites can be found by determining the highly populated interfaces in the ensemble of structures generated by global docking [START_REF] Hwang | Binding interface prediction by combining protein-protein docking results[END_REF][START_REF] Fernández-Recio | Identification of protein-protein interaction sites from docking energy landscapes[END_REF]. We implemented this approach by clustering the "interfacial" atoms in the low-energy docked poses. Although this method usually requires structures of the component proteins, we extended the approach to proteins with yet undetermined structures by docking multiple homology models. The extended method was applied to determining the interface in the Nef-Fyn(R96I)SH3 complex (PDB entry 1EFN). Ten models of the receptor (SH3 domain) and 2 models of the ligand (HIV-1 Nef protein) were constructed using the MODELLER program [START_REF] Eswar | Comparative protein structure modeling using MODELLER[END_REF] and based on homologous templates with 30-60% sequence identity (see Table S6 for the list of templates used). All possible receptor-ligand model pairs were docked using the approach developed for Other type of complexes [START_REF] Kozakov | How good is automated protein docking?[END_REF]. From each of the 20 docking runs, we selected the 1,500/20 = 75 lowest energy poses that were merged and clustered using RMSD as the distance metrics. The structures at the centers of these clusters were used to define interface atoms as atoms located within 5 Å of any atom of the partner protein. These interfacial atoms were then subjected to bottom-up hierarchial clustering using the Euclidian distance as the metrics. Clustering was terminated, i.e., neighboring clusters were not merged, if the minimal distance between a pair of their atoms was larger than the value of a separation parameter. The resulting clusters were ranked according to cluster population (i.e., the number of atoms in each cluster), and the largest cluster was considered to be the most probable prediction of the protein-protein interaction site. For comparison, we also predicted the interaction site by docking a single pair of homology models based on the templates with the highest sequence identity. In this case, a slightly larger value of the clustering separation parameter was used (1.35 Å rather than 1.30 Å). This change was due to the fact that a single docking run provided fewer interfacial atoms for hierarchical clustering, resulting in clusters that were too small. Therefore, the value of the cutoff parameter was increased to ensure that the relative population of the largest cluster was comparable to that obtained by merging the results from 20 docking runs. As shown in Fig. 4 D and E, docking of multiple homology models of the component proteins increased the accuracy of binding site prediction, compared with the result of using the maximum sequence identity models alone.

Application 6: Docking Flexible Peptides. The difficulty in docking short linear peptides is that their structure in solution is generally unknown and may be ill-defined. One possible solution is to dock a variety of peptide conformations, thus requiring multiple docking runs. We have recently developed an algorithm based on the use of structural templates extracted from the PDB with sequences that matched the known sequence motif in the peptide. These templates were docked individually using the FMFT algorithm. From each run, a number of low-energy poses were retained, the pooled peptide structures were clustered, and the highly populated cluster centers were reported as final models as in all applications of our docking algorithm.

Here we demonstrate this algorithm by docking the ace-PQQATDD peptide to the tumor necrosis factor receptor-associated factor 2 (TRAF2). For this peptide, the PXQ motif sequence known from the literature was extended to length 7 (PXQXXDD) and used to extract 316 structural templates from the PDB database. These templates were then used to model the target peptide. The models were aligned and clustered using the backbone RMSD as the distance measure, with 0.5 Å as the fixed clustering radius. Peptide structures corresponding to the centers of the 25 mostpopulated clusters were docked to the unbound receptor structure (chain A of PDB entry 1CA4).

The 250 lowest energy poses were retained from each docking run. The poses were merged and clustered using backbone RMSD as a distance measure with 3.5 Å as the fixed clustering radius. Cluster centers were ranked according to cluster populations and reported as final models. Docking results were evaluated using the backbone RMSD from the structure of the peptide in the native complex (chain A of PDB entry 1CZY). A near-native model of the protein-peptide complex was ranked fourth and had the backbone RMSD of 3.3 Å from the conformation in the X-ray structure (Table S7 and Fig. 5A). Note that docking only the most frequently occurring structural template provides less accurate models, as demonstrated in Fig. 5 B andC.

Conclusions

Extending the classical 3D Cartesian implementation of the FFT correlation approach to perform rotations in Fourier space without the need for recalculating the transforms has been a long-outstanding and extensively studied problem. The main difficulty in developing such methods is that, to achieve numerical efficiency, one can use only a moderate number of spherical basis functions to span the search space, and this may reduce the accuracy of energy evaluation. However, because we base model selection on the population of low-energy clusters rather than on energy values, minor deviations in energy generally do not affect the accuracy of final models. Here we present an elegant manifold FFT implementation of 5D search that is more than 10-fold faster than the traditional 3D approach. A major advantage of the method is that adding correlation function terms in the scoring function is computationally inexpensive, and hence the method works efficiently with very complex energy evaluation models, possibly including pairwise distance restraints that are difficult to deal with in traditional FFT-based docking. The improved efficiency implies that we can solve new classes of docking problems, including the docking of large ensembles of proteins rather than just a single protein pair, docking homology models, and flexible peptides that may have a large number of potential conformations. We note that the beta version of a code implementing the FMFT algorithm can be downloaded from https://bitbucket.org/abcgroup_midas/fmft_dock/, thus providing an opportunity for testing and using the method. In addition, we are in the process of adding FMFT as a new option to the server.

Materials and Methods

This section summarizes the implementation of the FMFT approach. For the mathematical details of the algorithm, see SI Materials and Methods.

The procedure starts with receptor-and ligand-associated components of each correlation term of the energy function being represented as sets of coefficients rðn, l, mÞ, lðn, l, mÞ that appear in the expansion shown as Eq. 4.

Here 1 ≤ n ≤ N, 0 ≤ l ≤ n -1, and -l ≤ m ≤ + l, where N governs the order at which the series is truncated. These coefficients, together with the translation range to be sampled (i.e., minimal and maximal distances between protein centers, calculated from the geometrical properties of the proteins), are submitted as input parameters to the program performing the FMFTbased sampling. To improve efficiency, two stages of FMFT sampling are being executed: The first one, performed with a maximal coefficient order N = 20 on a small FFT grid, is computationally inexpensive and provides a crude approximation of the energy landscape, which is then used to focus the search to the translation range potentially containing the energy minima, whereas the second one is executed with N = 30 on a full-sized FFT grid but performs the sampling only in the refined translation range, thus saving computational resources.

The actual sampling stage can be described as follows: After loading the input parameters, the program starts to iterate the allowed translation range in steps of 1 Å. For each translation step, the P nn1 P p r p ðn 1 , l 1 , m 1 Þl p ðn, l, m 2 ÞT jmj nln1 l1 ðzÞ product of coefficients and translation matrix elements is calculated, followed by a manifold FFT, which provides the values of energy score for all receptor-ligand orientations corresponding to a fixed distance between the centers of the two proteins. The resulting samples are located on the ðβ, γ, α′, β′, γ′Þ Euler angle grid with dimensions of 30 × 59 × 59 × 30 × 59 (or 16 × 30 × 30 × 16 × 30 for the loworder scan). K (on the order of 1,000 for a typical sampling run) lowest energy samples are retained for each translation step. After the entire translation range is processed, the low-energy samples from individual translation steps are merged and resorted by energy value to select the K lowest energy samples that are presented as the final results.

It is important to note here that the sampling of the S 2 × SOð3Þ manifold [in practice probed as ðSOð3Þ∖S 1 Þ × SOð3Þ], provided by the equispaced sampling of Euler angles, is inherently nonuniform. This becomes a significant problem if one seeks to obtain statistical information about the energy landscape of protein interaction, for example, to construct the partition function of the system. To battle this nonuniformity, a special procedure is used for the selection of low-energy scores. Specifically, once the 5D array of energy scores for a single translation step is acquired, the program starts selecting lowest-scoring conformations and ex-cluding the samples corresponding to the surrounding region from further consideration. Here the "surrounding region" is defined as the subset of elements fðx, yÞjxðβ, γÞ ⊂ S 2 , yðβ, γ, α′, β′, γ′Þ ⊂ SOð3Þg of the S 2 × SOð3Þ manifold, for which ðdist S 2 ðx, x min Þ < ΔÞ∧ðdist SOð3Þ ðy, y min Þ < ΔÞ, where Δ is a cutoff parameter chosen to be 6.0°, which is slightly less than the grid step of 360°=59 = 6.1°. This procedure ensures that the sampling explores a substantial fraction of the conformational space rather than producing structures very close to each other.

Fig. 2 .

 2 Fig. 2. Results of docking enzyme-inhibitor and domain-domain pairs. Bar heights represent the number of docking cases that fall into an appropriate category. (A) The number of hits among the 1,000 low-energy poses generated for enzyme-inhibitor complexes. (B) Ranking of final near-native models for enzyme-inhibitor complexes. (C) C α IRMSD of the final model for enzyme-inhibitor complexes (here only cases with both FMFT and PIPER producing a nearnative model were taken into account). (D) The number of hits among the 1,500 low-energy poses generated for domain-domain complexes. (E) Ranking of final near-native models for domain-domain complexes. (F) C α IRMSD of the final model for domain-domain complexes. As in C, only cases with both FMFT and PIPER producing a near-native model were taken into account.

Fig. 3 .

 3 Fig. 3. Docking of E2A and HPr proteins. (A) Model defined by the most populated cluster obtained without restraints. (B) Model defined by the most populated cluster obtained with restraints. A set of cyan cylinders represents one of the 20 restraints. (C) IRMSD versus energy score for docking without restraints. (D) IRMSD versus energy score for docking with restraints. Incorporation of experimental restraints substantially increased the population of the near-native cluster.

Fig. 4 .

 4 Fig. 4. Docking of structural ensembles. (A) Sampling the interaction energy landscape using a single E9 DNase domain structure and the first NMR model of IM9. The docking does not capture any near-native energy minimum. (B) Consensus energy values from the 80 pairwise dockings of four different X-ray structures of the E9 DNase domain to 20 NMR models of the IM9 protein. (C) Cartoon representation of the four E9 DNase domain and 20 IM9 structures used for docking, superimposed on the structure of the native complex (gray shade). (D) Binding site identification for the Nef-Fyn(R96I)SH3 complex obtained by docking the highest sequence identity models alone. (E) Using multiple homology models of the receptor and the ligand to identify the binding site for the Nef-Fyn(R96I) SH3 complex results in a more specific prediction.

Fig. 5 .

 5 Fig. 5. Docking of the ace-PQQATDD peptide to TRAF2. (A) Bound structure of the peptide (red) and the 3.3-Å model, ranked fourth (cyan). (B) Peptide backbone RMSD versus scoring function when docking the most common structural template alone. (C) Peptide backbone RMSD versus scoring function when using all 25 templates. Docking the ensemble substantially improves the results, and yields samples with less than 4.0-Å backbone RMSD.

E4288 | www.pnas.org/cgi/doi/10.1073/pnas.1603929113 Padhorny et al.

E4290 | www.pnas.org/cgi/doi/10.1073/pnas.1603929113 Padhorny et al.

ACKNOWLEDGMENTS. This work was supported by National Science Foundation (NSF) Grants AF 1527292 and DBI 1458509, NIH, National Institute of General Medical Sciences, Grants R35 GM118078 and R01 GM093147, Russian Scientific Foundation Grant 14-11-00877, and US Israel Binational Science Foundation Grant 2009418.


The second example is referred to the fact that, unlike classical Cauchy continua, second and higher order continua can respond to concentrated forces and generalized contact actions (see e

Formulation of the problem

The general case

The coordinates X , in the reference configuration, are those of the 2-dimensional body B. The internal energy density functional U ( G , ∇G ) depends not only on the strain G = F T F -I / 2 but also on its gradient ∇G , where F = ∇χ, χ is the placement function. In [START_REF] Mindlin | Micro-structure in linear elasticity[END_REF] a general form of the density of the strain energy functional of a linear isotropic second gradient elastic material is given, for the sake of simplicity, in indicial notation,

U ( G, ∇G ) = λ 2 G ii G j j + μG i j G i j + 4 α 1 G aa,b G bc,c + α 2 G aa,b G cc,b + 4 α 3 G ab,a G cb,c + 2 α 4 G ab,c G ab,c + 4 α 5 G ab,c G ac,b , (1)
where subscript j after comma indicates derivative with respect to X j and a general rule for index notation is the following: the subscript-indices of a symbol denoting a vector or a tensor quantity denote the components of that quantity. In the 2-dimensional case we have

U ( G, ∇G ) = ˜ U ( u ) = ( λ + 2 μ) u 2 1 , 1 + u 2 2 , 2 + μ u 2 1 , 2 + u 2 2 , 1 + 2 λu 1 , 1 u 2 , 2 + 2 μu 1 , 2 u 2 , 1 + 1 2 A u 2 1 , 22 + u 2 2 , 11 + 1 2 B u 2 1 , 11 + u 2 2 , 22 + C u 2 1 , 12 + u 2 2 , 12 + 2 D ( u 1 , 11 u 2 , 12 + u 2 , 22 u 1 , 12 ) + 1 2 ( A + B -2 C ) ( u 1 , 11 u 1 , 22 + u 2 , 11 u 2 , 22 ) + ( B -A -2 D ) ( u 1 , 12 u 2 , 11 + u 1 , 22 u 2 , 12 ) , (2)
where A = 2 α 3 + 2 α 4 + 2 α 5 , B = 8 α 1 + 2 α 2 + 8 α 3 + 4 α 4 + 8 α 5 .

(3)

C = 2 α 1 + α 2 + α 3 + 3 α 4 + 5 α 5 , D = 3 α 1 + α 2 + 2 α 3 , ( 4
)
where u is the displacement field, λ and μ are the Lamè coefficients and α i with i = 1 , 2 , 3 , 4 , 5 are the 5 second gradient constitutive parameters. Note that here we use the Lamè coefficients λ and μ to describe first-gradient, isotropic linear elasticity, but other choices could be made, e.g., the pair comprised of bulk modulus κ and shear modulus μ [START_REF] Federico | A transversely isotropic composite with a statistical distribution of spheroidal inclusions: A geometrical approach to overall properties[END_REF][START_REF] Federico | The linear elasticity tensor of incompressible materials[END_REF][START_REF] Federico | Considerations on incompressibility in linear elasticity[END_REF][START_REF] Hill | A self-consistent mechanics of composite materials[END_REF][START_REF] Walpole | Elastic behavior of composite materials: Theoretical foundations[END_REF][START_REF] Walpole | Fourth-rank tensors of the thirty-two crystal classes: Multiplication tables[END_REF], which is particularly convenient, e.g., when treating quasi-incompressibility.

In [START_REF] Mindlin | Micro-structure in linear elasticity[END_REF] , in order to have the positive definiteness of U , the following restrictions on the 7 constitutive parameters must be satisfied,

μ > 0 , 3 λ + 2 μ > 0 , -4 α 1 + 2 α 2 + 2 α 3 + 6 α 4 -6 α 5 > 0 , α 4 > α 5 , α 4 + 2 α 5 > 0 (5) 4 α 1 + α 2 + 4 α 3 + 2 α 4 + 4 α 5 > 0 , α 1 + α 2 < α 3 , 4 α 1 -2 α 2 -2 α 3 -3 α 4 + 3 α 5 > 0 .
In the 2-dimensional case the positive definiteness of U is implied again by the classical 2-dimensional restrictions

μ > 0 , λ + μ > 0 ,
and by the positive definiteness of the following matrix

⎛ ⎜ ⎜ ⎜ ⎜ ⎝ A 0 1 2 ( A + B -2 C ) 0 0 B -A -2 D 0 A 0 1 2 ( A + B -2 C ) B -A -2 D 0 1 2 ( A + B -2 C ) 0 B 0 0 2 D 0 1 2 ( A + B -2 C ) 0 B 2 D 0 0 B -A -2 D 0 2 D 2 C 0 B -A -2 D 0 2 D 0 0 2 C ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ .
Keeping this in mind, a classical variational procedure gives the following system of partial differential equations

∀ X i ∈ B, u 1 , 11 ( λ + 2 μ) + u 1 , 22 μ + u 2 , 12 ( λ + μ) = u 1 , 1111 B + u 1 , 2222 A + u 1 , 1122 ( A + B ) + ( u 2 , 1222 + u 2 , 1112 ) ( B -A ) -b ext 1 (6) u 2 , 22 ( λ + 2 μ) + u 2 , 11 μ + u 1 , 12 ( λ + μ) = u 2 , 2222 B + u 2 , 1111 A + u 2 , 1122 ( A + B ) + ( u 1 , 1222 + u 1 , 1112 ) ( B -A ) -b ext 2 , ( 7
)
and boundary conditions given ∀ X i ∈ ∂B from the following duality conditions

δu α t α -t ext α = 0 , δu α, j n j τ α -τ ext α = 0 , [ ∂ ∂ B ] δu α f α = [ ∂ ∂ B ] δu α f ext α , (8)
∂ B = m c=1 c , [ ∂ ∂ B ] = m c=1 V c ,
then the line and vertex-integrals of a generic field g ( X i ) are represented as follows,

∂B g ( X i ) = m c=1 c g ( X i ) , [ ∂ ∂ B ] g ( X i ) = m c=1 g X c i (9)
where X c i is the coordinate of the vertex V c . Moreover, the so called contact force t α , contact double force τ α and contact wedge force f α are defined, t α = S α j -T α jh,h n j -P ka T αh j P ah n j ,k (10)

τ α = T α jk n j n k (11) f α = T αhk V hk ( 12
)
where n i is the normal to the boundary ∂B, P ij is its tangential projector operator ( P i j = δ i jn i n j ), V is the vertex operator

V h j = ν l h n l j + ν r h n r j ,
where superscripts l and r refers (roughly speaking, left and right), respectively, to one and to the other sides that define a certain vertex-point V c ; ν is the external tangent unit vector. Stress and hyper stresses are defined,

S i j = ∂U ∂G i j , T i jh = ∂U ∂G i j,h . ( 13
)
2.2. Rectangles

The general case of straight lines

In the case of boundaries ∂B composed of straight-lines, the contact force in (10) , the contact double force in (11) and the contact wedge force (12) are

t α = S α j n j -T α jh,h + T αh j,h n j + T αh j,k n h n k n j , τ α = T α jk n j n k , f α = T αij V i j , (14)
that, in terms of the displacement fields, yield, t α = λu a,a n α + μu α, j n j + μu j,α n ju a,abb n α (

6 α 1 + 2 α 2 + 4 α 3 ) (15) -u a,aαk n k ( 6 α 1 + 2 α 2 + 4 α 3 + 2 α 4 + 8 α 5 ) -u α,aak n k ( 2 α 3 + 4 α 4 + 6 α 5 ) -u k,αaa n k ( 2 α 1 + 2 α 3 + 2 α 4 + 6 α 5 ) + u a,a jk n α n j n k ( 4 α 1 + 2 α 2 + 2 α 3 ) + u j,aak n α n j n k ( 2 α 1 + 2 α 3 ) + u α,abc n a n b n c ( 2 α 4 + 2 α 5 ) + u a,αbc n a n b n c ( 2 α 4 + 6 α 5 ) , τ α = u a,ab n α n b ( 4 α 1 + 2 α 2 + 2 α 3 ) + u a,bb n α n a ( 2 α 1 + 2 α 3 ) (16) + ( 2 α 1 + 2 α 3 ) u a,aα + u α,ab n a n b ( 2 α 4 + 2 α 5 ) + 2 α 3 u α,aa + u a,αb n a n b ( 2 α 4 + 6 α 5 ) .
We remark that the formulation expressed in ( 15) and ( 16) can also be used in the 3-dimensional case. This is the reason why (15) and ( 16) are expressed in terms of the 5 3-dimensional constitutive coefficients α i with i = 1 , 2 , 3 , 4 , 5 and not in terms of the 4 2-dimensional constitutive coefficients A, B, C and D .

In Fig. 1 we represent the scheme of a rectangle with side-names Q, R, S and T and vertex-names V 1 , V 2 , V 3 and V 4 .

Characterization of sides

The characterization of side S is done by setting n i = δ i 1 . Thus, from (15) with α = 1 , 2 , and from ( 16) with α = 1 , 2 , we have

t 1 = t S 1 = u 1 , 1 ( λ + 2 μ) + u 2 , 2 λ -Bu 1 , 111 -2 Du 2 , 222 - 1 2 ( A + B + 2 C ) u 1 , 122 -( B -A ) u 2 , 211 , (17)
t 2 = t S 2 = μ( u 1 , 2 + u 2 , 1 ) -( B -A ) u 1 , 112 -( B -A -2 D ) u 1 , 222 -Au 2 , 111 - 1 2 ( A + B + 2 C ) u 2 , 122 , (18)
τ 1 = τ S 1 = Bu 1 , 11 + 1 2 ( A + B -2 C ) u 1 , 22 + 2 Du 2 , 12 , ( 19
)
τ 2 = τ S 2 = ( B -A -2 D ) u 1 , 12 + Au 2 , 11 + 1 2 ( A + B -2 C ) u 2 , 22 .
(20)

The characterization of side Q is done by setting n i = -δ i 1 . Thus, from ( 15) with α = 1 , 2 , and from ( 16) with α = 1 , 2 , we have

t 1 = t Q 1 = -u 1 , 1 ( λ + 2 μ) -u 2 , 2 λ + Bu 1 , 111 + 2 Du 2 , 222 + 1 2 ( A + B + 2 C ) u 1 , 122 + ( B -A ) u 2 , 211 , (21)
t 2 = t Q 2 = -μ( u 1 , 2 + u 2 , 1 ) + ( B -A ) u 1 , 112 + ( B -A -2 D ) u 1 , 222 + Au 2 , 111 + 1 2 ( A + B + 2 C ) u 2 , 122 , (22)
τ 1 = τ Q 1 = Bu 1 , 11 + 1 2 ( A + B -2 C ) u 1 , 22 + 2 Du 2 , 12 , (23)
τ 2 = τ Q 2 = ( B -A -2 D ) u 1 , 12 + Au 2 , 11 + 1 2 ( A + B -2 C ) u 2 , 22 .
(24)

We remark that t Q 1 in (21) and t Q 2 in ( 22) are the opposite of t S 1 in (17) and of t S 2 in (18) , respectively, and that τ Q 1 in (23) and τ Q 2 in (24) are the same of τ S 1 in ( 19) and of τ S 2 in (20) , respectively.

The characterization of side R is done by setting n i = δ i 2 . Thus, from (15) with α = 1 , 2 , and from ( 16) with α = 1 , 2 , we have

t 1 = t R 1 = μ( u 1 , 2 + u 2 , 1 ) -( B -A ) u 2 , 122 -( B -A -2 D ) u 2 , 111 -Au 1 , 222 - 1 2 ( A + B + 2 C ) u 1 , 112 , (25)
t 2 = t R 2 = u 2 , 2 ( λ + 2 μ) + u 1 , 1 λ -Bu 2 , 222 -2 Du 1 , 111 - 1 2 ( A + B + 2 C ) u 2 , 112 -( B -A ) u 1 , 122 , (26)
τ 1 = τ R 1 = ( B -A -2 D ) u 2 , 12 + Au 1 , 22 + 1 2 ( A + B -2 C ) u 1 , 11 , (27) τ 2 = τ R 2 = B u 2 , 22 + 1 2 ( A + B -2 C ) u 2 , 11 +2D u 1 , 12 .
(28)

We remark that, because of isotropy, t R 1 in (25 Finally, the characterization of side T is done by setting n i = -δ i 2 . Thus, from (15) with α = 1 , 2 , and from ( 16) with α = 1 , 2 , we have

)
t 1 = t T 1 = -μ( u 1 , 2 + u 2 , 1 ) + ( B -A ) u 2 , 122 + ( B -A -2 D ) u 2 , 111 + Au 1 , 222 + 1 2 ( A + B + 2 C ) u 1 , 112 , (29) t 2 = t T 2 = -u 2 , 2 ( λ + 2 μ) -u 1 , 1 λ + Bu 2 , 222 + 2 Du 1 , 111 + 1 2 ( A + B + 2 C ) u 2 , 112 + ( B -A ) u 1 , 122 , (30) τ 1 = τ T 1 = ( B -A -2 D ) u 2 , 12 + Au 1 , 22 + 1 2 ( A + B -2 C ) u 1 , 11 , (31)
τ 2 = τ T 2 = B u 2 , 22 + 1 2 ( A + B -2 C ) u 2 , 11 +2D u 1 , 12 .
(32)

We remark that t T 1 in (29) and t T 2 in (30) are the opposite of t R 1 in (25) and of t R 2 in (26) , respectively, and that τ T 1 in (31) and τ T 2 in (32) are the same of τ R 1 in ( 27) and of τ R 2 in (28) , respectively.

Characterization of vertices

The last equation of ( 8) is reduced, because of (9) 2 to

[ ∂ ∂ B ] δu α f α -f ext α = δu α T αij V i j -f ext α V 1 + δu α T αij V i j -f ext α V 2 + δu α T αij V i j -f ext α V 3 + δu α T αij V i j -f ext α V 4 , ( 33
)
For vertex V 1 the side A has n j = -δ 1 j and ν i = δ i 2 and the side B has n j = δ 2 j and ν i = -δ i 1 so that

V i j V 1 = ν l i n l j + ν r i n r j V 1 = -δ i 2 δ 1 j -δ i 1 δ 2 j .
For vertex V 2 the side B has n j = δ 2 j and ν i = δ i 1 and the side C has n j = δ 1 j and ν i = δ i 2 so that

V i j V 2 = ν l i n l j + ν r i n r j V 2 = δ i 1 δ 2 j + δ i 2 δ 1 j .
For vertex V 3 the side C has n j = δ 1 j and ν i = -δ i 2 and the side D has n j = -δ 2 j and ν i = δ i 1 so that

V i j V 3 = ν l i n l j + ν r i n r j V 3 = -δ i 2 δ 1 j -δ i 1 δ 2 j .
For vertex V 4 the side D has n j = -δ 2 j and ν i = -δ i 1 and the side A has n j = -δ 1 j and ν i = -δ i 2 so that

V i j V 4 = ν l i n l j + ν r i n r j V 4 = δ i 1 δ 2 j + δ i 2 δ 1 j .
Thus, finally, the (33) yields

[ ∂ ∂ B ] δu α f α -f ext α = δu α -T α21 -T α12 -f ext α V 1 + δu α T α12 + T α21 -f ext α V 2 + δu α -T α21 -T α12 -f ext α V 3 + δu α T α12 + T α21 -f ext α V 4 , ( 34
)
where T α12 + T α21 , in terms of the displacement field, it is for α = 1

T 112 + T 121 = 2 Cu 1 , 12 + ( B -A -2 D ) u 2 , 11 + 2 Du 2 , 22 , (35)
and for α = 2 ,

T 212 + T 221 = 2 Cu 2 , 12 + ( B -A -2 D ) u 1 , 22 + 2 Du 1 , 11 .
(36)

Numerical simulations

Numerical data for the simulations that will be presented in this paper are here shown (see Fig. 1 )

L = 2 m, l = 1 m, μ = 10 MP a m, λ = 15 MP a m, ρ = 10 5 Kg/m 2 E = μ( 3 λ + 2 μ) λ + μ = 26 MP a m, ( 37
)
α 1 = El 2 m , α 2 = El 2 m , α 3 = 2 El 2 m , α 4 = El 2 m , α 5 = 1 2 El 2 m , l m = 10 cm, ( 38
)
and therefore

A = 7 El 2 m , B = 34 El 2 m , C = 21 2 El 2 m , D = 8 El 2 m .
With these data the positive definiteness of the strain energy functional is guaranteed.

The heavy sheet problem

We consider an heavy sheet appended at the top side R and constrained at sides Q and S to have null horizontal displacement. Thus, the kinematical restrictions are ( δu 2 ) R = 0 , ( δu 1 ) Q = 0 , ( δu 1 ) S = 0 .

(39) and also represented in Fig. 2 . Let the two partial differential Eqs. ( 6) and ( 7) be satisfied with the following external forces per unit area,

b ext 1 = 0 , b ext 2 = -ρg, (40)
and let the edge boundary conditions be as follows,

u 1 = 0 , t Q 2 = 0 , τ Q 1 = 0 , τ Q 2 = τ ext,Q 2 , ∀ X ∈ Q, ( 41
) t R 1 = 0 , u 2 = 0 , τ R 1 = 0 , τ R 2 = τ R,ext 2 , ∀ X ∈ R, ( 42
)
u 1 = 0 , t S 2 = 0 , τ S 1 = 0 , τ S 2 = τ ext,S 2 , ∀ X ∈ S, (43)
t T 1 = 0 , t T 2 = 0 , τ T 1 = 0 , τ T 2 = τ T,ext 2 , ∀ X ∈ T , ( 44
)
where the non-null external edge double forces are as follows,

τ ext,Q 2 = τ ext,S 2 = ( A + B -2 C ) ρg 2 ( λ + 2 μ) , τ R,ext 2 = τ T,ext 2 = ρgB ( λ + 2 μ) , ( 45
)
and let the only wedge conditions that are not implied by ( 41) -(44) be as follows,

f ext 2 V 3 = 0 , f ext 2 V 4 = 0 . ( 46
)
The analytical solution of this problem is achieved in [START_REF] Luca | Gedanken experiments for the determination of two-dimensional linear second gradient elasticity coefficients[END_REF] and here represented in terms of the displacement field,

u 1 = 0 , u 2 = ρg ( X 2 -l ) ( 3 l + X 2 ) 2 ( λ + 2 μ) . ( 47
)
The numerical simulations of this problem is shown in Fig. 3 a and have shown remarkable identification between exact analytical solution and the respective numerical simulation ( Fig. 3 b)

The bending problem

We consider the bending problem and constrain the whole side A to not displace in the horizontal direction and one of its point, the origin O , to have also null vertical displacement. Thus, the kinematical restrictions are ( δu 1 ) A = 0 , ( δu 2 ) O = 0 .

(48) and they are represented in Fig. 4 . It has to be remarked that such kinematical constrains are not of a general type. In fact, the (48) 2 is referred to a single point that is not a vertex of the domain. This means that the results will be reasonable only in the case of vertical null force at point O . Let the two partial differential Eqs. ( 6) and ( 7) be satisfied with the null external forces per unit area,

b ext 1 = 0 , b ext 2 = 0 , ( 49
)
and the edge boundary conditions be as follows,

1 = 0 , t Q 2 = 0 , τ Q 1 = 0 , τ Q 2 = τ ext,Q 2 , ∀ X ∈ Q, ( 50
) t R 1 = 0 , t R 2 = 0 , τ R 1 = 0 , τ R 2 = τ R,ext 2 , ∀ X ∈ R, ( 51
)
t 1 = t ext,S 1 , t S 2 = 0 , τ S 1 = 0 , τ S 2 = τ ext,S 2 , ∀ X ∈ S, ( 52
) t T 1 = 0 , t T 2 = 0 , τ T 1 = 0 , τ T 2 = τ T,ext 2 , ∀ X ∈ T , ( 53
)
where the non-null external edge force and double forces are as follows,

t ext,S 1 = 3 M ext X 2 2 l 3 , (54) τ ext,Q 2 = τ ext,S 2 = 3 M ext [ -( 5 λ + 8 μ) A + ( λ + 4 μ) B + 2 λC -( 4 λ + 8 μ) D ] 16 l 3 μ( λ + μ) , ( 55
) τ R,ext 2 = τ T,ext 2 = - 3 M ext [ ( λ + 2 μ) A + ( 3 λ + 2 μ) B -( 2 λ + 4 μ) C -( 4 λ + 8 μ) D ] 16 l 3 μ( λ + μ) . ( 56
)
and let the only wedge conditions that are not implied by ( 50) -( 53) be as follows,

f ext 2 V 1 = 0 , f ext 2 V 2 = 0 , f ext 2 V 3 = 0 , f ext 2 V 4 = 0 , f ext 1 V 2 = -f ext 1 V 3 = 3 M ext [ ( λ + 2 μ) ( A -B + 2 C ) + 4 μD ] 8 l 3 μ( λ + μ)
.

The analytical solution of this problem is achieved in [START_REF] Luca | Gedanken experiments for the determination of two-dimensional linear second gradient elasticity coefficients[END_REF] and here represented in terms of the displacement field, The numerical simulations of this problem are shown in Figs. 5 . The comparisons between the exact analytical solution and the numerical simulation have been done through the horizontal and vertical cuts of Fig. 6 and the respective numerical simulation is shown in Figs. 7 and8 .

u 1 = 3 M ext ( λ + 2 μ) X 1 X 2 8 l 3 μ( λ + μ) , u 2 = - 3 M ext λX 2 2 + ( λ + 2 μ) X 2 1 16 l 3 μ( λ + μ) , ( 57
)

The bending problem without double forces

We consider the same bending problem of the previous subsection, with the same kinematical restrictions (48) , the same external forces per unit area (49) , the same edge boundary conditions ( 50) -( 53) with the same edge force (54) but with null edge double forces

τ ext,Q 2 = τ ext,S 2 = 0 , τ R,ext 2 = τ T,ext 2 = 0 ,
and null wedge conditions

f ext 2 V 1 = 0 , f ext 2 V 2 = 0 , f ext 2 V 3 = 0 , f ext 2 V 4 = 0 , f ext 1 V 2 = -f ext 1 V 3 = 0 .
In this case we do not have an analytical solution but we make numerical simulations, that are shown in Figs. 9 and 10 . It can be remarked that the presence of double forces has a relatively strong influence on the numerical results.

The flexure problem

We consider the flexure problem and constrain the whole side C to displace in the vertical direction and one of the point of side A , the origin O , to have null horizontal and vertical displacements. Thus, the kinematical restrictions are ( δu 1 ) S = 0 , ( δu 1 ) O = 0 ( δu 2 ) O = 0 .

(58) C as follows,

t ext,S 1 = 0 , τ ext,Q 2 = τ ext,S 2 = 0 , τ R,ext 2 = τ T,ext 2 = 0 ,
and with the following system of wedge conditions

f ext 2 V 1 = 0 , f ext 2 V 2 = 0 , f ext 2 V 3 = -F , f ext 2 V 4 = 0 , f ext 1 V 2 = -f ext 1 V 3 = 0 .
Even in this modified bending case we do not have an analytical solution but we make numerical simulations, that are shown in Figs 15 . On the left-hand side we show the results of a numerical simulation of the analogous problem in the classical first gradient model and on the right-hand side a the results of a numerical simulation with the present second gradient model. It is immediately visible from Fig. 15 that the first gradient model is not adequate for wedge concentrated external forces. In Fig. 16 a convergence analysis is performed. From such a convergence analysis we deduce that first gradient models are not adequate to model concentrated external forces.

Experimental evidence of elastic second gradient contribution to the deformation energy

In [START_REF] Ferretti | Modelling the onset of shear boundary layers in fibrous composite reinforcements by second gradient theory[END_REF] it is shown that second-gradient energy terms allow the onset of internal shear boundary layers. These boundary layers are transition zones between two different shear deformation modes. In the same paper, on the one hand it is claimed that their existence cannot be described by a simple first-gradient model, and on the other hand that they are related to second-gradient material coefficients. In this section we show a result, on the pantographic structure of Fig. 17 , that makes explicit the experimental evidence of such a boundary layer. Besides, by using the second gradient model of dell 'Isola et al. (2016) , we show that it is possible to characterize the largeness of the boundary layer in terms of the second gradient coefficient (i.e., K II ) of that model, that means that a simple first gradient model (i.e., a model with K II = 0 ) is not sufficient to predict the correct experimental evidence.

U ( F , ∇F ) = 2 α=1 K e 2 F ab D α b -1 ( F ac D α c -1 ) + K II 2 F ab,c F ad,e D α b D α c D α d D α e F f g F f h D α g D α h - F ab F ac,d D α b D α c D α d F f g F f h D α g D α h 2 + K b 2 ⎡ ⎣ arccos F ab F ac D 1 b D 2 c F de F df D 1 e D 1 f F gh F gi D 2 h D 2 i - π 2 ⎤ ⎦ γ ( 73
)
where the two families of fibers are initially directed along the two orthogonal unit vectors D 1 and D 2 and where the material coefficients that have been used are K e = 0 . 134 MN/m, K p = 159 N/m and γ = 1 . 36 . We also remark that the bias test that is shown in Fig. 17 is accomplished by imposing a displacement, towards the direction parallel to the long side of the rectangle, of the short-side of the rectangle, that in turn is directed at π /4 with respect to the two orthogonal unit vectors D 1 and D 2 . In the deformed configuration the angle ϕ between the two families of fibers is not anymore at π /2. It is ϕ = arccos

F ab F ac D 1 b D 2 c F de F df D 1 e D 1 f F gh F gi D 2 h D 2 i ,
and the shear angle φ (or shear deformation in Fig. 18 ) is simply φ = ϕπ

2 .

If we evaluate the shear angle φ along the arc-length that is shown in Fig. 17 , it is almost zero near short side of the rectangle and reach a finite value by passing through a boundary layer. The largeness of such a boundary layer is related to the second gradient parameter K II .

In Fig. 18 we show that the optimal value for the second gradient constitutive parameter K II is K II = 0 . 0192 Nm . However, we also shown that different values of this parameter give a wrong largeness of the boundary layer. In particular we observe that a reduction 1/4 of the second gradient constitutive parameter K II give a smaller boundary layer and a magnification of 4 give a larger boundary layer. Finally, a numerical simulation in which the second gradient contribution of the strain energy (73) is assumed to vanish, i.e. K II = 0 , produce no-boundary layer. Besides, numerically instability is observed in this last case.

Conclusions

A perfect overlap of numerical simulations obtained with a commercial code and closed form solution of selected classical benchmark boundary value problems have been found and reported in this paper. The role of external double and wedge forces has also been presented. Besides, we show a mesh-independent behaviour of second gradient numerical solution with respect to the corresponding first gradient counterpart. Finally, we show an experimental bias test on a specific pantographic structure and extrapolate an internal boundary layer in terms of the shear angles across initially orthogonal fibers. A non-linear anisotropic model is also presented aimed to reproduce the shown experimental results. In particular, we exhibit comparisons between the numerical simulation of the proposed theoretical model and the experimental results in terms of the internal boundary layer. Such a comparison has permitted to identify the second gradient coefficient of this model.

  where b ext α , t ext α , τ ext α and f ext α are the external actions: b ext α is the external force per unit area and is applied on the whole 2-dimensional domain B; t ext α and τ ext α are the external force and double force (respectively) and are applied on (a part of) the one-dimensional boundary ∂B of the domain B; and f ext α is the external concentrated force applied on the set of points belonging to the boundary of the boundary [ ∂ ∂ B ] , so that the last integral can be also represented as the sum of the external works made by the concentrated forces acting on each vertex of the domain. In other words, if we define the boundary ∂B as the union of m regular parts c with c = 1 , . . . , m and [ ∂ ∂ B ] as the union of the corresponding m vertex-points V c with c = 1 , . . . , m,

Fig. 1 .

 1 Fig. 1. Picture of the 2-dimensional body B.

Fig. 2 .

 2 Fig. 2. Graphical representation of the kinematic restrictions for the heavy sheet problem.

Fig. 3 .

 3 Fig. 3. Vertical displacement of the heavy sheet problem. (a) Graphical representation of the numerical simulation and (b) comparison between the exact analytical solution and the respective numerical one through a given vertical cut.

Fig. 4 .

 4 Fig. 4. Graphical representation of the kinematic restrictions for the bending problem.

  u

Fig. 5 .Fig. 6 .

 56 Fig. 5. Numerical simulation for the bending problem. In particular, (a) horizontal and (b) vertical displacements are represented.

Fig. 7 .Fig. 8 .

 78 Fig. 7. Comparison between the exact analytical solution and numerical simulation through the horizontal cuts of Fig. 6 . In particular, (a) horizontal and (b) vertical displacements are represented.

Fig. 9 .

 9 Fig. 9. Comparison between the exact analytical solution with non null external double and wedge forces and numerical simulation with null external double and wedge forces through the horizontal cuts of Fig. 6 . In particular, (a) horizontal and (b) vertical displacements are represented.

Fig. 14 .

 14 Fig. 14. Comparison between the exact analytical solution and numerical simulation through the vertical cuts of Fig. 6 . In particular, (a) horizontal and (b) vertical displacements are represented.

Fig. 15 .

 15 Fig. 15. Comparison between the vertical displacement in the (a) first and (b) second gradient models in the case of an external wedge force.

Fig. 16 .

 16 Fig. 16. Comparison between the vertical displacement at the center of side S in the (a) first and (b) second gradient models in the case of an external vertical wedge force at V 3 .

Fig. 18 .

 18 Fig. 18. The boundary layer of the case related to the experiment that is shown in Fig. 17 is shown. The angles across the two families of fibers in the deformed configuration are evaluated by image analysis of the experimental result of Fig. 17 and by numerical simulations with different values of second gradient coefficients K II .

and they are represented in Fig. 11 . It has to be remarked that such kinematical constrain are not of a general type. In fact, the (58) 2,3 is referred to a single point that is not a vertex of the domain. This means that the results will be reasonable only in the case of null force at point O . Let the two partial differential Eqs. ( 6) and ( 7) be satisfied with null external forces per unit area (49) , let the edge boundary conditions be as follows,

where the non-null external edge force and double forces are as follows, and the only wedge conditions that are not implied by ( 59) -( 62) are as follows,

analytical solution of this problem is achieved in [START_REF] Luca | Gedanken experiments for the determination of two-dimensional linear second gradient elasticity coefficients[END_REF] and here represented in terms of the displacement field,

The numerical simulations of this problem are shown in Figs. 12 . The comparisons between the exact analytical solution and the numerical simulation have been done through the horizontal and vertical cuts of Fig. 6 and the respective numerical simulation is shown in Figs. 13 and 14 .

Convergence analysis for the wedge force problem

We consider again the bending problem with the same kinematical restrictions (48) (see also Fig. 4 ), the same external forces per unit area (49) , the same edge boundary conditions ( 50) -( 53) but with null edge forces and double forces at side