Charles Poli

Matthieu Bellec

Ulrich Kuhl

Fabrice Mortessagne

Henning Schomerus

Selective enhancement of topologically induced interface states

The recent realization of topological phases in insulators and superconductors has raised the prospects to advance robust quantum technologies. The desire to demonstrate the underlying topological features with a high level of control has given incentive to explore optical platforms for analogous realizations. Here we show that the functionality of optical systems can be enhanced by combining topological protection with non-hermitian symmetries that do not have an electronic counterpart. This is achieved by combining parity-time symmetric losses with a unique feature of topologically induced interface states, namely, that they break a sublattice symmetry. This property isolates the state from the losses and enhances its visibility both in the frequency and in the time domain.

1

A key ingredient for topological functionality are robust confined states that form at interfaces between regions with topologically distinct band structures. For electromagnetic waves, this can be realized in two dimensions by breaking symmetries in analogy to the quantum Hall effect [START_REF] Wang | Observation of unidirectional backscattering-immune topological electromagnetic states[END_REF][START_REF] Fang | Realizing effective magnetic field for photons by controlling the phase of dynamic modulation[END_REF][START_REF] Rechtsman | Photonic Floquet topological insulators[END_REF] or the quantum spin Hall effect [START_REF] Hafezi | Robust optical delay lines with topological protection[END_REF][START_REF] Khanikaev | Photonic topological insulators[END_REF][START_REF] Hafezi | Imaging topological edge states in silicon photonics[END_REF], while in one dimension one can employ lattice modulations [START_REF] Malkova | Observation of optical Shockleylike surface states in photonic superlattices[END_REF][START_REF] Kraus | Topological states and adiabatic pumping in quasicrystals[END_REF][START_REF] Kitagawa | Observation of topologically protected bound states in photonic quantum walks[END_REF]. The intrinsic robustness of these topological states raises the question whether they can be controlled and modified independently of the other states in the system. Here, we draw on passive non-hermitian parity-time(PT)symmetry [START_REF] Guo | Observation of PT -symmetry breaking in complex optical potentials[END_REF][START_REF] Rüter | Observation of parity-time symmetry in optics[END_REF][START_REF] Feng | Experimental demonstration of a unidirectional reflectionless parity-time metamaterial at optical frequencies[END_REF][START_REF] Regensburger | Parity-time synthetic photonic lattices[END_REF] to demonstrate the selective control and enhancement of a topologically induced state in a one-dimensional microwave set-up [START_REF] Bellec | Tight-binding couplings in microwave artificial graphene[END_REF][START_REF] Laurent | Localized modes in a finite-size open disordered microwave cavity[END_REF][START_REF] Bellec | Topological transition of Dirac points in a microwave experiment[END_REF][START_REF] Franco-Villafañe | First experimental realization of the Dirac oscillator[END_REF]. The state can be isolated from losses affecting all other modes in the system, which enhances its visibility in the temporal evolution of a pulse, even in presence of structural disorder.

A minimal one-dimensional model with a topological band structure is a chain of sites with alternating couplings (i.e. a dimer chain). This model was originally introduced in an electronic context by Su, Schrieffer and Heeger (SSH) [START_REF] Su | Solitons in polyacetylene[END_REF] to describe fractionalized charges in polyacetylene, which appear in the presence of a dimerization defect. Photonic systems provide a versatile platform to realize analogies of this situation. In absence of PT-symmetry, the topological defect state has been observed in a quantum walk scenario [START_REF] Kitagawa | Observation of topologically protected bound states in photonic quantum walks[END_REF] while a dimer chain with a non-topological defect has been realised in a wave guide array [START_REF] Keil | The random mass dirac model and long-range correlations on an integrated optical platform[END_REF], and the PT-symmetric variant of the chain without a defect has been discussed in [START_REF] Ramezani | Pt-symmetric Talbot effect[END_REF]. Here, we implement the SSH chain in presence of the defect and PT-symmetry by means of a set of identical coupled dielectric resonators placed in a microwave cavity. Localized absorptive losses enable us to address non-hermitian effects that do not have an electronic analogue.

RESULTS

Realization of dimer chains by coupled disk resonators. Fig. 1a depicts a chain of 21 microwave resonators with a central dimerization defect. We establish a one-dimensional tight-binding regime [START_REF] Bellec | Tight-binding couplings in microwave artificial graphene[END_REF], where the electromagnetic field is mostly confined within the resonators. For an isolated disk, only a single mode is important in a broad spectral range around the bare frequency ν b = 6.65 GHz. This mode spreads out evanescently, so that the coupling strength can be controlled by adjusting the separation distance between the resonators [START_REF] Bellec | Tight-binding couplings in microwave artificial graphene[END_REF]. The resulting system can be described by the following tight-binding equa-tions:

(α)      (ν -ν n )ψ n = t 2 ψ n-1 + t 1 ψ n+1 , n = -2, -4, . . . (ν -ν n )ψ n = t 1 ψ n-1 + t 2 ψ n+1 , n = -1, -3, . . . , inter- face { (ν -ν 0 )ψ 0 = t 2 (ψ -1 + ψ 1 ), n = 0, (1) (β)
     (ν -ν n )ψ n = t 2 ψ n-1 + t 1 ψ n+1 , n = 1, 3, . . . , (ν -ν n )ψ n = t 1 ψ n-1 + t 2 ψ n+1 , n = 2, 4, . . . ,
where n enumerates the resonators, with n = -10, -9, -8, . . . 10 and n = 0 for the central site (see Fig. 1b). The mode amplitude in the nth resonator is given by ψ n , t 1 and t 2 denote the alternating nearest-neighbor coupling strengths, while ν n is the resonance frequency of the nth isolated resonator. Without absorption, the resonance frequencies are equal to the bare frequency, ν n = ν b . Absorption is introduced on selected sites by depositing elastomer patches on top of the dielectric cylinders (see Fig. 1a). The losses give rise to a complex

b -t 1 -t 2 < ν < ν b -|t 1 -t 2 | and ν b + |t 1 -t 2 | < ν < ν b + t 1 + t 2 .
The DOS is not affected when the values of the couplings t 1 and t 2 are interchanged. Nevertheless, a topological distinction between these two situations (called hereafter α and β configurations) can be captured by a winding number associated to the Bloch wave functions (see [START_REF] Ryu | Topological origin of zero-energy edge states in particle-hole symmetric systems[END_REF][START_REF] Schomerus | Topologically protected midgap states in complex photonic lattices[END_REF] and Supplementary Section S.1). An interface between both configurations takes the form of a dimerization defect where two consecutive couplings are identical (red zone in Fig. 1b and Fig. 2b). The topological distinctiveness of the two phases leads to the formation of an exponentially localized midgap state at ν = ν b .

The corresponding wavefunction can be read off from equation ( 1), and takes the form

ψ n = (-t 1 /t 2 ) -|n|/2
for even n and ψ n = 0 for odd n. The midgap state is therefore confined to the sublattice with even index, which we will call the A-sublattice, while the sites with odd index are called the B-sublattice. The complementary state on the B-sublattice increases exponentially and is incompatible with the boundary conditions.

Dimerization defect. Fig. 2b shows the DOS measured for a 21-resonator SSH chain with a central dimerization defect, still without absorption. Twenty-one modes are observed within the spectral range of interest. Of these, 20 modes are extended over the whole system.

These modes group into two sets of 10 and correspond to the upper and the lower band of the infinite dimer chain. The bands remain separated by a 45 MHz gap. The topologically induced mode clearly sits the middle of the gap, at frequency ν b . We find that the intensity of the midgap state belonging to the B-sublattice is zero within experimental resolution, and thus confirm that the wavefunction is confined to the A-sublattice. The corresponding wavefunction intensity profile pertaining to the A sites is depicted in Fig. 2d (red dots). As expected, the intensity decays according to an exponential profile given by the theoretical result (shown in gray).

Selective enhancement. We now set out to enhance the visibility of the topologically induced state. Our approach rests on the realization that the topological features of the system extend to a staggered configuration of losses, obtained by depositing absorptive material on all B sites [START_REF] Schomerus | Topologically protected midgap states in complex photonic lattices[END_REF]. Both in the α and in the β configuration, the tight-binding system then still possesses a passive (PT )-symmetry, given by a reflection (P) at a point in the middle of a dimer, which maps the passive A sites onto the lossy B sites but leaves the couplings unchanged. This mapping corresponds to a time-reversal operation (T ) up to a constant complex frequency shift iγ. As a consequence of this symmetry, all extended states are uniformly suppressed by the losses (their complex resonance frequencies all acquire the same imaginary part -γ/2), while the topologically induced interface state remains pinned at ν = ν b and thus is unaffected by the losses (see Supplementary Section S.1).

The spectral analysis presented in Fig. 2c shows that the extended modes shift downwards in frequency and become broadened, while the overall spectral weight in the resulting continuous bands is reduced. These features are consistent with the frequency dressing on the B sites, ν b → ν b -γ -iγ. In contrast, the peak in the density of states associated to the zero-mode remains fixed at the bare frequency ν b , and its height and width is almost unchanged. As shown in Fig. 2e this mode remains well confined to the A-sublattice, and still displays an exponential intensity profile as one moves away from the defect site. (Under the same conditions non-topological defect states hybridize, thereby degrading their properties, see Supplementary Section S.2.)

The spectral analysis implies that the interface state is insensitive to the losses distributed on the B sites. It then should become the predominant mode during the time evolution.

To illustrate this feature, Fig. 3 shows the time evolution for both non-absorptive (a) and absorptive (b) cases, corresponding to a pulse launched on the defect site. Without absorption, diffraction and interferences spoil the propagation of the interface state, which cannot be discerned after 250 ns. On the contrary, adding losses drastically enhances the visibility of the topologically induced mode, which then dominates the propagation without any degradation.

Robustness against structural disorder. To probe the topological protection property of the midgap state, we introduce structural disorder by randomly distributing the intersite separations. In the tight-binding description, this corresponds to a random modification of the coupling strengths. We preserve the dimer structure by defining the couplings as t r 1,2 = t1,2 + t(W/2)ξ, where t1 = (3t 1 + t 2 )/4 = 31.5 MHz, t2 = (t 1 + 3t 2 )/4 = 20.4 MHz and t = (t 1 -t 2 )/2 = 11.1 MHz. Here W is the disorder strength varying from 0 to 0.95 and ξ is a random number uniformly distributed in the interval [-1, 1] (see Fig. 4a). As shown in Fig. 4b, for different values of W the zero-mode intensities exhibit a very similar profile: an approximately exponential decay on the A sites and insignificant intensity on the B sites.

Even though we do not perform any averaging over disorder realizations, the experimental profile is in remarkable agreement with the simple exponential ( t1 / t2 ) -|n| . Note that as t1 / t2 = 1.54 < t 1 /t 2 = 2.51 the present zero-modes are more extended compared to the situation considered previously (Fig. 2d). When absorption is added on B sites, we observe in Fig. 4c that the robustness of the mode persists. Due to the resonance frequency shift of the lossy resonators, the couplings are now slightly changed to an effective coupling ratio ( t1 / t2 ) eff = 1.92, which is taken into account in the theoretical profile (gray shaded curve).

DISCUSSION

The loss-induced selective enhancement of the topological interface state observed in this work exploits the unique structure of the zero-mode wavefunction, which is confined to a sublattice. Such sublattice-symmetry breaking is common in topologically induced states; for example, it also occurs for the 0th Landau level in the quantum Hall effect of massless relativistic particles [START_REF] Jackiw | Fractional charge and zero modes for planar systems in a magnetic field[END_REF], as well as in inhomogeneously strained graphene [START_REF] Guinea | Energy gaps and a zero-field quantum Hall effect in graphene by strain engineering[END_REF] and photonic analogues of deformed honeycomb lattices [START_REF] Rechtsman | Strain-induced pseudomagnetic field and photonic Landau levels in dielectric structures[END_REF][START_REF] Schomerus | Parity anomaly and Landau-level lasing in strained photonic honeycomb lattices[END_REF]. The combination of topological constraints and passive (PT )-symmetry therefore provides a generic concept which may be exploited in different settings. This also extends to atom-optical systems, where the defectless version of the passive SSH model has recently been realized in an optical lattice [START_REF] Atala | Direct measurement of the Zak phase in topological Bloch bands[END_REF], while the nonpassive case has been discussed in theoretical work [START_REF] Barnett | Edge-state instabilities of bosons in a topological band[END_REF]. We note that it is also possible to selectively suppress the interface state, which can be achieved either by placing the losses on the other sublattice or by interchanging the couplings (topological phases) on both sides of the defect. The enhancement or suppression of the state can therefore be used to detect the relative size of two coupling strengths. It is also attractive to replace the losses by amplification, e.g. in a layered structure of materials with different thickness, which could be used to realize a laser with a much simplified mode competition, or by nonlinear effects as occurring e.g. in chains of coupled quantum dot or quantum well exciton polaritons [START_REF] Walker | Exciton polaritons in semiconductor waveguides[END_REF].

METHODS

Microwave realization of tight-binding systems. The experimental setup and the tight-binding description of the microwave system are detailed in [START_REF] Bellec | Tight-binding couplings in microwave artificial graphene[END_REF]. We describe the main ingredients. The sites of the lattice are occupied by dielectric microwave resonators with a cylindrical shape (5 mm height, 8 mm diameter, and a refractive index of 6). The resonance frequency of an isolated resonator ν b is around 6.65 GHz and corresponds to the on-site energy of atoms in a tight-binding model. The dielectric cylinders are coupled by the evanescent electromagnetic field, the corresponding coupling strength t between two resonators is well described by a tight-binding-like hopping term; t depends on the separation d between resonators. Via a reflection measurement, one has access, at each site, to the local density of states and to the wavefunction intensity associated to each eigenfrequency. The density of states (DOS) is obtained by averaging the local density of states over all resonator positions. To obtain the time evolution of the pulse we measure the transmission between a source located at the interface site and a receiver successively placed at each site position.

By performing a Fourier transform, one obtains the temporal evolution of a pulse initiating from the defect site and propagating into the SSH chain. In all these experiments, we face an intrinsic on-site disorder of ∼ 0.15% in the values of ν b .

  resonance frequency ν n = ν ′ b -iγ, and also shift the real part of the bare frequency to ν ′ b ≈ ν b -γ. In our experiments, γ ≃ 40 MHz, while the separations d 1 = 12 mm and d 2 = 15 mm correspond to couplings t 1 = 37.1 MHz and t 2 = 14.8 MHz, respectively. Defect-free chains. As a reference situation, we first consider a defectless SSH chain without absorption. As depicted in Fig. 2a, the density of states (DOS) is characterized by a band structure with two bands separated by a finite gap of size 2|t 1 -t 2 | = 45 MHz (gray zone). The extended states occupy bands in the range ν

FIG. 1 .FIG. 2 . 10 FIG. 3 .FIG. 4 .

 121034 FIG.1. Complex SSH chain. a. Picture of the experimental microwave realization of the complex SSH chain. The lattice is composed of 21 identical coupled dielectric cylindrical resonators (5 mm height, 8 mm diameter, and a refractive index of 6) sandwiched between two metallic plates (note that the top plate is not shown). To implement dimer chains, the resonators are separated byd 1 or d 2 with d 1 < d 2 , i.e. coupling t 1 > t 2 .A central dimerization defect is introduced. The defect creates an interface state at zero energy, a "zero-mode", whose visibility is enhanced by means of absorptive patches placed on one of the two sublattices. The resulting wavefunction intensity is superimposed onto the chain. b. Schematic of the complex SSH chain, with A and B-sublattices indicated in white and gray respectively. The strong (weak) coupling strength is represented by a thick (thin) line. In our system, the couplings can be controlled by varying the resonator spacings.The topologically induced zero-mode appears at the interface (red) between α configuration (with strong intradimer coupling) and β configuration (with weak intradimer coupling).

Acknowledgments. CP and HS acknowledge support of EPSRC via grant EP/J019585/1.

Author Contributions. CP and HS developed the concept and theory, MB and FM performed the experiments and analyzed the data. All the authors participated in the discussions and in writing the paper.

Competing financial interests. The authors declare no competing financial interests.


Julian Böhm

Matthieu Bellec

Fabrice Mortessagne

Ulrich Kuhl
email: ulrich.kuhl@unice.fr

Sonja Barkhofen

Stefan Gehler

Hans-Jürgen Stöckmann

Iain Foulger

Sven Gnutzman

Gregor Tanner

Experiments on 'quantum' search and directed transport in microwave artificial graphene

Keywords: numbers: 03.67.Ac, 72.80.Vp, 03.67.Lx, 42.50.Md

A series of quantum search algorithms has been proposed recently providing an algebraic speed-up compared to classical search algorithms from N to √ N where N is the number of items in the search space. In particular, devising searches on regular lattices have become popular extending Grover's original algorithm to spatial searching. Working in a tight-binding setup, it could be demonstrated theoretically, that a search is possible in the physically relevant dimensions 2 and 3 if the lattice spectrum possess Dirac points. We present here a proof of principle experiment implementing wave search algorithms and directed wave transport in a graphene lattice arrangement. The idea is based on bringing localized search states in resonance with an extended lattice state in an energy region of low spectral density -namely at or near the Dirac point. The experiment is implemented using classical waves in a microwave setup containing weakly coupled dielectric resonators placed in a honeycomb arrangement, i.e. artificial graphene. We furthermore investigate the scaling behavior experimentally using linear chains.

Introduction.-Currently, one of the most fruitful branches of quantum information is the field of quantum search algorithms. It started with Grover's work [1] describing a search algorithm for unstructured databases which has been implemented experimentally in NMR [2,3] and in optical experiments [4,5]. More recently, spatial quantum search algorithms have been proposed based on quantum walk mechanism [6,[START_REF]Quantum Walks and Search Algorithms[END_REF]. All these algorithms can achieve up to quadratic speed-up compared to the corresponding classical search. For quantum searches on generic d-dimensional lattices certain restrictions have been observed, however, depending on whether the underlying quantum walk is discrete [START_REF] Aharonov | [END_REF] or continuous [9]. While effective search algorithms for discrete walks on square lattices have been reported for d ≥ 2 [START_REF] Ambainis | Proceedings of the Sixteenth Annual ACM-SIAM Symposium on Discrete Algorithms, SODA '05[END_REF][START_REF] Abal | [END_REF], continuous-time quantum search algorithms on the same lattice show speed-up compared to the classical search only for d ≥ 4 [12].

Experimental implementations of discrete quantum walks need time stepping mechanisms such as laser pulses [13][14][15][16][17][18]. By switching to a continuous-time evolution based, for example, on tight-binding coupling between sites one can avoid time-discretization in an experiment. It has been shown in [19] that continuous-time quantum search in 2D is indeed possible when performed near the Dirac point in graphene or more generally for lattices with cone structure in the dispersion relation [20], i.e., a linear growth of the density of states (DOS). This effect adds a new dimension to the material properties of graphene [21,[START_REF] Wong | Carbon Nanotube and Graphene Device Physics[END_REF] with potential applications in sensing and detection as well as directed charge carrier transport. This may provide new ways of channeling intensity and information across lattices and between distinct sites, such as for single-molecule sensing as described in [START_REF] Schedin | [END_REF]24]. In this Letter, we present the first experimental realization of a 2D search in a tight-binding setup based on a microwave experiment using artificial graphene as discussed in [25][26][27].

In the following, we recapitulate briefly the theory of quantum searching on graphene and describe the experimental setup. We then demonstrate both searching as well as directed transport in graphene-like lattice structures. To explore the scaling behavior with the number of sites, we consider a linear chain demonstrating explicitly the √ N scaling. Theoretical background.-All quantum search algorithms, starting from Grover's search on an unstructured database, are based on the same principle: the system is setup by bringing an extended state of the (unperturbed) system in resonance with a localized state originating from a perturbation, thereby forming an avoided crossing in the spectrum of resonances. One then uses the 2-level dynamics at the avoided crossing in order to rotate the system from the extended state into the localized state thus 'finding' the position of the perturbation [28,29]. The subtleties in setting up such a search lie in i) choosing an unperturbed system with eigenstates extending uniformly across all sites; ii) finding a suitable perturbation which carries localized state and iii) working in an energy range with a low density of states making it possible to isolate the two-level crossing from the rest of the spectrum.

In the following, we will focus on continuous-time walks on regular, finite lattices in a tight-binding setup. The entries in the database are represented by the N lattice sites and associate orthonormal states {|i } N i=1 spanning an N dimensional Hilbert (search) space. The system is described by a tight-binding Hamiltonian H 0 with associate eigenstates {|e } N e=1 . The extended eigenstates most suitable for a search are those corresponding to reciprocal lattice-vectors at k = 0 or near the band edge. Introducing a perturbation of the form H = H 0 +W with W supporting a localized state at the 'marked site' being in resonance with a uniform extended state. The interaction at the avoided crossing is controlled by the overlap integral | i|e | 2 ≈ N -1 , which leads to an energy gap of order 1/ √ N . The search is started by preparing the system in an extended state |e . The state then evolves according to |ψ(t) = e -iHt/h |e describing a rotation taking place predominantly in the two-dimensional subspace spanned by |i , |e . Performing a measurement after the system has evolved into |i at a time t c = O( √ N ) completes the search.

For generic lattices, the number of states with energy less than E in d dimensions scales typically as N E d/2 . One thus expects the first state above the ground state at an energy

E 1 = O(N -2/d
). An energy separation between the states at the avoided crossing (being at a distance O(N -1/2 )) and the rest of the spectrum is thus only possible for d ≥ 4 in the large N limit. Running the search at a Dirac point overcomes this problem; here the energy scaling is reduced to E 1 = O(N -1/d ) and searching becomes possible for d ≥ 2 [19,20]. In contrast to [19], the perturbation W introduced in this letter is in terms of coupling an additional site to the lattice state |i with the perturber's on-site energy tuned to an extended state close to the Dirac point.

Experiments on graphene.-The experimental microwave setup is shown in the inset of Fig. 1. A metallic plate supports ceramic cylinders of height h = 5 mm and radius r = 4 mm which have a high index of refraction (n ≈ 6) thus acting as resonators. Microwaves are coupled into the system via antennas, which excite the first transverse electric resonance, called TE 1 , at ν 0 ≈ 6.65 GHz of the resonators. The system is closed from above by a metallic top plate at a distance of h p = 16 mm (not shown). We form a 'graphene flake' by positioning 216 resonators in a hexagonal lattice, see inset of Fig. 1, thus creating artificial graphene [25][26][27]. The top plates holds a loop antenna which can be positioned arbitrarily in the xy plane above each resonator. It couples via the magnetic field into the TE 1 mode. For details on the experimental setup and its relations to a tight-binding model refer to Ref. [27]. The reflection of the moving antenna at the center of the resonator is proportional to the intensity of the eigenmodes

|Ψ(r n )| 2 ,
where n labels the resonator, thus allowing for a determination of the local density of states (LDOS) and by integration the DOS as well. A transmission measurement from the fixed antenna 1 to the movable antenna 2 yields also the phases. Thus, by means of a Fourier transform, we determine the time dependent propagator

P (r n , t) = |F T (S 12 (r n , ν))| 2
, where S 12 is the measured transmission spectrum taken at every resonator position r n of the movable antenna [30,[START_REF] Poli | Selective enhancement of topologically induced interface states[END_REF]. The DOS for the graphene flake close to the Dirac point is shown in Fig. 1.

One can clearly identify two isolated states, which are extended lattice states. The lower one is denoted by ν l and the upper one by ν u . We will use the frequencies of these states near the Dirac point as working points for our search algorithm.

In the next step, we attach two perturbations at the side of the flake, see the inset of Fig. 1. The perturbation in the foreground is a single resonator with eigenfrequency ν rs close to ν u , the one in the background consists of a dimer, that is, two strongly coupled resonators with its lower frequency ν dim adjusted to ν l (see also dashed lines in Fig. 1). We have chosen this dimer configuration for having a parameter to fine-tune the frequency ν dim . The perturbations induce new resonance states interacting with the lattice states ν l and ν u . The time dependence of the search dynamics can now be demonstrated by Fourier transforming S 12 in a small frequency window around ν l or ν u . Fig. 2(a,c) shows the thus obtained initial state at t = 0 for the two different incident states. The state shown in Fig. 2(a) involves only frequencies close to ν l and ν dim and in Fig. 2(c) only frequencies close to ν u and ν s are included. After some time t = T dim (t = T sr ) the dimer (the single resonator) is illuminated and the perturbation is thus found [see Fig. 2(b,d)]. The search times T dim and T sr are slightly different due to differences in the coupling of the perturbation to the corresponding extended states.

By working near the Dirac point, the search can be extended to an arbitrary number of resonators N in principle; in praxis, the size of our model system is limited due to the overall absorption (quality facotr arround 1000). A detailed analysis of how the search time scales with N is thus not possible here. We will demonstrate this in a slightly different setup using a linear chain at the end of the Letter.

Our experiment demonstrates that the effect of spatial searching can be achieved in a graphene-like setup with tight-binding interaction. This is encouraging as the model here differs significantly from the theoretical study presented in [19] due to different boundary conditions, experimental uncertainties and absorption. Our results point towards a completely new set of applications in the actual carbon material -graphene -where the limitations due to absorption are less severe. In addition to forming the basis of a storage device with fast searching facility as demonstrated above, the results presented in Fig. 2 can also be interpreted in the framework of directed electron transport. One could use the device as a sensitive switch, where current will be directed either to the upper (dimer) or lower (single resonator) port depending on the Fermi energy.

Spatial quantum search can also be used for communication and 'quantum state transfer' [19,29,[START_REF]Quantum State Transfer and Network Engineering[END_REF]. The lattice is perturbed by two equivalent resonators, whose eigenfrequencies are tuned to a single eigenfrequency of the unperturbed lattice. The perturbers are attached at two different positions to the lattice and interact only via the lattice state. Launching an initial pulse in the lattice state will illuminate at the same time the two resonators but with reduced brightness. If we, however, prepare the initial state in one of the perturbing resonators and let the system evolve thereafter, the pulse will actually travel from the initial resonator via the lattice to the second resonator. The evolution is presented in Fig. 3. This time the upper resonators [Fig. 3(a)] is illuminated initially, we then go to a state living both in the lattice and on the two perturbing resonators [Fig. 3(b)]. Then the other resonator lights up [Fig. 3(c)] and (nearly) the total amplitude is transferred from one resonator to the other. This opens the way towards directed signal transfer and control in graphene; of course, graphene can not yet be manipulated on a single atom level as would be necessary for making use of the effects described in this letter. Our results may serve as incentive to trigger research efforts.

Experiment on linear chains.-In order to show the √ N scaling behavior of the search time in our experimental setup, we go to a simpler, quasi one-dimensional system, a linear chain. Searching is possible here only for small N , as the distance between neighboring resonances scales like N -2 and the number of resonances will flood the avoided crossing eventually for large N [19]. However, for small N < N cut , the effect can be nicely demonstrated as shown below. A photograph of the unperturbed chain with N = 11 resonators is depicted in Fig. 4(a); that is well below the cut-off which is N cut = 27 for the setup shown here. The reflection spectrum 1-|S 11 | 2 of this unperturbed chain is given in Fig. 4(b) as a black line (central spectrum), where we observe the 11 resonances of the regular chain. The central frequency for a chain with an odd number of resonators is always at the eigenfrequency of the single resonator and thus does not depend on N . The experimentally measured lattice mode which we use for the wave search realization is superimposed on the photograph. It is the central mode which for an odd number of states is given by |Ψ(r n )| = 1/ (N + 1)/2 for n odd and 0 otherwise. Additionally the reflection spectrum measured on a single state is presented (red, lower spectrum). The eigenfrequency matches the central resonance of the unperturbed chain. When attaching a perturber resonator to the chain, we find an additional resonance interact with the lattice state resulting in a resonance splitting (blue, top spectrum). The propagator is calculated by Fourier transforming the frequency range marked by the dotted vertical lines in Fig. 4(b) containing only central mode and the perturber state.

Fig. 5(a) depicts the temporal behavior of the intensity of the whole chain (black line) and of the perturber resonator (red line). The initial state (t = 0) is adjusted to the maximal amplitude on the lattice. We observe an oscillatory behavior corresponding to a beating between the two states. The search time T s is given by half the beating time T . The insets show the corresponding intensities of the wave function at the indicated times. We perform the same kind of measurements for different number of Conclusion.-We have demonstrated for the first time an experimental realization of 'quantum' wave search on an (artificial) graphene lattice. The search is facilitated by bringing a lattice state in resonance with a localized perturber state at an avoided crossing. The effects range from searching perturber states, to addressing particular sites using a frequency scan or initiating a switching behavior and by transferring signals between sites without knowing their positions. This is a proof of principle experiment limited mainly by losses and absorption. Reducing the resonance widths further, i.e., obtaining a better quality factor Q, can be realized using coupled supra conducting cavities (quality factors of about Q ≈ 10 8 are possible). One can also think of modeling graphene using elastodynamic systems based on silicon (Q ≈ 10 7 ).

  FIG. 1.(color online) The inset is a photograph of the graphene flake including the supporting metallic plate and the perturber resonators at the boundary (white arrows). The graph shows the DOS of the unperturbed flake, i.e. without single resonator and dimer attached. The resonance frequencies of the single resonator νsr and the dimer ν dim are marked by the dashed vertical lines.

FIG. 2 .

 2 FIG. 2. (color online) (a) and (c): The initial lattice state at t = 0 for two different initial frequency ranges is shown. (c) for the range including the dimer frequency ν dim and the lower lattice state ν l and (d) the range including νs and νu. (b) and (d) show the illuminated perturber state at the search time t = T dim and Tsr, respectively. See Supplemental Material at ... for a video showing the full dynamics (DimerSearchOn-Graphene.mp4 and MonomerSearchOnGraphene.mp4).

TimeFIG. 3 .

 3 FIG. 3. (color online) Intensity propagation in a graphene flake with two equal perturbing resonators attached. The initial state at t = 0 is localized at the upper resonator and is transferred via a lattice state to the state sitting on the opposite resonator. See Supplemental Material at ... for a video showing the full dynamics (StateTransferOnGraphene.mp4).

FIG. 4 .

 4 FIG. 4. (color online) (a) Experimental setup showing the linear chain with N =11 resonators without perturber. The top plate and antenna are not shown. The wave function corresponding to the central resonance which we use for the wave search is superimposed. (b) Reflection measurement of a single resonator (red, bottom), a regular chain with 11 resonators (black, center), and of the single resonator being attached to the regular chain (blue, top). The dotted lines mark the frequency range used for the Fourier transform to obtain the propagations. The baseline of the spectra are shifted correspondingly.

FIG. 5 .

 5 FIG. 5. (color online) (a) The time dependence of the intensity of the chain with 11 resonators (black line) and of the perturber (red line) is shown. The corresponding intensity distribution for specific times are shown as insets. See Supplemental Material at ... for a video showing the full dynamics (DiskSearchOnLinearChain.avi). (b) The beating time as a function of N ′ is shown. The dashed line corresponds to a linear increase whereas the solid line marks a square root dependance. The errorbars indicate the time difference between different recurrences.

Acknowledgements.-We would like to acknowledge inspiring discussion with Klaus Richter.


Didier Devaurs
email: devaurs@laas.fr

Amarda Shehu
email: amarda@gmu.edu

Thierry Siméon

Juan Cortés
email: jcortes@laas.fr

Sampling-based Methods for a Full Characterization of Energy Landscapes of Small Peptides

Keywords: energy landscape, peptide, sampling-based method

Obtaining accurate representations of energy landscapes of biomolecules such as proteins and peptides is central to structure-function studies. Peptides are particularly interesting, as they exploit structural flexibility to modulate their biological function. Despite their small size, peptide modeling remains challenging due to the complexity of the energy landscape of such highly-flexible dynamic systems. Currently, only samplingbased methods can efficiently explore the conformational space of a peptide. In this paper, we suggest to combine two such methods to obtain a full characterization of energy landscapes of small yet flexible peptides. First, we propose a simplified version of the classical Basin Hopping algorithm to quickly reveal the meta-stable structural states of a peptide and the corresponding low-energy basins in the landscape. Then, we present several variants of a robotics-inspired algorithm, the Transition-based Rapidly-exploring Random Tree, to quickly determine transition state and transition path ensembles, as well as transition probabilities between meta-stable states. We demonstrate this combined approach on the terminally-blocked alanine.

I. INTRODUCTION

Global thermodynamic and kinetic properties of molecules can be extracted from an analysis of their conformational energy landscapes [START_REF] Wales | Energy Landscapes: Applications to Clusters, Biomolecules and Glasses[END_REF]. In particular, obtaining an accurate representation of a molecule's energy landscape is a significant first step to conducting detailed structure-function studies for bio-molecules of central importance in the cell, such as proteins and peptides [START_REF] Shehu | Probabilistic search and optimization for protein energy landscapes[END_REF].

In this work, we focus on small peptides. Despite their modest size, they represent in many ways a more challenging setting than larger proteins. Peptides exhibit high structural flexibility, which enables them to recognize different molecular partners in the cell, and thus to modulate their biological function [START_REF] Paine | Prediction of the native conformation of a polypeptide by a statistical-mechanical procedure. III. Probable and average conformations of enkephalin[END_REF]. Contrary to proteins, that are often characterized by a unique native state and a funnel-shaped energy landscape, peptides are characterized by several meta-stable structural states; their energy landscape may contain a multitude of competitive low-energy basins.

The existence of multiple local minima in a molecule's energy landscape makes it particularly challenging to map this landscape and reconstruct all the functionally-important regions in it. Experimental methods, such as X-ray crystallography or nuclear magnetic resonance (NMR), cannot reveal such maps, as they can uncover only few structures at best [START_REF] Amodeo | Conformational sampling of bioactive conformers: A low-temperature NMR study of 15 N-Leu-Enkephalin[END_REF]. It is therefore the task of computational techniques to obtain detailed representations of energy landscapes.

Currently, only sample-based representations of the energy landscape can be afforded. Even for small peptides, the space of possible conformations is vast, and the effective degrees of freedom needed to represent a conformation are numerous. The high dimensionality of the space is accompanied by a complex (non-linear, non-convex) expression for the conformational energy, which is the result of competing local and non-local interatomic interactions. Probing this landscape is therefore very computationally-costly. Currently, only stochastic optimization techniques provide the right balance between accuracy and computational efficiency [START_REF] Wales | Energy Landscapes: Applications to Clusters, Biomolecules and Glasses[END_REF], [START_REF] Shehu | Probabilistic search and optimization for protein energy landscapes[END_REF].

Obtaining a representation of an energy landscape can be divided into two sub-problems: (1) determining meta-stable structural states (i.e. local energy minima); [START_REF] Shehu | Probabilistic search and optimization for protein energy landscapes[END_REF] computing transition paths between the identified states. Both can be addressed by achieving a dense yet efficient sampling of the conformational space.

In this paper, we propose to combine two sampling-based techniques to obtain a full characterization of energy landscapes of small yet highly-flexible peptides. First, we utilize our own variant of the classical Basin Hopping algorithm [START_REF] Wales | Global optimization by basin-hopping and the lowest energy structures of Lennard-Jones clusters containing up to 110 atoms[END_REF] to sample local minima in a peptide's energy landscape, and characterize low-energy basins in this landscape. Local minima are then organized via density-based clustering to reveal meta-stable structural states. Second, we use several variants of a robotics-inspired algorithm, the Transition-based Rapidly-exploring Random Tree (T-RRT) [START_REF] Jaillet | Randomized tree construction algorithm to explore energy landscapes[END_REF], to map out the connectivity between these states, thus completing the reconstruction of the peptide's energy landscape. In particular, we propose a new variant of T-RRT allowing to compute transition path ensembles and transition probabilities between meta-stable structural states. To illustrate this, and as a proofof-concept, we present results obtained on the terminallyblocked alanine, which is a frequent benchmark for studies in theoretical physical chemistry.

II. METHODS

The work presented in this paper is motivated by recent studies showing that robotics-inspired sampling-based algorithms can be a good basis for efficient conformational sampling and exploration in computational structural biology [START_REF] Gipson | Computational models of protein kinematics and dynamics: Beyond simulation[END_REF], [START_REF] Al-Bluwi | Motion planning algorithms for molecular simulations: A survey[END_REF]. The Transition-based RRT (T-RRT) algorithm is an example of such algorithms [START_REF] Jaillet | Sampling-based path planning on configuration-space costmaps[END_REF], [START_REF] Jaillet | Randomized tree construction algorithm to explore energy landscapes[END_REF]. It is based on the Rapidly-exploring Random Tree (RRT) [START_REF] Lavalle | Rapidly-exploring random trees: progress and prospects[END_REF], a popular path planning algorithm that can tackle complex problems in high-dimensional spaces. RRT has been successfully used in various disciplines, such as robotics, aerospace, computer animation, and computational structural biology. T-RRT is an extension of RRT involving a probabilistic transition test based on the Metropolis criterion. Like Metropolis Monte Carlo methods [START_REF] Frenkel | Understanding Molecular Simulations: From Algorithms to Applications[END_REF], it applies small moves to the molecular system; but, instead of generating a single path over the search space, it constructs a tree, providing a more efficient exploration. Moreover, the tree construction is intrinsically biased toward unexplored regions of the space, and favors expansions in lowenergy areas.

This work is partly based on an existing variant of T-RRT, called Multi-T-RRT [START_REF] Devaurs | A multi-tree approach to compute transition paths on energy landscapes[END_REF], and on a new extension to it, called Multi-T-RRT with Cycles. These algorithms are combined with an in-house variant of the Basin Hopping algorithm involving a simplified minimization process. Together they are used to obtain a complete representation of the energy landscape of a small yet flexible peptide. Details of these methods are presented in the next sections.

A. Basin Hopping

The Basin Hopping (BH) algorithm is a popular method for sampling local minima of an energy landscape. It was originally introduced to obtain the Lennard-Jones minima of small atom clusters [START_REF] Wales | Global optimization by basin-hopping and the lowest energy structures of Lennard-Jones clusters containing up to 110 atoms[END_REF]. Recently, BH has gained new attention to predict protein structure [START_REF] Olson | Evolutionary-inspired probabilistic search for enhancing sampling of local minima in the protein energy surface[END_REF], and to find intermediate structures of chemical reactions [START_REF] Kim | Efficient basin-hopping sampling of reaction intermediates through molecular fragmentation and graph theory[END_REF]. The method consists of repeatedly applying a structural perturbation followed by an energy minimization, which yields a trajectory of minima. The result is a (discrete) coarse-grained representation of the energy landscape that can be seen as a collection of interpenetrating staircases.

Our variant of BH differs from the classical one in that it does not involve local, gradient-based minimizations, but relies on simple Monte-Carlo-based (MC-based) minimizations. Our implementation of BH (that is presented in Algorithm 1) follows a random restart procedure performing several rounds, each one starting from a conformation randomly sampled in the search space. Every round builds a trajectory of minima by performing a succession of structural perturbations followed by MC-based minimizations. Every MC-based minimization starts from a conformation obtained by performing a large-amplitude perturbation of the minimum reached at the previous step, or from the random sample, in the first step. An MC-based minimization is an iterative succession of small-amplitude perturbations. The search goes on, accepting new conformations if the energy decreases or the Metropolis criterion is satisfied, until a given number of consecutive rejections is reached. Every MC-based minimization produces a low-energy conformation that we call a "minimum" in a minor abuse of language. Every round produces what we call a milestone: the minimum (along the trajectory) having the lowest energy.

All the milestones (or the minima) produced by BH have to be grouped to provide a comprehensible list of metastable structural states. This clustering can be done in several ways. In this work, we perform a density-based clustering requiring the user to define the distance threshold between two clusters [START_REF] Kim | Density-based clustering of small peptide conformations sampled from a molecular dynamics simulation[END_REF]. Distance computations are based on Root-Mean-Square Deviation (RMSD) and involve dihedral angles or atom coordinates. q ← doLargeAmplitudePerturbation(q) 7 q ← doMonteCarloBasedMinimization(q) 8 addMinimum(T , q) 9 addTrajectory(L, T ) 10 return L

B. Multi-T-RRT with Cycle-Addition

The Rapidly-exploring Random Tree (RRT) algorithm is a well-known path planning method that can deal with complex problems by performing an efficient exploration, even in high-dimensional search spaces [START_REF] Lavalle | Rapidly-exploring random trees: progress and prospects[END_REF]. Starting from an initial conformation q init , RRT performs its exploration by iteratively constructing a tree T that tends to rapidly expand over the conformational space C. The nodes and edges of T correspond to states (i.e. molecular conformations) and small-amplitude moves between states, respectively. At each iteration of the tree construction, a conformation q rand is randomly sampled in C. Then, an extension toward q rand is attempted, starting from its nearest neighbor q near , in T . This means performing a linear interpolation between q near and q rand , at a distance equal to the extension step-size δ, from q near . If the extension succeeds, a new conformation q new is added to T and an edge is built between q near and q new . The criteria on when to stop the exploration can be reaching a given target conformation q goal , a given number of nodes in the tree, a given number of iterations, or a given running time.

The Transition-based RRT (T-RRT) algorithm is a variant of RRT developed to explore a conformational space while taking the conformational energy into account [START_REF] Jaillet | Sampling-based path planning on configuration-space costmaps[END_REF], [START_REF] Jaillet | Randomized tree construction algorithm to explore energy landscapes[END_REF]. T-RRT extends RRT by integrating a stochastic transition test favoring the exploration of low-energy regions of the space. As the Metropolis criterion typically used by Monte Carlo methods, this transition test is used to accept or reject a candidate state. Here, the transitionTest presented in Algorithm 2 is used to evaluate the local move from q near to q new based on their respective energies E i and E j . More precisely, a downhill move (E j ≤ E i ) in the energy landscape is always accepted. An uphill move is accepted or rejected based on the probability e -(Ej -Ei) / (K•T ) (where K is the Boltzmann constant), which decreases exponentially with the energy variation E j -E i . In that case, the level of selectivity of the transition test is controlled by the temperature T , which is an adaptive parameter of the algorithm. Low temperatures limit the expansion to gentle slopes of the energy landscape, and high temperatures enable it to climb steep slopes. After each accepted uphill move, T is decreased to avoid over-exploring high-energy regions: it is divided by 2 (Ej -Ei) / energyRange(G) , where energyRange(G) is the energy difference between the highest-energy and the lowest-energy conformations in the Algorithm 2: transitionTest (G, E i , E j ) input : the current temperature T ; the temperature increase rate Trate; the Boltzmann constant K output: true if the transition is accepted, false otherwise

1 if Ej ≤ Ei then return True 2 if e -(E j -E i ) / (K•T ) > 0.5 then 3 T ← T / 2 (E j -E i ) / energyRange(G) ; return True 4 else 5 T ← T • 2 T rate ; return False Algorithm 3: Multi-T-RRT input : the conformational space C; the extension step-size δ; the energy function E : C → R; the initial conformations q k init , k = 1..n output: the tree T 1 for k = 1..n do 2 T k ← initTree(q k init ) 3 while not stoppingCriteria({T k | k = 1..n}) do 4 T ′ ← chooseNextTreeToExpand() 5 q rand ← sampleRandomConformation(C) 6 q ′ near ← findNearestNeighbor(T ′ , q rand ) 7 qnew ← extend(q ′ near , q rand , δ) 8 if qnew = null and 9 transitionTest(T ′ , E(q ′ near ), E(qnew)) then 10 addNewNode(T ′ , qnew) 11 addNewEdge(T ′ , q ′ near , qnew) 12 (T ′′ , q ′′ near ) ← findNearestTree(qnew) 13 if distance(qnew, q ′′ near ) ≤ δ then 14 T ← merge(T ′ , qnew, T ′′ , q ′′ near ) ; n ← n -1 15 return T
graph G built so far. After each rejected uphill move, T is increased to facilitate the exploration and avoid being trapped in a local energy minimum: it is multiplied by 2 Trate , where T rate ∈ (0, 1] is the temperature increase rate.

The Multi-T-RRT algorithm is a multiple-tree variant of T-RRT [START_REF] Devaurs | A multi-tree approach to compute transition paths on energy landscapes[END_REF]. Instead of building a single tree rooted at some initial conformation, the idea is to build n trees rooted at n given conformations q k init , k = 1..n. The pseudo-code of the Multi-T-RRT is presented in Algorithm 3. At each iteration, a tree T ′ is chosen for expansion in a round-robin fashion. Then, an extension is attempted toward a randomly sampled conformation q rand , starting from its nearest neighbor q ′ near , in T ′ . If the extension succeeds, the new conformation q new is added to T ′ , and connected to q ′ near . Then, we search for the conformation q ′′ near , which is the closest to q new within all trees other than T ′ . If the distance between q new and q ′′ near is less than or equal to the extension step-size δ, T ′ is linked to and merged with T ′′ , the tree to which q ′′ near belongs. In that case, the number of trees is decreased by 1. The space exploration continues until all trees are merged into a single one or another stopping criterion (number of nodes, number of expansions, running time) is met.

A drawback of the Multi-T-RRT is that it can produce "unnatural" transition paths (cf. Section III). To address this issue, we propose a new extension to T-RRT, based on a cycleaddition procedure. Starting from the tree produced by the Multi-T-RRT, the idea is to allow the space exploration to continue, and to add edges leading to the creation of cycles. This enables us to construct a graph from which several paths can be extracted between two given conformations. The pseudo-code of this T-RRT with Cycles is shown in Algorithm 4. It differs from T-RRT only in that, after every successful extension, an edge is added between q new and each conformation in the graph G that is not already connected to q new , and whose distance to q new is less than or equal to the extension stepsize δ. The stopping criteria involve the number of nodes and of expansions, as well as the running time. Finally, the Multi-T-RRT with Cycles consists of consecutively running the Multi-T-RRT and the T-RRT with Cycles.

As several paths may exist in the graph G between two given conformations, a quality criterion is required to compare paths. This is achieved by associating weights with the edges of G, based on the notion of mechanical work [START_REF] Jaillet | Sampling-based path planning on configuration-space costmaps[END_REF]. More precisely, the weight of the directed edge connecting q i and q j is equal to max{ 0 , E(q j ) -E(q i ) } , i.e. to the positive energy variation between q i and q j . This constitutes the amount of energy that has to be added to the molecule for the transition from q i to q j to happen. Note that using the mechanical work as a quality criterion requires to create two directed edges between q i and q j , instead of creating a single undirected edge. Finally, given two conformations in G, the best (directed) path linking them is defined as the one minimizing the mechanical work, and is obtained using Dijkstra's algorithm.

In this work, we use the mechanical work as a pathquality criterion because it has been shown that T-RRT tends to generate paths having a low mechanical work [START_REF] Jaillet | Sampling-based path planning on configuration-space costmaps[END_REF]. Furthermore, in many situations, the mechanical work can assess the quality of a path better than a simple criterion such as the integral of the cost along the path [START_REF] Jaillet | Sampling-based path planning on configuration-space costmaps[END_REF]. Other criteria could be considered, such as the minimum resistance (linked to the MaxFlux algorithm [START_REF] Huo | The MaxFlux algorithm for calculating variationally optimized reaction paths for conformational transitions in many body systems at finite temperature[END_REF], [START_REF]Direct computation of long time processes in peptides and proteins: reaction path study of the coil-to-helix transition in polyalanine[END_REF]) or the maximum flux [START_REF] Zhao | Maximum flux transition paths of conformational change[END_REF]. However, finding which criterion is the most relevant is out of the scope of this paper. Furthermore, it has been shown that, in the case of the terminally-blocked alanine, transition paths produced using different criteria are often similar [START_REF] Zhao | Maximum flux transition paths of conformational change[END_REF].

C. Force Field

To compute conformational energy values, we use an inhouse implementation of the AMBER parm96 force-field with an implicit representation of the solvent using the Generalized Born approximation. All results are obtained on an Intel Core i5 processor at 2.6 GHz with 8 GB of memory.

III. RESULTS AND DISCUSSION

A. Terminally-Blocked Alanine

As a proof of concept, we use our methods to explore the energy landscape of the terminally-blocked alanine. This molecule is an alanine residue acetylated in its N-terminus and methylamidated in its C-terminus: Ace-Ala-Nme (cf. Fig. 1). Despite its small size, it is a common test-model because of its complex energy landscape characterized by several local minima connected by multiple pathways involving various transition states [START_REF] Bolhuis | Reaction coordinates of biomolecular isomerization[END_REF], [START_REF] Chodera | Long-time protein folding dynamics from short-time molecular dynamics simulations[END_REF], [START_REF] Okumura | Temperature and pressure dependence of alanine dipeptide studied by multibaric-multithermal molecular dynamics simulations[END_REF], [START_REF] Gaigeot | Unravelling the conformational dynamics of the aqueous alanine dipeptide with first-principle molecular dynamics[END_REF], [START_REF] Velez-Vega | Kinetics and reaction coordinate for the isomerization of alanine dipeptide by a forward flux sampling protocol[END_REF].

The exploration of the conformational space involves an internal-coordinate representation of the terminally-blocked alanine, assuming constant bond lengths and bond angles. The conformational parameters are the dihedral angles {φ, ψ, χ} of the Ala residue, χ of the Ace capping, χ of the Nme capping, ω of the Ace-Ala peptide bond, and ω of the Ala-Nme peptide bond (cf. Fig. 1). As the peptide bond torsions are known to undergo only small variations, the ω angles are allowed to vary only up to 10 • from the planar trans conformation. The φ and ψ dihedral angles of the terminally-blocked alanine are critical, for their flexibility allows internal hydrogen bonds to form. Consequently, the energy landscape of the terminally-blocked alanine is often projected on these two coordinates to facilitate analysis [START_REF] Chodera | Long-time protein folding dynamics from short-time molecular dynamics simulations[END_REF], [START_REF] Okumura | Temperature and pressure dependence of alanine dipeptide studied by multibaric-multithermal molecular dynamics simulations[END_REF], [START_REF] Velez-Vega | Kinetics and reaction coordinate for the isomerization of alanine dipeptide by a forward flux sampling protocol[END_REF]. In a similar way, we compute a Ramachandran map by projecting the energy landscape of the terminally-blocked alanine on the 2-dimensional space defined by the (φ, ψ) coordinates (cf. Fig. 2). This 2-D map is generated by varying both dihedral angles with a 5 • step and finding the lowest-energy conformation corresponding to each (φ, ψ) pair. This local minimum is determined using a random restart procedure performing multiple, independent minimizations. Each minimization starts from a randomly sampled conformation, and generates new conformations by perturbing any dihedral angle (except φ and ψ) of the current minimum. A new conformation is accepted if its energy is lower than the current minimum, otherwise it is rejected. The minimization stops when a maximum number of consecutive rejections is reached. The six local minima visible on the (φ, ψ) Ramachandran map (cf. Fig. 2) correspond to six known meta-stable states of the terminally-blocked alanine, namely the C 5 , P II , α R , α P , C ax 7 , α L states [START_REF] Chodera | Long-time protein folding dynamics from short-time molecular dynamics simulations[END_REF], [START_REF] Okumura | Temperature and pressure dependence of alanine dipeptide studied by multibaric-multithermal molecular dynamics simulations[END_REF]. To determine the energy values and the (φ, ψ) coordinates of these minima, the exhaustive exploration of the (φ, ψ) projection of the energy landscape is refined with a 1 • step in the relevant areas. Results are reported in Table I.

B. Local Energy Minima of the Terminally-blocked Alanine

We now show that, using Basin Hopping, it is possible to quickly find the energy minima of the terminally-blocked alanine, as identified by the (very costly) exhaustive exploration of the (φ, ψ) space presented previously. The parameters of BH are set as follows: we perform 100 rounds, each round executing 10 MC-based minimizations whose maximal number of consecutive rejections is set to 10. Therefore, every run of BH produces 1000 minima and 100 milestones, in about 2 s (on average over 100 runs). The large-amplitude perturbations initiating MC-based minimizations affect only the φ or ψ dihedral angle. The small-amplitude perturbations performed during an MC-based minimization affect only one, randomly chosen angle among all dihedral angles. The clustering procedure used to group the milestones is based on an RMSD similarity measure involving the coordinates of the carbon atoms of the backbone of the terminally-blocked alanine, and on a similarity threshold of 0.25 Å. This method consistently produced six clusters, after each one of the 100 runs, and the cluster representatives fitted well the six target minima (see Fig. 3). To evaluate the precision of the method, we compute the average (and standard deviation) of the coordinates of these representatives, and compare them to the (φ, ψ) coordinates of the target minima. Results (not reported here due to space limitations) show that precision is very high. For each coordinate of any minimum, the difference between the average of the values returned by BH and the target value is never greater than 0.5 • (which is the error inherent to the resolution of the exhaustive search), and the standard deviation of the values returned by BH is never greater than 0.3 • .

We have also used BH to obtain a description of the catchment basins on the energy landscape of the terminallyblocked alanine. This enables us to characterize the lowenergy basins to which the meta-stable states belong, as well as the transition state ensembles between them. To do that, the parameters of BH are set as follows: we perform 100 rounds, each round executing 10 MC-based minimizations whose maximal number of consecutive rejections is set to 100. Therefore, a run of BH produces 1000 minima, which takes about 25 s. Clustering these minima (and not the milestones, as previously done) produces again the same six clusters. This allows grouping the associated starting points of the MC-based minimizations into six clusters, representing six catchment basins (see Fig. 4). A visual analysis of the result shows that these clusters provide a good characterization of the lowenergy basins (where the clusters do not overlap) and of the transition state ensembles between them (where the clusters overlap). These regions could be automatically characterized using segmentation techniques, but this is out of the scope of this paper.

C. Transition Paths of the Terminally-blocked Alanine

We now show that T-RRT can quickly produce many transition paths between meta-stable states of the terminallyblocked alanine. The conformations used as input for Multi-T-RRT are the energy minima produced by BH. The conformational distance is based on the L 2 norm in the (φ, ψ) space. The extension step-size δ is set to 0.1 rad, meaning that the maximal angular variation between two conformations is about 5.7 • . The Boltzmann constant being 1.987 • 10 -3 kcal/mol/K, by setting the initial temperature to 70 K, we impose the probability of accepting an energy increment of 0.1 kcal/mol to be around 50% at the beginning of the exploration. Finally, the temperature increase rate T rate is set to 0.1.

Starting with six trees, Multi-T-RRT quickly returns a single tree (in about 0.1 s on average over 100 runs) from which paths connecting the minima are extracted. These paths are then projected on the (φ, ψ) space and plotted on the Ramachandran map of the terminally-blocked alanine. We observe that some of these paths may look "unnatural": in the instance illustrated by Fig. 5 (top), the transition path from P II to α L goes through almost all the other minima. Using T-RRT with Cycles allows solving this issue, thanks to the alternative paths created in the graph. For instance, Fig. 5 (bottom) shows a graph obtained after letting T-RRT with Cycles run for 100 s. This method produces more direct transition paths between all pairs of minima. In fact, the longer the running time of T-RRT with Cycles, the more direct paths are found thanks to the increasing coverage of the space by the graph and to the increasing number of cycles in the graph. While remaining acceptable for a user, a running time of 100 s enables us to obtain a reasonably dense coverage of the space.

To get an idea of the diversity of the transition paths of the terminally-blocked alanine, we run both variants of Multi-T-RRT 100 times and aggregate all the extracted paths on the same Ramachandran map, as illustrated by Fig. 6. We can see that most low-energy and medium-energy regions of the landscape are covered by the transition paths. However, in the case of Multi-T-RRT, the paths going through the medium-energy area corresponding to φ ∈ [-100 • , -50 • ] and ψ ∈ [0 • , 50 • ] are under-represented in comparison to the case of Multi-T-RRT with Cycles. Finally, we observe that no transition path goes through the energetic barrier corresponding to φ ∈ [100 Multi-T-RRT with Cycles can be used to compute transition probabilities between all pairs of energy minima of the terminally-blocked alanine. This is done by running it 1000 times and counting how many runs yield a graph from which a direct transition path can be extracted between a given pair of minima. More precisely, we consider that a run produces a direct transition path between two minima if the best path in the graph between them does not go through another minimum. By varying the maximal running time of the algorithm, it is possible to determine which transitions are most or least likely. The results of this experiment are reported in Fig. 7. It appears that the most likely transitions involve the pairs (C 5 , P II ) and (α R , α P ), as already shown in [START_REF] Chodera | Long-time protein folding dynamics from short-time molecular dynamics simulations[END_REF]. Direct transition paths are easily found between the minima in these pairs, even when runtime is low. If we had no other information about the energy landscape, we could infer that these minima are not separated by high-energy barriers or other low-energy basins. On the other hand, it requires a longer time for the transition probabilities associated with (C ax 7 , α L ) to reach 1, meaning that high-energy barriers probably separate them. Furthermore, since transitions from α P to α R are very likely, it takes some time before the algorithm is able to produce transition paths from α P to C 5 or P II that do not go through α R . That is why the corresponding transition probabilities increase slowly with runtime. This is similar for the transitions going from C 5 to α P . Interpreting the other curves is difficult: it is impossible to know how many high-energy barriers or low-energy basins separate the corresponding minima without additional information.

D. Transition States of the Terminally-blocked Alanine

From the transition paths computed by Multi-T-RRT, we can extract some information about the transition state ensembles of the terminally-blocked alanine. First, we define a transition state as the conformation having the highest energy along a given path. Our experiment consists of running Multi-T-RRT 1000 times (which takes about 120 s on average) and aggregating on the Ramachandran map the transition states found along all direct transition paths between any pair of minima. To obtain ensembles, we group the transition states using the same clustering procedure as when grouping the minima produced by BH. Among all the clusters, we reject those containing less than 100 states (i.e. 10% of the number of runs) as being not significant. We also reject clusters containing more than 1000 states because they cannot correspond to a single transition state ensemble. Finally, we choose as cluster representatives the conformations having the lowest energy within each cluster. After repeating the whole procedure 100 times, we observed that it consistently produced four clusters, i.e. four transition state ensembles (as illustrated by Fig. 8).

Within each transition state ensemble, the conformation with the lowest energy corresponds to a saddle point on the energy landscape. Our experiment shows that the terminallyblocked alanine features exactly four energy saddle-points (see Fig. 8). To assess the precision of our procedure at determining the (φ, ψ) coordinates and the energy values of these saddle points, we compute averages and standard deviations across the 100 trials we performed. Results are reported in Table II.

From their standard deviations, we can infer that it is difficult to accurately estimate the coordinates of the saddle points. However, since the standard deviation of the energy values is very low (less than 0.05 kcal/mol), we can speculate that this is mainly due to the flatness of the saddles.

IV. CONCLUSION

In this paper, we have presented a methodology to explore and obtain a complete representation of the energy landscape of small flexible peptides. The methodology combines a variant of Basin Hopping (BH) with several extensions of a robotics-inspired algorithm: the Transition-based Rapidlyexploring Random Tree (T-RRT). The paper makes several algorithmic contributions. First, by implementing a simplified version of the BH algorithm, where local, gradient-based energy minimization steps are replaced by simple Monte-Carlobased minimization steps, we achieve a relevant exploration of the energy landscape yielding numerous samples around energy minima. This leads to a quick determination of metastable structural states and low-energy basins containing them. In addition, the multiple-tree version of T-RRT is very fast at generating transition paths between meta-stable states. This provides us with a good description of the transition path and transition state ensembles. Finally, by using an anytime variant of the Multi-T-RRT algorithm, we are able to estimate transition probabilities between meta-stable states.

The aim of this paper was to demonstrate that combining sampling-based algorithms such as BH and T-RRT allows quickly obtaining a full representation of the energy landscape of a small yet highly-flexible peptide. Studying the terminallyblocked alanine served as a proof-of-concept in this paper. Our next targets are bigger peptides, starting with the met-Enkephalin (Tyr-Gly-Gly-Phe-Met). Despite the great interest it has sparked due to its dynamical structure and its potential role in pain inhibition, alcoholism, and cancer treatment, this peptide is currently poorly understood [START_REF] Sauriyal | Extending pharmacological spectrum of opioids beyond analgesia: Multifunctional aspects in different pathophysiological states[END_REF]. The preliminary tests we have performed on the met-Enkephalin show that the methods we propose can remain efficient with bigger peptides. Indeed, when running the Multi-T-RRT using as input height conformations corresponding to local energy minima of the met-Enkephalin (produced by BH), we obtain a graph connecting them in 25 s (on average over 100 runs). If we disregard the solvent in our force field, a solution can be obtained in about 2 s (on average over 100 runs). As part of our ongoing work, we also have to ensure that the clustering procedures involved in our methods remain effective when dealing with bigger peptides.

Directions of future work include exploiting the graph produced by a single run of the Multi-T-RRT with Cycles to describe transition path and transition state ensembles. This will allow us to make better use of computational resources, as opposed to aggregating paths extracted from several runs of the Multi-T-RRT. In addition, Markov-based transition-step analysis can be conducted on the graph produced by one or more runs of T-RRT. This analysis can allow estimating stabilities of each computed state, and provide a rigorous basis for the designation of a state as stable or semi-stable.

Algorithm 1 : 3 T

 13 Basin Hopping input : the conformational space C the number of rounds nbRounds the number of Monte Carlo minimizations nbM C output: the list of trajectories of minima L 1 L ← φ 2 for r = 1..nbRounds do ← φ ; q ← sampleRandomConformation(C) 4 for m = 1..nbM C do 5 if m > 1 then

6

 6

Algorithm 4 :

 4 T-RRT with Cycles input : the conformational space C; the extension step-size δ; the energy function E : C → R; the tree built by the Multi-T-RRT T output: the graph G 1 G ← initGraph(T ) 2 while not stoppingCriteria(G) do 3 q rand ← sampleRandomConformation(C) 4 qnear ← findNearestNeighbor(G, q rand ) 5 qnew ← extend(qnear, q rand , δ) 6 if qnew = null and 7 transitionTest(G, E(qnear), E(qnew)) then 8 addNewNode(G, qnew) 9 addNewEdge(G, qnear, qnew) 10 for qcan ∈ G \ {qnew} do 11 if distance(qnew, qcan) ≤ δ and 12 noEdgeBetween(qnew, qcan) then 13 addNewEdge(G, qnew, qcan) 14 return G

Fig. 1 .

 1 Fig. 1. Internal-coordinate representation of the terminally-blocked alanine.

Fig. 2 .

 2 Fig.2. Projection of the energy landscape of the terminally-blocked alanine on the (φ, ψ) internal coordinates. For reference, the lowest energy value is set to 0 kcal/mol. Names used for the local minima are taken from[START_REF] Chodera | Long-time protein folding dynamics from short-time molecular dynamics simulations[END_REF],[START_REF] Okumura | Temperature and pressure dependence of alanine dipeptide studied by multibaric-multithermal molecular dynamics simulations[END_REF].

Fig. 3 .

 3 Fig. 3. Results of the exploration of the energy landscape of the terminallyblocked alanine produced by Basin Hopping. The minima generated by the Monte-Carlo-based minimizations are represented by grey disks, and the selected milestones by yellow disks. Among them, the red disks are the cluster representatives. For reference, the labeled black circles show the six (target) minima identified by the exhaustive search.

Fig. 4 .

 4 Fig. 4. Characterization of the catchment basins on the energy landscape of the terminally-blocked alanine, obtained with Basin Hopping. The colored disks represent the starting points of the MC-based minimizations. The color of a disk identifies the cluster to which the minimum produced by the corresponding MC-based minimization belongs.

Fig. 5 .

 5 Fig. 5. Graphs generated by two variants of T-RRT exploring the conformational space of the terminally-blocked alanine. Top: graph produced by Multi-T-RRT (which stops when all trees are merged) in about 0.1 s. Bottom: graph produced by T-RRT with Cycles after running for 100 s. The edges of the graph are represented by thin grey lines. The transition path from P II to α L (extracted from the graph) is highlighted in black.

Fig. 6 .

 6 Fig. 6. Transition paths generated by 100 runs of Multi-T-RRT (top) and 100 runs of Multi-T-RRT with Cycles (bottom) exploring the conformational space of the terminally-blocked alanine. The maximal running time was set to 100 s for the latter.

Fig. 7 .

 7 Fig.7. Transition probabilities (expressed as percentages) between all pairs of minima of the terminally-blocked alanine, in relation to the running time of Multi-T-RRT with Cycles (in seconds, on a logarithmic scale). Probabilities are computed on the basis of 1000 runs. In each plot, the title indicates the initial conformation of the transition paths extracted from the graph, and the five curves correspond to the possible goal conformations (see color legend).

Fig. 8 .

 8 Fig. 8. Transition state ensembles (represented by black dots) of the terminally-blocked alanine, obtained by aggregating the results of 1000 runs of Multi-T-RRT, and clustering the extracted transition states (i.e. the highestenergy conformations along all direct transition paths). Transition states rejected by the clustering procedure appear as grey dots. The red disks correspond to the cluster representatives, i.e. the lowest-energy conformations in all clusters; they are energy saddle points.

TABLE I .

 I ENERGY AND (φ, ψ) COORDINATES OF THE LOCAL MINIMA OF THE TERMINALLY-BLOCKED ALANINE. NAMES ARE TAKEN FROM[START_REF] Chodera | Long-time protein folding dynamics from short-time molecular dynamics simulations[END_REF],[START_REF] Okumura | Temperature and pressure dependence of alanine dipeptide studied by multibaric-multithermal molecular dynamics simulations[END_REF]. FOR REFERENCE, THE LOWEST ENERGY IS SET TO 0 KCAL/MOL.

		C5	PII	αR	αP	C ax 7	αL
	φ ( • )	-141	-62	-57	-140	44	43
	ψ ( • )	163	150	-53	-72	-107	61
	E (kcal/mol)	0	0.44	1.11	1.52	3.27	3.48

TABLE II .

 II ENERGY AND (φ, ψ) COORDINATES OF THE SADDLE POINTS OF THE TERMINALLY-BLOCKED ALANINE, AVERAGED OVER 100 TRIALS OF A PROCEDURE BASED ON THE MULTI-T-RRT. STANDARD DEVIATION IS ALSO GIVEN FOR THE COORDINATES. STANDARD DEVIATION OF THE ENERGY WAS LESS THAN 0.05 KCAL/MOL.

	φ ( • )	67 ± 3	1 ± 2	0 ± 3	69 ± 3
	ψ ( • )	128 ± 5	92 ± 2	-88 ± 3	-13 ± 4
	E (kcal/mol)	6.2	6.8	6.9	7

ACKNOWLEDGMENT

This work has been partially supported by the French National Research Agency (ANR) under project ProtiCAD (project number ANR-12-MONU-0015).