Juan Cortés
email: jcortes@laas.fr

Sophie Barbe
email: sbarbe@insa-toulouse.fr

Monique Erard
email: monique.erard@ipbs.fr

Thierry Siméon

Encoding Molecular Motions in Voxel Maps

Understanding life at the atomic level requires the development of new methodologies, able to overcome the limitations of available experimental and computational techniques for the analysis of processes involving molecular motions. With this goal in mind, we develop new methods, combining robotic path planning algorithms and molecular modeling techniques, for computing large-amplitude motions. This paper builds on these new methods, and introduces voxel maps as a computational tool to encode and to represent such motions. Voxel maps can be used to represent relative motions of two molecules, as well as conformational changes in macromolecules. We investigate several applications and show results that illustrate the interest of such representation. In particular, voxel maps are used to display channels into proteins, to analyze protein-ligand specificity, and to represents protein loop and domain motions.

I. INTRODUCTION

Nowadays, the dynamic nature of biological macromolecules, opposed to the static picture provided by X-ray crystallography, is generally accepted. Furthermore, it has been shown that flexibility plays key roles in molecular interactions such as protein-ligand [START_REF] Carlson | Protein flexibility is an important component of structurebased drug discovery[END_REF], [START_REF] Cavasotto | The challenge of considering receptor flexibility in ligand docking and virtual screening[END_REF] and proteinprotein docking [START_REF] Janin | Assessing predictions of protein-protein interaction: the CAPRI experiment[END_REF], [START_REF] Ehrlich | The impact of protein flexibility on protein-protein docking[END_REF]. Unfortunately, and despite great advances achieved in the last years [START_REF] Katona | Raman-assisted crystallography reveals end-on peroxide intermediates in a nonheme iron enzyme[END_REF], [START_REF] Schanda | Protein folding and unfolding studied at atomic resolution by fast two-dimensional nmr spectroscopy[END_REF], an atomicresolution structural description of slow-timescale (largeamplitude) molecular motions is out of reach for currently available experimental methods [START_REF] Henzler-Wildman | Dynamic personalities of proteins[END_REF]. Computational methods are therefore necessary to complement experimentation.

Molecular dynamics (MD) [START_REF] Leach | Molecular Modeling: Principles and Applications[END_REF], [START_REF] Schlick | Molecular Modeling and Simulation -An Interdisciplinary Guide[END_REF] is the most widely used computational method to simulate molecular motions. MD is an appropriate method to analyze motions taking place in a short timescale (up to some nanoseconds). However, it is too computationally expensive for routine simulations of large-amplitude motions of macromolecules. In some cases, MD simulations can be accelerated by the introduction of artificial forces [START_REF] Izrailev | Steered molecular dynamics[END_REF]. Nevertheless, devising such forces may require prior knowledge of the particular problem, and they can excessively bias the resulting trajectories. Simulation methods based on Monte Carlo (MC) algorithms [START_REF] Leach | Molecular Modeling: Principles and Applications[END_REF], [START_REF] Schlick | Molecular Modeling and Simulation -An Interdisciplinary Guide[END_REF] have been developed to overcome the limitations of MD. Such methods present however a major drawback for computing large-amplitude motions since the conformational exploration tends to get trapped into the many local minima of the complex molecular energy landscape. Alternatives to MD and MC simulations have been proposed using very different methods such as iterative NMA calculations [START_REF] Mouawad | Motions in hemoglobin studied by normal mode analysis and energy minimization: evidence for the existence of tertiary T-like, quaternary R-like intermediate structures[END_REF], [START_REF] Jeong | A method for finding candidate conformations for molecular replacement using relative rotation between domains of a known structure[END_REF], or structurally constrained conformational exploration using models from rigidity theory [START_REF] Wells | Constrained geometric simulation of diffusive motion in proteins[END_REF].

Our method builds on path planning algorithms [START_REF] Choset | Principles of Robot Motion: Theory, Algorithms, and Implementations[END_REF], [START_REF] Lavalle | Planning Algorithms[END_REF], originally developed in the field of robotics. Such algorithms are efficient tools for exploring constrained high-dimensional spaces. In the recent years, path-planning-based methods have been successfully applied for investigating different problems in structural biology such as: protein-ligand access and docking [START_REF] Singh | A motion planning approach to flexible ligand binding[END_REF], [START_REF] Apaydin | Stochastic conformational roadmaps for computing ensemble properties of molecular motion[END_REF], [START_REF] Cortés | A path planning approach for computing large-amplitude motions of flexible molecules[END_REF], protein and RNA folding [START_REF] Amato | Using motion planning to map protein folding landscapes and analyze folding kinetics of known native structures[END_REF], [START_REF] Apaydin | Stochastic roadmap simulation: an efficient representation and algorithm for analyzing molecular motion[END_REF], [START_REF] Tang | Using motion planning to study RNA folding kinetics[END_REF], protein loop motions [START_REF] Cortés | Geometric algorithms for the conformational analysis of long protein loops[END_REF], domain motions [START_REF] Kirillova | An NMA-guided path planning approach for computing large-amplitude conformational changes in proteins[END_REF], and motions of pairs of α-helices in transmembrane proteins [START_REF] Enosh | Prediction and simulation of motion in pairs of transmembrane α-helices[END_REF].

This paper recalls our approach for computing molecular motions (Section II-B) and introduces voxel maps as a new and general computational tool to encode them (Section II-C). Section III illustrates the potential interest of the method on several structural biology problems. The presented results show how voxel maps can effectively represent relative motions of two molecules, as well as conformational changes in proteins. In Section III-A, voxel maps are applied to identify channels in proteins, by exploring and encoding possible motions of a single atom between the active site and the surface. The second application (Section III-B) addresses protein-ligand interactions. Voxel maps permit to reflect differences between the access/exit pathways of different ligands to the active site of a protein. Finally, Section III-C deals with the representation of conformational changes due to loop and domain motions.

II. METHODOLOGY A. Overview

The method presented in this paper builds on the two-stage approach proposed in [START_REF] Cortés | A path planning approach for computing large-amplitude motions of flexible molecules[END_REF] for computing large-amplitude molecular motions. The first and main stage consists in a geometric processing of the strongest molecular constraints (no atom overlaps, no bond breaking). Fast geometric computation [START_REF] Ruiz De Angulo | BioCD: An efficient algorithm for self-collision and distance computation between highly articulated molecular models[END_REF] combined with efficient path planning algorithms [START_REF] Cortés | Disassembly path planning for complex articulated objects[END_REF] permits our method to generate large-amplitude motions of flexible molecules with very low computational cost. In the second stage, results of the geometric exploration can be refined and analyzed using classic molecular modeling tools (e.g. energy evaluation/minimization).

The voxel-map representation described below can be seen as an intermediate layer between the two stages. It permits to arrange the information obtained from the exploration of a high-dimensional space (the molecular conformational space) into a simple three-dimensional data structure. The choice of the three dimensions and the size of voxels depends on the application (see Section III). In addition to the information structuring, voxel maps permit a visual analysis of the results of the conformational exploration.

B. Exploring geometrically feasible motions

The conformational search method applied in this work (described in more detail in [START_REF] Cortés | A path planning approach for computing large-amplitude motions of flexible molecules[END_REF]) is based on a mechanistic modeling of molecules [START_REF] Zhang | A new method for fast and accurate derivation of molecular conformations[END_REF]. Groups of bonded atoms form the bodies of the mechanism, which are linked by articulations corresponding to bond torsions. These torsions are the molecular degrees of freedom. The atoms are represented by rigid spheres with (a percentage of) van der Waals radii. These spheres cannot overlap. Additional distance and orientation constraints can be imposed between elements of this mechanistic model in order to simulate interactions such as hydrogen bonds.

The technique applied to explore feasible motions of the mechanistic molecular model is derived from the Rapidlyexploring Random Trees (RRT) algorithm [START_REF] Lavalle | Rapidly-exploring random trees: Progress and prospects[END_REF]. The basic principle of RRT is to incrementally grow a random tree, rooted at a given initial conformation q init , to explore the search space for finding feasible paths. The exploration process is illustrated in Figure 1-b. At each iteration, the tree is expanded toward a randomly sampled conformation q rand . This random sample is used to simultaneously determine the tree node to be expanded and the direction in which it is expanded. Given a distance metric in the search space, the nearest node q near in the tree is selected. The expansion of q near in the direction of q rand leads to the generation of a a) b) d) c) q q q q init near new rand Fig. 1. Two-dimensional illustration of a voxel map construction. Given a geometric model of the molecules (a), a path planning algorithm is used to explore the subset of the conformational space that is reachable from the initial conformation satisfying motion constraints. b) Illustration of an expansion step of an RRT-like path planning algorithm. c) Search tree resulting from the geometric exploration. d) Voxel map associated with the search tree.

new node q new and a feasible local path. The key feature of the RRT expansion strategy is to bias the exploration toward unexplored regions before uniformly covering the space.

The algorithmic variant used in this work is the ML-RRT [START_REF] Cortés | Disassembly path planning for complex articulated objects[END_REF]. This variant considers two sets of conformational parameters: active and passive. Active parameters are randomly sampled, while passive parameters are treated separately, and only when they hinder the motion of active parts. The main advantage of ML-RRT is a higher efficiency for dealing with high-dimensional problems thanks to the decoupled treatment of parameter subsets, which favors the exploration of active parameters. In the applications presented in this paper, the passive parameters correspond to the torsion angles of the protein side-chains, while the active parameters correspond to the other variables: the location and the internal torsions of the ligand, and the torsion angles of flexible protein backbone segments.

C. Putting search trees into voxel maps

Nodes and edges of the RRT search tree are embedded in a high-dimensional space (the conformational space of the molecular model). The idea is to arrange this information into a lower-dimensional data structure for facilitating further analysis. The voxel-map representation has been chosen for different reasons: (1) it is a simple and regular structure, which facilitates operations such as nearest-neighbor search, (2) it is three-dimensional, which permits a visual display.

The three dimensions of the voxel map correspond to three variables of interest for the particular problem. These three variables may be a subset of the conformational parameters of the mechanical model (e.g. position of the reference frame of a molecule, three selected bond torsions). They can also be chosen to encode information obtained from the conformation (e.g. position of one atom, center of mass of one domain). The voxel size also depends on the application. For instance, the resolution may be chosen depending on the motion amplitude, and on the cost of the ulterior treatment (geometric and/or energetic analysis). Note that increasing the voxel resolution does not affect computational efficiency.

The process for generating the voxel map from the RRT search tree is very simple. It is illustrated in Figure 1-d. Each time a node q new is generated, it is added to the corresponding voxel. The procedure is slightly more complex for the edges. When an edge goes across voxels that do not contain the nodes it connects, an intermediate conformation along the edge is kept into each new voxel. Thus, each voxel contains a list of conformations corresponding to the RRT nodes and edges that are projected on it.

Voxels can be labeled in different ways. For instance, values can be assigned depending on the chronological order of generation. In this way, the voxel map can be used to display the regions of the space that are reached first during the conformational exploration (see Section III-B). Other labeling procedures can be devised based on conformational or geometric features, such as the distance between catalytic residues (see the protein loop example in Section III-C.1). The voxel labeling can also be made according to energy evaluation. An energetic analysis of the conformations associated with voxels may provide very useful information about the conformational energy landscape (see the example of protein domain motions in Section III-C.2).

III. APPLICATIONS

This section presents several examples illustrating voxel map computations for different types of problems involving molecular motions. The goal is to show the interest of the proposed approach. Therefore, the section does not provide in-depth explanations about the considered systems and the biological importance of the results.

The method has been implemented within our software prototype BioMove3D. PyMOL [START_REF] Delano | The PyMOL molecular graphics system[END_REF] has been used for viewing molecular models and voxel maps. The computing times given below correspond to tests run on a single AMD Opteron 148 processor at 2.6 GHz.

A. Findings channels into proteins

The most straightforward application of the method is the search and representation of channels into proteins. The channels are searched using the ML-RRT algorithm to explore feasible motions of a ball of arbitrary radius inside the protein model. This algorithm permits to directly treat all side-chain flexibility with a low computational cost. In this type of application, the three variables used for the voxelmap representation are the position parameters of the moving ball. The voxel resolution can be chosen in relation to the ball size.

The chosen benchmark for this application is cytochrome P450. The in/out channels of this enzyme have been recently characterized in [START_REF] Cojocaru | The ins and outs of cytochrome P450s[END_REF]. The authors of the cited study applied the technique CAVER [START_REF] Petřek | CAVER: A new tool to explore routes from protein clefts, pockets and cavities[END_REF] for computing channels between the buried active site and the surface in different crystal structures of P450. CAVER is based on the construction of a vertex-weighted graph from a discrete three-dimensional grid model of the protein. The weights are computed from the distance to the protein atoms, the lowest weights corresponding to nodes with the highest clearance. A variant of Dijkstra's algorithm is applied to search the shortest lowcost paths. CAVER considers static structures, and performs systematic exploration of the protein interior. Molecular flexibility can only be indirectly treated by applying the technique to a set of structures (e.g. samples of a molecular dynamics simulation). The results described below show several advantages of the voxel map method over CAVER.

The structure represented in Figure 2 corresponds to bacterial P450-BM3 (PDBid 1JPZ). The voxel map in the figure represents the channels found by our method for a moving ball of radius 1.2 Å. Three of these channels, W, 2b and 2f, were also found by CAVER. However, channel 2d was not reported for this structure. Nevertheless, this channel was observed in a small number of other P450-BM3 structures [START_REF] Cojocaru | The ins and outs of cytochrome P450s[END_REF]. A further analysis of channel 2d shows that the exit of the moving ball requires slight motions of some side-chains, in particular those of residues Leu20 and Leu29. This result shows the interest of considering sidechain motions when computing channels in proteins.

Also note that very similar channels were obtained from several runs of the voxel map computation with the moving ball initially located at different positions in the enzyme active site. Such a reliability shows the low sensibility to the starting point, which is another advantage over CAVER 1 .

Finally note that, even considering side-chain flexibility, the computing time required to construct the voxel map representation of the channels remains low (only 5 minutes).

B. Analyzing ligand access/exit pathways

The proposed approach can be applied to carry out an analysis of access/exit pathways of ligands (or substrates/products) to the active site of a protein. Such pathways can play important roles in the protein activity and specificity, particularly for proteins presenting a deep binding pocket. For instance, they may be one factor in enzyme enantioselectivity [START_REF] Guieysse | A structurecontrolled lipase enantioselectivity investigated by a path planning approach[END_REF].

The example illustrated in Figure 3 corresponds to the interaction between Candida antarctica lipase B (CALB) and (R,S)-enantiomers of 4-methyloctanoic acid [START_REF] Franssen | The interaction between candida antarctica lipase and branched chain fatty acids: a kinetic and molecular modelling study[END_REF]. Our method has been applied to explore possible motions of the R and S substrates from its catalytic position. The model considers the flexibility of the substrate and all protein sidechains. The voxel maps computed in this case represent the feasible positions reached by the substrate center of mass during the conformational exploration. Voxel resolution is 0.1 Å. This representation shows significant differences between both enantiomers. First, the bottom of the voxel map is narrower for the (S)-enantiomer than for the (R)-enantiomer, which shows that possible motions of the (S)-enantiomer are more constrained in this region. Voxels have been colored depending on the chronological order of generation. Darkblue voxels correspond to the positions that are reached first. For the (R)-enantiomer (Figure 3-c), dark-blue voxels reach the middle-part of the map. However, such voxels are concentrated in the bottom part of the map for the (S)enantiomer (Figure 3-d). This means that, only considering geometric constraints, the access/exit of the (R)-enantiomer can be faster. These results tend to indicate that the topology of CALB active site is better suited for facilitating the reaction with (R)-4-methyloctanoic acid, which correlates with experimental results on the enzyme enantioselectivity showing a significant preference for this enantiomer [START_REF] Franssen | The interaction between candida antarctica lipase and branched chain fatty acids: a kinetic and molecular modelling study[END_REF].

Concerning computational performance, the construction of the voxel map (including the ML-RRT exploration process) took less than 7 minutes for each enantiomer. Note that, aiming to get a very good coverage of the space into the active site cavity, the ML-RRT search tree expansion process yielded voxel maps with more than 25,000 and 15,000 voxels respectively for the (R)-and (S) enantiomers. This result shows that even such an exhaustive exploration is computationally fast.

C. Representing loop/domain motions

Our path-planning-based approach can be applied to compute large-amplitude internal molecular motions such as loop and domain motions [START_REF] Cortés | Geometric algorithms for the conformational analysis of long protein loops[END_REF], [START_REF] Cortés | A path planning approach for computing large-amplitude motions of flexible molecules[END_REF], [START_REF] Kirillova | An NMA-guided path planning approach for computing large-amplitude conformational changes in proteins[END_REF]. Integrating the voxelmap representation in this approach provides a new tool for analyzing such deformations.

1) Loop motions: The example illustrated in Figure 4 concerns the "WDP loop" in Yersinia protein tyrosine phosphatase (PTPase). The movement of the WDP loop plays a central role in the PTPase-mediated catalytic process [START_REF] Stuckey | Crystal structure of Yersinia protein tyrosine phosphatase at 2.5 Å and the complex with tungstate[END_REF], [START_REF] Hu | Dynamics of the WPD loop of the Yersinia protein tyrosine phosphatase[END_REF]. An open conformation of this loop permits the substrate access to the protein active site. Then, the WDP loop is required to adopt a closed conformation that brings the catalytic residue Asp356 to a specific location for proteinsubstrate interaction. Starting from the open conformation of PTPase [START_REF] Stuckey | Crystal structure of Yersinia protein tyrosine phosphatase at 2.5 Å and the complex with tungstate[END_REF] (PDB ID: 1YPT), the ML-RRT algorithm was applied to explore the mobility of the WDP loop (residues 352-361). A voxel map obtained from this exploration is represented in the right part of Figure 4. It displays the positions reached by the Cα atom of the middle loop residue Glu357. Voxels size is 0.5 Å, and colors have been assigned depending on the distance between the referred atom and the Cα of Val407, which is located on the bottom of the binding pocket. These atoms were also chosen in [START_REF] Hu | Dynamics of the WPD loop of the Yersinia protein tyrosine phosphatase[END_REF] to measure the WDP loop gating during molecular dynamics simulations. The distance in the initial crystal structure is 17 Å. The minimum and maximum distances obtained through the conformational exploration are 11 Å and 22 Å respectively. The WDP loop reaches conformations very similar to the one in the closed structure [START_REF] Stuckey | Crystal structure of Yersinia protein tyrosine phosphatase at 2.5 Å and the complex with tungstate[END_REF] (PDB ID: 1YTN), with Cα RMSD below 1 Å. The voxel map shows that the loop can adopt more open conformations, as also suggested by molecular dynamic simulations [START_REF] Hu | Dynamics of the WPD loop of the Yersinia protein tyrosine phosphatase[END_REF]. Interestingly, the voxel map presents a marked pipe-like shape. Such a shape indicates that the WDP loop in PTPase is "mechanically" designed to perform opening-closure motions, while lateral-motions are not likely. This mechanical predisposition may explain the rapid opening-closure WDP loop motions reported in [START_REF] Hu | Dynamics of the WPD loop of the Yersinia protein tyrosine phosphatase[END_REF]. Finally note that constructing the voxel map required about 5000 iterations of the ML-RRT algorithm, with a computing time of 5 minutes.

2) Domain motions: Similarly, the proposed approach can be applied for analyzing protein domain motions. The example presented here concerns the N-Oct-3 POU domain. This molecule, represented in Figure 5, comprises two distinct, highly conserved sub-domains, termed POUs and POUh, connected by a flexible linker [START_REF] Latchman | POU family transcription factors in the nervous system[END_REF]. It has been shown that, due to its remarkable plasticity, the N-Oct-3 POU domain can adopt different configurations and corresponding homodimerization patterns, with the POUs and POUh sub-domains acting as sensors for the distinct underlying structures which characterize the respective DNA targets [START_REF] Alazard | Identification of the 'NORE' (N-Oct-3 responsive element), a novel structural motif and composite element[END_REF]. The conformational transitions of the N-Oct-3 POU domain are currently being investigated by combining experimental techniques (SAXS) and the computational approach presented in this paper. Figure 6 shows preliminary results of this work. The voxel map represents possible locations of POUh with respect to POUs (relative positions of the centers of mass) considering the flexibility of the long linker between them (18 residues are considered to be fully flexible and 13 have limited flexibility). All the protein sidechains are also potentially flexible (i.e. they move if they hinder backbone motions). Maximum and minimum values of the radius of gyration have been integrated as constraints for the conformational exploration. This structural information has been derived from SAXS experiments coupled to molecular mechanics [START_REF] Alazard | Fine-tuning of intrinsic N-Oct-3 POU domain allostery by regulatory DNA targets[END_REF]. The voxels size is 2 Å. The wideness of the explored region reflects the high flexibility of the linker. Voxels have been labeled using an energetic scoring. The conformations associated with each voxel have been clustered into significantly different sets. Then, one conformation of each cluster has been energy minimized 2 imposing constraints on the backbone atom positions in order to remain within the voxel. The voxel color is assigned according to the resulting lowest-energy conformation. The goal of this work is to compute a complete conformational map of the protein enabling the identification of low-energy conformations attainable from the initial structure, and possible transition pathways between them. Comparison with free and pre-bound models of N-Oct-3 POU constructed from experimental data and molecular modeling tools will allow us to validate this novel approach.

IV. CONCLUSIONS

We have presented a general framework for computing and representing molecular motions. The basic principle of the approach is to apply path-planning-based algorithms to explore feasible motions of mechanistic molecular models. The efficiency of such conformational exploration permits to attain large-amplitude motions with low computational cost. The voxel-map representation facilitates ulterior treatment of the explored conformations and permits direct visual interpretation of results. Besides, voxel maps could be used to bias or to focus the exploration to specific regions of 2 AMBER ff03 force field [START_REF] Duan | A point-charge force field for molecular mechanics simulations of proteins[END_REF] has been used for the energetic analysis. the conformational space, and could permit to device more accurate metrics, considering motion feasibility, for improved path-planning algorithms.

Voxel maps can represent relative motions of two molecules. Such a representation displays the geometric suitability of a protein presenting a narrow, deep binding site for interacting with different ligands. When applied to explore molecular deformations, voxel maps can provide a global representation of the conformational space of protein loops and protein domains undergoing large-amplitude motions.

First results highlight the potential of the proposed approach. Voxel maps can be seen as a general tool that, combined with other computational and experimental methods, will help to investigate the importance of flexibility and motion in molecular interactions.

Fig. 2 .

 2 Fig.2. Voxel-map representation of channels connecting the active site and the surface in bacterial P450-BM3. Channel identifiers follow the nomenclature in related work[START_REF] Cojocaru | The ins and outs of cytochrome P450s[END_REF].

Fig. 3 .

 3 Fig. 3. a)-b) Models of the (R,S)-enantiomers of 4-methyloctanoic acid in the active site of CALB [33]. c)-d) Voxel maps representing locations of the center of mass of the (R,S)-enantiomers reachable from the catalytic position. Voxels size is 0.1 Å, and colors indicate the chronological order of generation.

Fig. 4 .

 4 Fig. 4. (Right) Structure of PTPase (PDB ID: 1YPT), with the WDP loop in black color. (Left) Detail of the WDP loop and voxel map displaying the explored locations of the Cα atom of Glu357. Voxel colors indicate the distance between Cα atoms of Glu357 and Val407. The most open and closed loop conformations are represented in red and blue respectively.

Fig. 5 .

 5 Fig. 5. Structure of N-Oct-3 POU domain.

Fig. 6 .

 6 Fig. 6. Two views of the voxel map representing geometrically feasible motions of POUh with respect to POUs. The voxels display the relative position of the centers of mass of both domains. Voxels size is 2 Å, and colors have been assigned depending on the energies of the associated conformations. The black voxel indicates the initial conformation.

Results reported in[START_REF] Cojocaru | The ins and outs of cytochrome P450s[END_REF] indicate that CAVER presents a significant sensibility to the staring point.

ACKNOWLEDGMENTS

The authors thank Maurice Franssen, David Guieysse and Magali Remaud for providing models and information about CALB. This work has been partially supported by the Région Midi-Pyrénées under projets ITAV ALMA and CTP AMOBIO, and by the French National Agency for Research (ANR) under project NanoBioMod.