Adrien Morel
email: adrien.morel@cea.fr

Gaël Pillonnet

Pierre Gasnier

Elie Lefeuvre

Adrien Badel
email: adrien.badel@univ-smb.fr

Frequency tuning of piezoelectric energy harvesters thanks to a short-circuit synchronous electric charge extraction

Keywords: Energy harvesting, Nonlinear electrical interface, Piezoelectricity, Damping control, Frequency tuning, Synchronous Electric Charge Extraction

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Wireless sensor nodes designed to facilitate monitoring, risk prevention, and more generally improve comfort in our day-to-day life, are increasingly present in our everyday environment. Most of these devices require a fairly small amount of energy to sense, process, then send an information. This makes energy harvesting an attractive alternative to electrochemical batteries, as it could lead to devices combining compactness, long lifetime, and robustness to harsh environments [START_REF] Penella | A review of commercial energy harvesters for autonomous sensors[END_REF]. In closed and confined environments, where the thermal gradients and solar radiations are relatively weak, the major energy source can be found under the form of mechanical vibrations [START_REF] Vullers | Micropower energy harvesting[END_REF]. This may be the case in biomedical applications, for instance, to power implantable medical devices which are placed inside the human body. Here, compactness as well as autonomy are undeniably vital [START_REF] Olivo | Energy harvesting and remote powering for implantable biosensors[END_REF].

Among electromechanical transducers, piezoelectric materials are of interest because of their good trade-off between power density and technological maturity [START_REF] Roundy | Improving power output for vibration-based energy scavengers[END_REF]. In order to amplify the mechanical strain applied on such materials, they are usually deposited on a linear mechanical oscillator having a relatively high quality factor. The harvested power from such electromechanical system is optimal when the damping induced by the energy extraction is equal to the mechanical intrinsic damping in the structure, and when the ambient vibration frequency matches the resonant frequency of the mechanical oscillator [START_REF] Mitcheson | Tuning the Resonant Frequency and Damping of an Electromagnetic Energy Harvester Using Power Electronics[END_REF]. However, the harvested power drastically decreases when a slight mismatch appears between those two frequencies. This mismatch may be caused by many reasons: inaccurate design of the harvester leading to a difference between the expected and real resonant frequencies, aging of the mechanical oscillator, change in the ambient vibration spectrum, among others [START_REF] De Godoy | Effect of parametric uncertainties on the performances of a piezoelectric energy harvesting device[END_REF].

In order to face this challenge, researchers have proposed mechanical solutions based on nonlinear oscillators to enlarge the harvesting frequency bandwidth. The nonlinearities are usually induced thanks to nonlinear springs [START_REF] Boisseau | Nonlinear H-shaped springs to improve efficiency of vibration energy harvesters[END_REF], buckled beam configurations [START_REF] Ando | A bistable buckled beam based approach for vibrational energy harvesting[END_REF], or magnetic forces [START_REF] Lin | The magnetic coupling of a piezoelectric cantilever for enhanced energy harvesting[END_REF], leading, in the most common cases, to softening, hardening, and bistable behaviors [START_REF] Harne | A review of the recent research on vibration harvesting via bistable systems[END_REF]. Nonetheless, these solutions are not flawless, as the hysteretic frequency response of such oscillators may lead to low displacement orbits, hindering the benefits of this approach [START_REF] Daqaq | On the role of nonlinearities in vibratory energy harvesting: A critical review and discussion[END_REF]. Furthermore, the shapes of those harvesters' power frequency responses are highly dependent on the vibrations' amplitude.

Recently, some researchers proposed to adapt the dynamics of the piezoelectric energy harvester using the electrical interface connected to the piezoelectric material. Indeed, thanks to the indirect piezoelectric effect, the voltage across the piezoelectric material induces a force on the mechanical oscillator. Depending on its phase and amplitude, this force has an influence on both the harvester resonant frequency and damping [START_REF] Badel | Nonlinear Conditioning Circuits for Piezoelectric Energy Harvesters[END_REF]. The phase and amplitude of the piezoelectric force depending on the equivalent resistive and reactive components of the electrical load connected to the piezoelectric material, this phenomenon can also be understood as an electromechanical impedance matching, as explained extensively in [START_REF] Liang | Impedance Modeling and Analysis for Piezoelectric Energy Harvesting Systems[END_REF].

Based on these principles, W.-J. Wu et al. first proposed in 2006 to adapt a linear capacitive load connected to the harvester, which allows to electrically tune the stiffness of the harvester, and impacts its electromechanical resonant frequency [START_REF] Wu | Tunable resonant frequency power harvesting devices[END_REF]. Further studies of the capacitive compensation have been proposed later, in 2009 [START_REF] Liao | Optimal parameters and power characteristics of piezoelectric energy harvesters with an RC circuit[END_REF] and 2013 [START_REF] Ahmed-Seddik | Self-powered resonant frequency tuning for piezoelectric vibration energy harvesters[END_REF]. A. Badel et al. proposed in 2014 a nonlinear interface called the Frequency Tuning Synchronous Electric Charge Extraction (FTSECE), which consisted in adapting the amount of electric charges collected every semi-period of vibration, and the phase of the harvesting events [START_REF] Badel | Wideband Piezoelectric Energy Harvester Tuned Through its Electronic Interface Circuit[END_REF]. Y. Cai et al. proposed in parallel a similar approach, analyzing this strategy using equivalent impedances methods [START_REF] Cai | System design of a time-controlled broadband piezoelectric energy harvesting interface circuit[END_REF]. In 2015, P. H. Hsieh et al. proposed to add a tunable switching delay to the Parallel Synchronized Switch Harvesting on Inductance (P-SSHI) [START_REF] Hsieh | Improving the scavenged power of nonlinear piezoelectric energy harvesting interface at off-resonance by introducing switching delay[END_REF]. Thereafter, B. Zhao et al. proposed to combine the idea of adding a tunable switching delay with the synchronized multiple bias-flip (SMBF) concept [START_REF] Zhao | Phase-variable Control of Parallel Synchronized Triple Bias-Flips Interface Circuit towards Broadband Piezoelectric Energy Harvesting[END_REF], which consists in optimizing the charges inversion processes thanks to multiple voltage-changing steps [START_REF] Liang | Synchronized bias-flip interface circuits for piezoelectric energy harvesting enhancement: A general model and prospects[END_REF]. In 2017, our group developed a novel approach to electrically tune the resonant frequency of piezoelectric energy scavengers and enlarge their harvesting bandwidth: the N-SECE, which consists in harvesting the energy at a different frequency than the vibration's frequency [START_REF] Morel | A unified N-SECE strategy for highly coupled piezoelectric energy scavengers[END_REF]. Recently, we proposed another concept, the SC-SECE, which relies on a tunable short-circuit phase [START_REF] Morel | Short Circuit Synchronous Electric Charge Extraction (SC-SECE) Strategy for Wideband Vibration Energy Harvesting[END_REF].

In this paper, we propose a thorough and extended analysis of the SC-SECE, and we compare it with traditional electrical interfaces (Full bridge rectifier, SECE), and other strategies to enlarge the bandwidth (i.e. using nonlinear mechanical oscillators and other tunable electrical interfaces). SC-SECE allows to adapt both the resonant frequency and the damping of a linear Piezoelectric Energy Harvester (PEH), thanks to a short-circuit phase whose start angular time 𝜙 " and angular duration 𝛥𝜙 are both adjustable. First, we propose a model of the voltage and displacement waveforms associated with the SC-SECE, and derive an analytical expression of the harvested power, depending on the harvesting frequency, the harvester characteristics, and the two electrical parameters 𝜙 " and 𝛥𝜙. Thereafter, we numerically analyze our model in order to determine the power frequency responses of PEH associated with the electrical interface implementing the SC-SECE strategy. In a last part, we experimentally verify the validity of our analytical predictions, using a highly coupled piezoelectric harvester connected to an electrical interface implementing the SC-SECE strategy. The experimental results are then commented, compared to the theoretical ones, and the proposed SC-SECE strategy is finally compared to other wideband solutions.

Electromechanical system model

Piezoelectric harvester modelling

A linear piezoelectric energy harvester is usually made of a piezoelectric material deposited upon a mechanical resonator. Around its first resonant frequency, such a system can be modeled by a single mechanical degree of freedom system, as depicted in Fig. 1 [START_REF] Wang | Frequency Analysis of Vibration Energy Harvesting Systems[END_REF].

Figure 1 -Linear piezoelectric harvester made of a cantilever beam and its equivalent single degree of freedom modelling

The mechanical resonator can be, for instance, a cantilever beam of stiffness 𝐾 ", whose free end has been fixed to an inertial mass 𝑀, and subject to frictional forces modelled by a linear damping coefficient 𝐷. Due to the indirect piezoelectric effect, the piezoelectric material exerts a force on the mechanical resonator, which is, in first approximation, proportional to the voltage across the piezoelectric element, 𝑣 0 . Under an external acceleration 𝛾, the second order differential equation governing the displacement 𝑥 of the inertial mass is then given by (1).

-Mγ = Mx + Dx + K ;< x + αv ? [START_REF] Penella | A review of commercial energy harvesters for autonomous sensors[END_REF] Where 𝛼 is the force factor, and is dependent on both the mechanical resonator's geometry and the beam and piezoelectric material's properties. The piezoelectric material, under the strain induced by the resonator oscillation, generates electrical charges. These charges are either stored in the material dielectric capacitance 𝐶 0 , or harvested by the electrical interface. The differential equation governing the electrical dynamics of the harvester is given by [START_REF] Vullers | Micropower energy harvesting[END_REF]. αx = C ? v ? + i ? [START_REF] Vullers | Micropower energy harvesting[END_REF] Where 𝑖 0 accounts for the current drawn by the electrical interface. The electromechanical model introduced by (1) and ( 2) is accurate as long as the mechanical linearity hypothesis is true, the dielectric losses in the piezoelectric material are small and the influences of higher mechanical modes remain negligible. Under a harmonic excitation, and as we assume that the mechanical oscillator has a relatively important quality factor, we assume that the inertial mass displacement is only a sinusoidal signal of magnitude 𝑋 F , as expressed by (3).

x = X H 𝑐𝑜𝑠 ωt = X H 𝑐𝑜𝑠 (θ) [START_REF] Olivo | Energy harvesting and remote powering for implantable biosensors[END_REF] Where 𝜔 is the angular frequency of the harmonic excitation. In order to facilitate analytical developments, several dimensionless variables have been introduced in previous works [START_REF] Badel | Nonlinear Conditioning Circuits for Piezoelectric Energy Harvesters[END_REF], and are summarized in [START_REF] Roundy | Improving power output for vibration-based energy scavengers[END_REF].

ω S = K ;< /M Ω H = ω/ω S k H W = α W /(K ;< C ? ) Q H = K ;< M/D (4)
Where 𝜔 S is the short-circuit angular resonant frequency of the mechanical resonator, and 𝛺 F is the normalized ambient vibration frequency, 𝑘 F W is the modified coupling factor of the harvester, and corresponds to the ratio of the maximum electrostatic energy stored in the piezoelectric material divided by the maximum elastic energy stored in the mechanical resonator under quasi static excitation. 𝑄 F accounts for the quality factor of the mechanical resonator. The maximum harvestable power expression with a linear PEH excited under a harmonic acceleration has been analytically determined in previous works [START_REF] Badel | Nonlinear Conditioning Circuits for Piezoelectric Energy Harvesters[END_REF][START_REF] Arroyo | Comparison of electromagnetic and piezoelectric vibration energy harvesters: Model and experiments[END_REF] and is given by ( 5). This limit can only be reached when the harvester is working at resonance, and when the damping induced by the electrical interface is equal to the mechanical damping. [START_REF] Lefeuvre | Piezoelectric Energy Harvesting Device Optimization by Synchronous Electric Charge Extraction[END_REF] as an alternative to Standard Energy Harvesting (SEH) interfaces, as it allows to harvest more energy for lowly coupled / highly damped piezoelectric harvesters, when the product 𝑘 F W 𝑄 F < 𝜋/4. Furthermore, contrarily to SEH and Synchronized Switch Harvesting on Inductance (SSHI) interfaces, its efficiency does not depend on the voltage across the storage element, and hence does not require a Maximum Power Point Tracking (MPPT) loop working constantly to optimize the harvested power [START_REF] Guyomar | Recent progress in piezoelectric conversion and energy harvesting using nonlinear electronic interfaces and issues in small scale implementation[END_REF][START_REF] Wang | Dimensionless Optimization of Piezoelectric Vibration Energy Harvesters with Different Interface Circuits[END_REF]. In the standard SECE strategy, the interface circuit is not connected to the piezoelectric element during most of the vibration's period. Most of the time, the piezoelectric harvester works in open-circuit configuration, and all the charges generated, thanks to the piezoelectric effect, are stored in the dielectric capacitance 𝐶 0 . When the voltage reaches an extremum, an inductance is briefly connected to the piezoelectric material, forming a LC resonator. All the electrostatic energy stored in the piezoelectric material is then quickly transferred in the inductance. When the voltage becomes 0V, meaning that all the energy has been transferred from the piezoelectric material to the electrical circuit, the inductance is disconnected from the piezoelectric material, the energy stored in the inductance is transferred in a storage capacitor and the harvester restarts to operate in open-circuit configuration. The energy harvested during a single harvesting event can be calculated using [START_REF] De Godoy | Effect of parametric uncertainties on the performances of a piezoelectric energy harvesting device[END_REF].

P H]^= Mγ W 8D = γ W Q H M 8 ω S (
W defe = 1 2 C ? v H W (6)
Where 𝑣 F stands for the piezoelectric voltage extremum value. The piezoelectric voltage and inertial masse displacement waveforms, as well as a typical interface circuit implementing the SECE strategy are shown in Fig. 2 and Fig. 3. If the SECE strategy is efficient when associated with lowly coupled / highly damped PEHs, it has however been proven that for highly coupled / lowly damped piezoelectric harvesters (𝑘 F W 𝑄 F > 𝜋/4), the SECE strategy draws too much power from the mechanical resonator [START_REF] Morel | Regenerative synchronous electrical charge extraction for highly coupled piezoelectric generators[END_REF]. This overdamps the mechanical resonator, leading to a reduced inertial mass displacement amplitude, and low harvested energy performances. Furthermore, the energy harvested with the SECE interface remains fairly limited when the piezoelectric harvester is excited outside the resonance.

SC-SECE strategy

The proposed SC-SECE strategy's aim is to overcome the standard SECE weaknesses thanks to the introduction of a tunable short-circuit phase. This short-circuit phase allows to limit the amount of energy converted from the mechanical resonator to the piezoelectric material, hence avoiding any overdamping effect when 𝑘 F W 𝑄 F > 𝜋/4 [START_REF] Morel | A unified N-SECE strategy for highly coupled piezoelectric energy scavengers[END_REF]. First, let's define 𝜙 " ∈ [0, 𝜋] as the angle at which the short-circuit phase starts. In order not to waste any energy, 𝜙 " is associated as well with the energy harvesting events, which does not necessarily happen on the voltage extrema anymore (corresponding to 𝜙 " = 0). We define as well 𝛥𝜙 ∈ [0, 𝜋] as the angular time spent in short-circuit condition. The piezoelectric voltage and inertial mass displacement waveforms are shown in Fig. 2, and a systemlevel implementation of the SC-SECE is shown in Fig. 3. Compared to the standard SECE system's implementation [START_REF] Lefeuvre | Piezoelectric Energy Harvesting Device Optimization by Synchronous Electric Charge Extraction[END_REF], the SC-SECE only requires an additional switch and a few control blocks. These additional requirements do not cost much space as they can all easily be integrated. The choice of the electrical parameters (𝜙 m , 𝛥𝜙) can either be done manually or could be automated thanks to an ultralow-power algorithm such as a gradient algorithm. From (2), the analytical expression of the piezoelectric voltage 𝑣 0 during a vibration's period [𝜙 " , 𝜙 " + 2𝜋) can be expressed as a piecewise function of 𝑥, as shown in [START_REF] Boisseau | Nonlinear H-shaped springs to improve efficiency of vibration energy harvesters[END_REF].

v ? = 0,
Since 𝑣 0 is a periodic function, we can find its associated Fourier series, expressed by [START_REF] Ando | A bistable buckled beam based approach for vibrational energy harvesting[END_REF].

v ? = a S + a z 𝑐𝑜𝑠 (nθ) | z}~+ b z 𝑠𝑖𝑛 (nθ) | z}~ (8)
Where 𝑎 S , 𝑎 ' , 𝑏 ' are the piezoelectric voltage series coefficients. Due to the filtering effect of the linear mechanical resonator, we assume that only the piezoelectric voltage fundamental component has an impact on the harvester dynamics. This assumption greatly simplifies the analytical calculations of the piezoelectric voltage. The fundamental 𝑣 0 ~ of the piezoelectric voltage 𝑣 0 can be expressed as follows:

v ? ~= a ~𝑐𝑜𝑠 θ + b ~𝑠𝑖𝑛 (θ) (9)
Applying Fourier series theory on [START_REF] Boisseau | Nonlinear H-shaped springs to improve efficiency of vibration energy harvesters[END_REF], 𝑎 ~ and 𝑏 ~ can be found by solving [START_REF] Harne | A review of the recent research on vibration harvesting via bistable systems[END_REF]. [START_REF] Harne | A review of the recent research on vibration harvesting via bistable systems[END_REF] Injecting ( 3) in [START_REF] Harne | A review of the recent research on vibration harvesting via bistable systems[END_REF] leads to the following Fourier coefficient expressions:

a ~= α X H π C ? π -Δϕ + 𝑠𝑖𝑛 2ϕ ; + 2Δϕ 2 + 𝑠𝑖𝑛 2ϕ ; 2 + 2 𝑐𝑜𝑠 ϕ ; + Δϕ 𝑠𝑖𝑛 (ϕ ; ) b ~= - α X H π C ? 𝑐𝑜𝑠 ϕ ; + 𝑐𝑜𝑠 ϕ ; + Δϕ W (11)
Thus, combining ( 9) and ( 11) leads to the expression of the fundamental 𝑣 0 ~ of the piezoelectric voltage in the frequency domain:

v ? ~= α x C ? ε … + jε ‡ ε … = 1 - Δϕ π + 𝑠𝑖𝑛 2ϕ ; + 2Δϕ 2π + 𝑠𝑖𝑛 2ϕ ; 2π + 2 𝑐𝑜𝑠 ϕ ; + Δϕ 𝑠𝑖𝑛 ϕ ; π ε ‡ = 𝑐𝑜𝑠 ϕ ; + 𝑐𝑜𝑠 ϕ ; + Δϕ W π (12)
Where 𝑥, 𝜀 ‰ , 𝜀 Š are the displacement written in the frequency domain, and the normalized damping and stiffness induced by the electrical interface, respectively. Combining (1) and ( 12) in the frequency domain, and isolating 𝑥 leads to the inertial mass displacement expression. Normalizing the obtained expression using the dimensionless variables introduced by (4), we get [START_REF] Liang | Impedance Modeling and Analysis for Piezoelectric Energy Harvesting Systems[END_REF].

x = ω S W γ 1 -Ω H W + k H W ε … + j Ω H /Q H + k H W ε ‡ (13)
Thus, the magnitude of the inertial mass displacement is expressed by [START_REF] Wu | Tunable resonant frequency power harvesting devices[END_REF].

X H = ω S W γ 1 -Ω H W + k H W ε … W + Ω H /Q H + k H W ε ‡ W (14)
The harvested energy, traditionally given by expression [START_REF] De Godoy | Effect of parametric uncertainties on the performances of a piezoelectric energy harvesting device[END_REF] for the standard SECE, becomes [START_REF] Liao | Optimal parameters and power characteristics of piezoelectric energy harvesters with an RC circuit[END_REF] in the case of the SC-SECE strategy.

W df‹defe = 1 2 C ? v ? ϕ ; W = α W X H W 2 C ? 𝑐𝑜𝑠 ϕ ; + 𝑐𝑜𝑠 ϕ ; + Δϕ W (15)
The power harvested thanks to the SC-SECE can thus be derived:

P df‹defe = ω π W df‹defe = ω α W X H W 2 π C ? 𝑐𝑜𝑠 ϕ ; + 𝑐𝑜𝑠 ϕ ; + Δϕ W (16)
Combining ( 14) and ( 16), we can get the following power expression:

P df‹defe = k H W Mγ W ω 𝑐𝑜𝑠 ϕ ; + 𝑐𝑜𝑠 ϕ ; + Δϕ W 2 π K ;< 1 -Ω H W + k H W ε … W + Ω H /Q H + k H W ε ‡ W (17)
This power expression can be normalized with respect to the maximum harvestable power, whose expression is given by [START_REF] Mitcheson | Tuning the Resonant Frequency and Damping of an Electromagnetic Energy Harvester Using Power Electronics[END_REF]. The normalized power 𝑃 '•ŽF ∈ [0,1] is given by [START_REF] Cai | System design of a time-controlled broadband piezoelectric energy harvesting interface circuit[END_REF]. 𝑃 '•ŽF gives an insight about how efficient the mechanical-to-electrical conversion is.

P z••H = P df‹defe P H]^= 4 k H W Ω H 𝑐𝑜𝑠 ϕ ; + 𝑐𝑜𝑠 ϕ ; + Δϕ W Q H π 1 -Ω H W + k H W ε … W + Ω H /Q H + k H W ε ‡ W (18)

Numerical analysis of the tuning parameters influences

Based on the previously derived analytical expression of the harvested power (18), we compute the power frequency responses of the proposed SC-SECE strategy. For every vibration's frequency, the parameters (𝜙 " , 𝛥𝜙) are optimized in order to maximize the harvested power. For low coupling harvesters, the power frequency responses of the SECE and SC-SECE strategies are almost superimposed. Thus, in this case, there is little to no interest to include a shortcircuit after the harvesting event. As predicted, both the SECE and SC-SECE strategies are more efficient than the SEH. As the coupling increases, it can be observed that the maximum power of the standard SECE strategy decreases, which is due to the overdamping phenomena extensively explained in [START_REF] Morel | A unified N-SECE strategy for highly coupled piezoelectric energy scavengers[END_REF][START_REF] Morel | Regenerative synchronous electrical charge extraction for highly coupled piezoelectric generators[END_REF]. In another hand, the SC-SECE maximum power remains equal to the maximum harvestable power 𝑃 F'' . Furthermore, the bandwidth of the SC-SECE becomes more and more important as the coupling is increased. The product of the maximum normalized power 𝑃 " = 𝑚𝑎𝑥(𝑃 '•ŽF ) times 𝑄 F is a good indicator of how much power can be harvested from a PEH associated with a particular electrical strategy [START_REF] Badel | Nonlinear Conditioning Circuits for Piezoelectric Energy Harvesters[END_REF]. Indeed, 𝑃 " gives an insight of the mechanical-to-electrical conversion efficiency, while the net maximum power that can be harvested is directly proportional to the quality factor of the mechanical resonator 𝑄 F , as shown by [START_REF] Mitcheson | Tuning the Resonant Frequency and Damping of an Electromagnetic Energy Harvester Using Power Electronics[END_REF]. This product has been computed in Fig. 7 as a function of 𝑘 F W , and for various mechanical quality factor 𝑄 F . Here as well, we can observe that for lowly coupled / highly damped piezoelectric harvesters (when the product 𝑘 F W 𝑄 F is small), both the SECE and SC-SECE allow to harvest more power than the SEH strategy. When the coupling gets higher, the maximum harvested power with the SECE interface starts decreasing, while both the SC-SECE and the SEH strategy allow to harvest the maximum harvestable power (𝑃 '•ŽF = 1). This confirms the trends that appears in Fig. 4, Fig. 5 and Fig. 6. For the highly coupled PEH introduced in next section (𝑄 F = 20, 𝑘 F W = 0.4), the SC-SECE allows to harvest up to 2.6 times more energy than with the SECE strategy. Fig. 8 shows the product of the normalized maximum harvested power 𝑃 " with the quality factor of the mechanical resonator 𝑄 F and with the normalized harvesting bandwidth 𝛥𝛺/𝛺 -, where 𝛥𝛺 corresponds to the normalized frequency range where the harvested power is greater than half of the maximum harvested power, and 𝛺 -the central frequency of this frequency range. This product is computed as a function of the squared coupling coefficient 𝑘 F W and with various mechanical quality factors 𝑄 F . We can observe that 𝑃 " 𝑄 F 𝛥𝛺/𝛺 -is generally increasing with 𝑘 F W . We can however notice that for high 𝑘 F W values, this product stops increasing for the SECE interface, because of the overdamping effect. For high values of 𝑄 F and 𝑘 F W , this product may even start decreasing in the case of SEH. In another hand, in the case of the SC-SECE, this product gets always greater with a higher 𝑘 F W , and is always more important than with any other strategy. This confirms the potential of the SC-SECE to both optimize the harvested power 𝑃 " from lowly and highly coupled PEH, while drastically increasing the power harvesting bandwidth 𝛥𝛺/𝛺 -.

Experimental validation

Experimental setup

Experimentations have been conducted in order to validate the proposed SC-SECE strategy. The experimental setup shown in Fig. 9 consists of an electromagnetic shaker, which converts the electrical signal delivered by the DSpace real-time control and acquisition board and amplified by a power amplifier into an acceleration. The input voltage delivered by the DSpace board is regulated in order to maintain a constant shaker's acceleration amplitude, thanks to an internal control loop. A highly coupled piezoelectric harvester made with 10x5x0.5mm 3 plates of PZN-PT single crystals by Microfine Company (Singapore) poled along the [011] axis and exhibiting a very high transverse coupling coefficient 𝑘 -~ of 0.92 along the [100] axis has been designed, and is depicted in Figure 10. This harvester has been fixed on the electromagnetic shaker, and its electrodes are connected to an electrical circuit made of discrete components which allows to obtain the voltage waveforms of the proposed SC-SECE strategy. Here, the circuit only allows to extract energy, without storing it, the stored energy optimization as well as the interface topology choice being beyond the scope of this paper. The electrical signals controlling the two tunable parameters (𝜙 " , 𝛥𝜙) are sent from the DSpace board, and the voltage across the piezoelectric harvesters, as well as the positions of the harvester and its inertial mass are recorded through the DSpace board and stored on a computer for each parameter couple combination, and for each vibration frequency. The experimentations are run for 30 values of 𝜙 " between 0 and 𝜋, 30 values of 𝛥𝜙 between 0 and 𝜋, and 100 vibration's frequencies 𝑓 ™š› between 90Hz and 140Hz. The piezoelectric energy harvester's parameters have been summarized in Table 1.

Experimental results

At first, the shaker has been excited under a relatively low acceleration amplitude 𝛾 = 0.03𝑔. Some experimental voltages and displacement waveforms are shown in Fig. 11 for various (𝜙 " , 𝛥𝜙). 𝑣 ,•'žŽ•Ÿ represents the voltage signal controlling the short-circuit transistor (S 2 in Fig. 3) of the SC-SECE interface. When 𝑣 ,•'žŽ•Ÿ is set high, the piezoelectric element is short-circuited, and when 𝑣 ,•'žŽ•Ÿ is set low, the piezoelectric element works in open-circuit configuration. We can observe that the displacement amplitude is sinusoidal, which validates the assumption formulated in (3). Even though the higher mechanical modes were considered negligible in the theoretical part, we can see their impact on the piezoelectric voltage waveforms, on which ripples can be observed. However, these influences remain fairly limited, and does not impact much the harvesting performances of the SC-SECE, as it will be shown further. -, WS For various vibration frequency, the extracted power from the piezoelectric harvester has been determined experimentally with various (𝜙 " , 𝛥𝜙), as shown in Fig. 12. The maximum harvested power with the SC-SECE as well as with the optimal tunable parameters (𝜙 " •0ž , 𝛥𝜙 •0ž ) under a small acceleration of 0.03G is shown in Fig. 11. We can see that the experimental results are in good agreement with the harvested power predictions given by ( 18), even if some noticeable differences can be observed due to higher mechanical modes excitations. This is due to the non-linear treatment of the piezoelectric voltage which stimulates higher harmonics of the mechanical system. Under small vibrations of 0.03G, the maximum harvested power reaches 6𝜇𝑊. The proposed SC-SECE strategy allows to harvest more than 3𝜇𝑊 (which corresponds to 50% of the maximum harvested power) between 90Hz and 140Hz, which makes the harvesting bandwidth 𝛥𝛺/𝛺 , as large as 50% of the short-circuit resonant frequency of the harvester. We can observe a bifurcation of the optimal parameters when 𝛺 F ∈ [0.98, 1.15]. Indeed, on this frequency range, two distinct couples of parameters (𝜙 " , 𝛥𝜙) are leading to two voltage waveforms that share the same (and optimal) first harmonic characteristics (phase and amplitude) even though their spectral contents are different. Experimentally, one of the solutions (straight line) is most of the time a bit more optimal than the other one (dashed line), due to the higher voltage harmonics which have a small impact on the harvested power but have not been considered in the analytical model, due to the first harmonic approximation.

We ran the same experiment with a higher acceleration amplitude, 𝛾 = 0.08𝐺, as shown in Fig. 13. Figure 13 -Experimental (crosses) and theoretical (straight lines) power frequency response of the proposed SC-SECE under small vibrations, with 𝛾 = 0.08𝐺, and optimal couples (𝜙 " , 𝛥𝜙).

Overall, the experimental results remain in good agreement with the theoretical model. The maximum harvested power under 0.08G is 35𝜇𝑊, which is a bit lower than what could be predicted by [START_REF] Mitcheson | Tuning the Resonant Frequency and Damping of an Electromagnetic Energy Harvester Using Power Electronics[END_REF]. Indeed, under higher acceleration, the quality factor of the mechanical resonator 𝑄 F decreases, and is here around 15. The harvesting bandwidth 𝛥𝛺/𝛺 , remains around 50%, from 90Hz to 140Hz. On this frequency band, the harvested power is always greater than 15𝜇𝑊.

The same experiment has been finally run with 𝛾 = 0.2𝐺. The results are shown below, in Fig. 14. Under this stronger vibration, the experimental results remain in relatively good agreement with the theoretical predictions, even though we start to observe a shift of the mechanical resonant frequency of the harvester, which seems to have decreased from 100Hz to 96Hz. This is particularly apparent when comparing the experimental and theoretical (𝜙 " , 𝛥𝜙). This shift of the mechanical resonant frequency is due to some mechanical nonlinearities: as the vibration amplitude gets larger, the global stiffness of the beam decreases (softening behavior).

The maximum harvested power under 0.2𝐺 is 185𝜇𝑊, which is slightly lower than what could be predicted by [START_REF] Mitcheson | Tuning the Resonant Frequency and Damping of an Electromagnetic Energy Harvester Using Power Electronics[END_REF]. Once again, the harvesting bandwidth 𝛥𝛺/𝛺 , is close to 50%, while the harvested power is greater than 92𝜇𝑊 between 90Hz and 140Hz. Those experimental results, summarized in Fig. 15, prove that the harvesting bandwidth as well as the power frequency responses of the harvester do not change much with the vibration's amplitude, which confirms the robustness and potential of the SC-SECE approach.

Comparison with other state-of-the-art wideband solutions

Previous approaches for enlarging the harvesting bandwidth of energy harvesters usually consist in replacing the linear mechanical resonator with a nonlinear one. These nonlinearities can for instance be induced thanks to magnets repulsion and attraction forces, to nonlinear springs, or buckled beam structures. It has been proven in the last decade that such nonlinear harvester can drastically enhance the harvesting bandwidth compared to linear PEHs [START_REF] Harne | A review of the recent research on vibration harvesting via bistable systems[END_REF]. Recently, new electrically-based approaches have been investigated such as the FTSECE [START_REF] Badel | Wideband Piezoelectric Energy Harvester Tuned Through its Electronic Interface Circuit[END_REF] and the N-SECE [START_REF] Morel | A unified N-SECE strategy for highly coupled piezoelectric energy scavengers[END_REF]. In this section, we compare the experimental results obtained with those different electrical and mechanical approaches with the proposed SC-SECE approach.

𝑆𝐹𝑜𝑀 §¨= 𝐸 ª¨ 𝛥𝛺 𝛺 , 𝑄 F 𝜌 ¬-- (19)
Where 𝐸 ª¨= 𝑃 '™ /𝑃 H]^ corresponds to the energy harvester effectiveness criterion, 𝑃 '™ being the average harvested power over the considered bandwidth, and 𝑃 H]^ being the maximum available power as defined by [START_REF] Mitcheson | Tuning the Resonant Frequency and Damping of an Electromagnetic Energy Harvester Using Power Electronics[END_REF]. 𝛥𝛺/𝛺 , corresponds to the normalized harvesting bandwidth, as defined previously. 𝑄 F is the quality factor of the harvester and is a good indicator of what the maximum available power 𝑃 H]^ is, as formulated in [START_REF] Mitcheson | Tuning the Resonant Frequency and Damping of an Electromagnetic Energy Harvester Using Power Electronics[END_REF]. 𝜌 ¬--is the effective mass density, defined as 𝜌 ¬--= 𝑀/𝑉, 𝑀 being the inertial mass of the harvester, and 𝑉 the total volume which can be approximated as the sum of the PEH's inertial mass volume 𝑉 F'"" and the motion space 𝑉 F•žš•' . In the case of our highly coupled PEH associated with the SC-SECE interface, 𝐸 ª¨= 0.8, 𝛥𝛺/ 𝛺 , =0.48, 𝑄 F = 20, 𝜌 ¬--= 5700𝑘𝑔. 𝑚 ‹-, leading to a 𝑆𝐹𝑜𝑀 §¨ of about 43776𝑘𝑔. 𝑚 ‹-. Figure 16 compares the 𝑆𝐹𝑜𝑀 §¨ of state-of-the-art wideband solutions with the proposed SC-SECE.

Figure 16 -𝑆𝐹𝑜𝑀 §¨ of mechanical [START_REF] Stanton | Nonlinear dynamics for broadband energy harvesting: investigation of a bistable piezoelectric inertial generator[END_REF][START_REF] Arrieta | A piezoelectric bistable plate for nonlinear broadband energy harvesting[END_REF][START_REF] Mann | Investigations of a nonlinear energy harvester with a bistable potential well[END_REF][START_REF] Liu | Novel piezoelectric bistable oscillator architecture for wideband vibration energy harvesting[END_REF][START_REF] Liu | A wideband integrated piezoelectric bistable generator: experimental performance evaluation and potential for real environmental vibrations[END_REF][START_REF] Masana | Relative performance of a vibratory energy harvester in mono-and bistable potentials[END_REF][START_REF] Ashraf | A wideband, frequency up-converting bounded vibration energy harvester for a low-frequency environment[END_REF][START_REF] Sebald | Simulation of a Duffing oscillator for broadband piezoelectric energy harvesting[END_REF] and electrical [17, 22, 39, 40, this paper] solutions

Comparison of the SC-SECE with other wideband electrical solutions

Among the electrical solutions, the most important 𝑆𝐹𝑜𝑀 §¨ is exhibited by the FTSECE design proposed in [START_REF] Badel | Wideband Piezoelectric Energy Harvester Tuned Through its Electronic Interface Circuit[END_REF]. However, these results are only theoretical, and have not been experimentally validated. The N-SECE proposed in [START_REF] Morel | A unified N-SECE strategy for highly coupled piezoelectric energy scavengers[END_REF] exhibits a lower 𝑆𝐹𝑜𝑀 §¨, but only requires a single parameter tuning which is much easier to implement. The FTSECE proposed in [START_REF] Cai | A Piezoelectric Energy-Harvesting Interface Circuit with Fully Autonomous Conjugate Impedance Matching, 156% Extended Bandwidth, and 0.38µW Power Consumption[END_REF], while being selfpowered and autonomous, still exhibits a relatively important 𝑆𝐹𝑜𝑀 §¨. This design is a proof that an ultra-low-power implementation of such electrical wideband solutions is feasible. Among the electrical solutions which have been experimentally validated, the proposed SC-SECE has the greatest 𝑆𝐹𝑜𝑀 §¨, which confirms the potential of the SC-SECE compared to other electrical approaches.

General comparison of electrical solutions with wideband mechanical solutions

As it can be seen in Fig. 16, electrically-based wideband solutions 𝑆𝐹𝑜𝑀 §¨ are substantially higher than the mechanically-based solutions. This is mainly due to two reasons:

• The effective mass densities 𝜌 of the electrical solutions are more important than their mechanical counterparts. Indeed, with equivalent inertial mass, linear harvesters' volumes are usually smaller than nonlinear ones. • The energy harvester effectiveness criterion 𝐸 ª¨ is generally greater for electrical solutions.

However, the comparison with wideband mechanical solutions would not be fair without highlighting some limitations of this analysis:

This approach does not consider the power consumption of the control circuit (except for [START_REF] Cai | A Piezoelectric Energy-Harvesting Interface Circuit with Fully Autonomous Conjugate Impedance Matching, 156% Extended Bandwidth, and 0.38µW Power Consumption[END_REF]). Indeed, the control circuit as well as the 3-D MPPT loop that aims to control both 𝜙 " and 𝛥𝜙 should consume as few power as possible in order to preserve the efficiency of this approach and not waste all the harvested energy. Nonetheless, we are confident that designing this kind of lowpower interface is possible, since low-power self-powered SECE have already been designed in previous works [START_REF] Cai | A Piezoelectric Energy-Harvesting Interface Circuit with Fully Autonomous Conjugate Impedance Matching, 156% Extended Bandwidth, and 0.38µW Power Consumption[END_REF][START_REF] Hehn | A fully autonomous integrated interface circuit for piezoelectric harvesters[END_REF][START_REF] Morel | A Shock-Optimized SECE Integrated Circuit[END_REF]. 3-D MPPTs have also been implemented in prior art for DC/DC converters, exhibiting energy consumptions as low as 𝜇𝑊𝑠 [START_REF] Rawy | An 88% Efficiency 0.1-300-µW Energy Harvesting System with 3-D MPPT Using Switch Modulation for IoT Smart Nodes[END_REF].

• The electrical solutions have been proposed to harvest energy from mono-frequency vibrations, and to track the vibration's frequency with the harvester's resonant one. In the case of multi-frequency harmonic excitations, our circuit, because of its nonlinearity, would not exhibit the frequency response previously shown in the experimental part. Studying the behavior of the SC-SECE under wideband excitations is out of the scope of this paper but may be considered in future researches.

In another hand, we can also underline two main advantages the electrical solutions present and that are not pointed up by the 𝑆𝐹𝑜𝑀 §¨:

• The electrical solutions do not induce any high orbit, nor low orbit behaviors contrary to nonlinear mechanical oscillator for whom several more or less energy-generating behaviors can be observed for the same sinusoidal excitation. If rightly tuned, the electrical solutions will always extract the same amount of energy, which is predictable and hence offer more robustness than the nonlinear solution. • The electrical solutions bandwidth remains constant with any vibration's amplitude. This is not the case of the nonlinear mechanical solutions, whose frequency responses shapes are dependent on the excitation amplitude.

Conclusion

In this paper, we proposed a new electrical strategy called SC-SECE. This strategy allows to electrically tune both the resonant frequency and the induced damping of a highly coupled linear piezoelectric harvester. Thanks to the first harmonic analysis, we first proposed an analytical expression of the harvested power, and we determined numerically the power frequency responses of piezoelectric energy harvesters associated with the SC-SECE strategy. Next, we verified experimentally our analytical model on a fabricated highly coupled piezoelectric harvester made with PZN-PT plates. Under very low vibration's amplitude of 0.03G, we have been able to harvest more than 4𝜇𝑊 on a relatively large frequency band, between 90Hz and 140Hz. Increasing the acceleration up to 0.2G did not change much the reachable harvesting bandwidth, which remains of approximatively 50Hz. In a last section we compared the results obtained on this electrically tunable strategy with nonlinear harvesters. The reachable figure of merit seems promising, as long as the consumption of the interface circuit remains considerably smaller than the harvested power. Our future work will focus on the design of a dedicated ultra-low power ASIC realizing this strategy, in order to push forward the proposed concept, and turn it into a self-powered interface circuit.

Figure 2 -Figure 3 -

 23 Figure 2 -Voltage and mechanical displacement waveforms associated with a. the standard SECE and b. the proposed SC-SECE

Figures 4 ,

 4 5 and 6 present the SECE (dashed lines), SEH (dots), and SC-SECE (straight lines) power frequency responses for PEH exhibiting a low coupling (𝑘 F W = 0.01), medium coupling (𝑘 F W = 0.08), and high coupling (𝑘 F W = 0.4) respectively, for a given 𝑄 F = 20.

Figure 4 -Figure 5 -Figure 6 -

 456 Figure 4 -Power frequency responses of a lowly coupled PEH, with 𝑄 F = 20 and 𝑘 F W = 0.01

Figure 7 -

 7 Figure 7 -𝑃 " 𝑄 F as a function of 𝑘 F W and for 𝑄 F = 10 in blue, 𝑄 F = 20 in red, and 𝑄 F = 50 in yellow.

Figure 8 -

 8 Figure 8 -Product 𝑃 " 𝑄 F 𝛥𝛺/𝛺 -as a function of 𝑘 F W and for 𝑄 F = 10 in blue, 𝑄 F = 20 in red, and 𝑄 F = 50 in yellow.

Figure 9 -Figure 10 -

 910 Figure 9 -Experimental setup including a highly coupled PEH

Figure 11 -

 11 Figure 11 -Experimental waveforms of the proposed SC-SECE with various (𝜙 " , 𝛥𝜙) combinations: a. -¡ , W b. 0, ¡ ¢ c. 0, W £ d.

Figure 12 -

 12 Figure 12 -Experimental (crosses) and theoretical (straight lines) power frequency response of the standard SECE (green) and of the proposed SC-SECE (black) under small vibrations, with 𝛾 = 0.03𝐺, and optimal couples (𝜙 " , 𝛥𝜙).

Figure 14 -

 14 Figure 14 -Experimental (crosses) and theoretical (straight lines) power frequency response of the proposed SC-SECE under small vibrations, with 𝛾 = 0.2𝐺, and optimal couples (𝜙 " , 𝛥𝜙).

Figure 15 -

 15 Figure 15 -Experimental variations of 𝑃 F'' , 𝑄 F , 𝛥𝛺/𝛺 , , and 𝑓 ", with the acceleration amplitude 𝛾

5. 1

 1 Figure of merit classifying wideband harvesting solutions W. Liu et al. proposed a figure of merit called Systematic Figure of Merit with Bandwidth information (𝑆𝐹𝑜𝑀 §¨) in order to classify and categorize wideband piezoelectric harvesters [30]. The expression of this Figure of Merit (FoM) is reminded below:

Table 1 .

 1 Highly coupled PEH characteristics

	Parameters names	Parameters values	Units
	𝑘 W	0.32	-
	𝑘 F W	0.48	-
	𝑄 F	20	-
	𝐶 0	1.5	nF

Acknowledgments

This work was supported by the French Alternative Energies and Atomic Energy Commission.


Antoine Charpentier

David Mignon

Sophie Barbe

Juan Cortés

Thomas Schiex
email: thomas.schiex@inra.fr

Thomas Simonson

David Allouche
email: david.allouche@inra.fr

Juan Cortes

†miat

Variable Neighborhood Search with Cost Function Networks To Solve Large Computational Protein Design Problems

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Computational protein design (CPD) is an important developing tool in biotechnology, with possible applications in drug design, green chemistry, biomaterials, bioenergy and beyond. It is also an attractive approach to understand protein structure-function relationships.

The method starts from a particular protein sequence and structure and amino acid positions are iteratively selected and allowed to mutate so as to optimize a desired property, such as protein stability or ligand binding. In recent years, CPD has led to proteins with optimized thermostability and solubility, [START_REF] Korkegian | Computational thermostabilization of an enzyme[END_REF][START_REF] Diaz | Computational design and selections for an engineered, thermostable terpene synthase[END_REF] novel ligand-binding, [START_REF] Rudicell | Enhanced Potency of a Broadly Neutralizing HIV-1 Antibody In Vitro Improves Protection against Lentiviral Infection In Vivo[END_REF][START_REF] Roberts | Computational Design of a PDZ Domain Peptide Inhibitor that Rescues CFTR Activity[END_REF][START_REF] Verges | Computer-aided engineering of a transglycosylase for the glucosylation of an unnatural disaccharide of relevance for bacterial antigen synthesis[END_REF] novel enzyme activity, [START_REF] Smith | Novel proteins: from fold to function[END_REF][START_REF] Jiang | De novo computational design of retro-aldol enzymes[END_REF][START_REF] Khare | Computational redesign of a mononuclear zinc metalloenzyme for organophosphate hydrolysis[END_REF] and completely redesigned sequences. [START_REF] Kuhlman | Design of a novel globular protein fold with atomic-level accuracy[END_REF] CPD methods are characterized by an energy function, a description of conformational space, and a search method to explore sequences and conformations. If many amino acid positions are allowed to mutate, the problem becomes enormous, and so efficient energy functions and powerful search methods are required. Many applications rely on a molecular mechanics energy function and a simple implicit solvent model, such that the energy function takes the form of a sum over amino acid pairs. While recent works try to better represent the flexibility of proteins at the cost of increased computational complexity, [START_REF] Hallen | Dead-end elimination with perturbations (DEEPer): A provable protein design algorithm with continuous sidechain and backbone flexibility[END_REF] conformational space is usually defined by a fixed protein backbone and a discrete library of rotamers for each amino acid side chain. The unfolded state is usually treated implicitly, through a set of amino acid chemical potentials. Even with these drastic simplifications, the space to explore can be huge and the optimization problem remains NP-hard. Thus, with n mutating positions, twenty allowed amino acid types, an average of 12 rotamers per side chain type as in the Tuffery rotamer library, 11 and a total of N amino acids, the number of possible states is 12 N × 20 n .

The goal of CPD is usually to identify a small collection of low energy solutions, which can be experimentally tested. For some applications, it is important to sample large numbers of low energy states, or to enumerate all the states within a given energy window. This allows thermodynamic properties to be estimated, such as ligand binding free energies or acid/base constants. This can be achieved using ε-guaranteed algorithms such as TopN 12 and K * , [START_REF] Roberts | Computational Design of a PDZ Domain Peptide Inhibitor that Rescues CFTR Activity[END_REF]13 later improved with CFN algorithms [14][15][16] or extended to affinity-based sequence optimization. 17 Or it can be achieved using Monte Carlo based approaches [18][START_REF] Villa | Adaptive landscape flattening in amino acid sequence space for the computational design of protein:peptide binding[END_REF][START_REF] Polydorides | Monte Carlo simulations of proteins at constant pH with generalized Born solvent, flexible sidechains, and an effective dielectric boundary[END_REF][START_REF] Villa | Comparing pairwise-additive and manybody generalized Born models for acid/base calculations and protein design[END_REF] with only asymptotic guarantees in the limit of infinite sampling. For other applications, a few solutions are sufficient, which can then be refined or modified manually using expert knowledge and/or experiments. To evaluate the quality of sampled sequences and conformations, an important goal is to identify, when possible, the global minimum energy sequence and conformation, referred to as the GMEC. The GMEC is an important candidate solution in its own right, and it also serves as a reference, indicating the energy range that should be sampled by good solutions.

With a given energy function and conformational space, the quantity and quality of sampled states are determined by the exploration method. Several classes of methods exist. Simulated annealing with Monte Carlo (MC) is very common; by performing multiple runs one expects to sample a collection of low energy solutions, and to reach the GMEC in the limit of a very long run. [START_REF] Van Laarhoven | Simulated annealing: Theory and applications[END_REF] However, in practical situations, there are no guarantees and few theoretical results on convergence and the nature of the sampled ensemble. Related methods with similar properties include genetic algorithms [START_REF] Chowdry | An object-oriented library for computational protein design[END_REF] and a steepest descent heuristic (SDH [START_REF] Wernisch | Automatic protein design with all atom force fields by exact and heuristic optimization[END_REF][START_REF] Schmidt Am Busch | Computational protein design: software implementation, parameter optimization, and performance of a simple model[END_REF][START_REF] Simonson | Computational protein design: The proteus software and selected applications[END_REF] ) with multiple restarts. A second class of methods uses MC or Replica Exchange MC (REMC). [START_REF] Villa | Adaptive landscape flattening in amino acid sequence space for the computational design of protein:peptide binding[END_REF][START_REF] Yang | Computational methods for protein design and protein sequence variability: biased Monte Carlo and replica exchange[END_REF][START_REF] Druart | Protein:ligand binding free energies: a stringent test for computational protein design[END_REF][START_REF] Mignon | Comparing three stochastic search algorithms for computational protein design: Monte Carlo, Replica Exchange Monte Carlo, and a multistart, steepestdescent heuristic[END_REF] These methods efficiently sample very large numbers of states according to a Boltzmann distribution, which provides rigorous thermodynamic properties. Although convergence to the Boltzmann distribution is only guaranteed in the infinite limit, [START_REF] Grimmett | Probability and random processes[END_REF] it is common to observe numerical convergence over long runs for many smalland medium-sized problems. Assuming convergence, one then obtains tight confidence intervals for average properties such as rotamer probabilities, binding free energies or sequence profiles.

A third class of exploration methods aims to provably identify the GMEC and possibly states within a small energy window above it. Cost function networks [START_REF] Allouche | Computational protein design as a cost function network optimization problem[END_REF][START_REF] Traoré | A new framework for computational protein design through cost function network optimization[END_REF][START_REF] Allouche | Computational protein design as an optimization problem[END_REF][START_REF] Traoré | Fast search algorithms for computational protein design[END_REF] are one framework that can provably find the GMEC. However, the problem being NP-hard, success is not guaranteed within a reasonable computing time. In practice, state-of-the-art CFN algorithms use branch-and-bound search and tight incremental bounds to solve large design problems. The search is separated in two phases: in a first phase, increasingly good solutions are found until the optimal solution is found. In a second phase, search is proving that this solution is optimal (no better solution exists).

When the search is interrupted because of a cpu-time limit, one usually says that optimality was not proven: for sure the last phase was not finished (possibly not even started). With the Rosetta energy function [START_REF] O'meara | A Combined Covalent-Electrostatic Model of Hydrogen Bonding Improves Structure Prediction with Rosetta[END_REF][START_REF] Alford | he Rosetta All-Atom Energy Function for Macromolecular Modeling and Design[END_REF] and the Dunbrack rotamer library, [START_REF] Shapovalov | A Smoothed Backbone-Dependent Rotamer Library for Proteins Derived from Adaptive Kernel Density Estimates and Regressions[END_REF] the full design problem, where all amino acids were allowed to mutate, could be solved for proteins of about 100 amino acids. [START_REF] Simoncini | Guaranteed Discrete Energy Optimization on Large Protein Design Problems[END_REF] Comparable performances were obtained with a molecular force field and an EEF1 solvent model. [START_REF] Traoré | A new framework for computational protein design through cost function network optimization[END_REF] In recent experiments with the same force field but a different solvent model and rotamer library, the GMEC was usually found and proven for problems where 20 positions could mutate and up to 96 other could change rotamers. [START_REF] Mignon | Comparing three stochastic search algorithms for computational protein design: Monte Carlo, Replica Exchange Monte Carlo, and a multistart, steepestdescent heuristic[END_REF] For larger problems, CFN algorithms were not able to prove optimality.

In this work, we propose a new approach for large protein design problems that combines the capacity of CFN algorithms with a so-called Variable Neighborhood Search (VNS): a stochastic optimization algorithm originally defined in Operations Research. [START_REF] Mladenović | Variable Neighborhood Search[END_REF] The idea is to select a subset of amino acid positions and identify the GMEC for these using CFN algorithms, with the others held fixed. A new subset is then chosen, and so on, until the energy cannot be improved in a reasonable time. We will refer to the subset of positions as the mutator subset or simply the mutators. Given a current selection of mutators, its neighborhood designates its design space: the collection of designs that can be reached by changing the amino acid identity or rotamer of one or more mutators. The mutator optimizations are done with CFN. With a mutator number of one or two, VNS can be seen as MC in the low temperature limit. We will typically use larger mutator numbers, and will compare several methods to choose the mutators, which can be deterministic, probabilistic, or structure-based.

Below, we describe the CPD context and the VNS procedure, including the CFN framework and the VNS search. Next, we evaluate several variants of the VNS method on nine small proteins from the SH3, SH2, and PDZ families. The VNS variants are compared to each other and to results obtained here and earlier with SDH, MC, and REMC. [START_REF] Mignon | Comparing three stochastic search algorithms for computational protein design: Monte Carlo, Replica Exchange Monte Carlo, and a multistart, steepestdescent heuristic[END_REF] Finally, we discuss the potential impact of the method in CPD.

Theoretical framework

The CPD problem

In this work, we assume a "pairwise decomposable" energy function, [START_REF] Samish | The Framework of Computational Protein Design[END_REF] a fixed protein backbone, a discrete library of amino acid side chain rotamers. The pairwise energy function describes the energy of the molecular system of interest as a sum of terms, each of which depends on the relative positions of at most two atoms or atom groups. In the CPD case, this reduces to pairs of side-chain conformations. In this context, CPD typically searches for sequences that optimize the folding energy, namely the difference between the energies of the folded and unfolded states. The unfolded state is usually not explicitly modeled but captured by a sum of fitted terms, each of which depends on the amino acid type at one position. With these assumptions, a sequence-conformation is defined by a vector r whose components r i are integer numbers that indicate the side chain type and rotamer at amino acid position i. The folding energy of a sequence-conformation can be written as:

E(r) = E bb + i E(r i ) + i<j E(r i , r j ) (1)
E bb is the energy of the protein backbone. E(r i ) describes intra side chain and side chain-backbone interactions, and contains a contribution from the unfolded state. E(r i , r j ) is a side chain-side chain interaction. All these terms can be pre-computed and stored in an energy matrix. In some applications, the backbone term E bb may also include contributions from side chains that are not allowed to move or mutate.

Cost Function Networks

Cost function networks were introduced in Constraint Programming for discrete optimization. The problem is to minimize a cost function that depends on many discrete variables, [START_REF] Meseguer | Handbook of Constraint Programming[END_REF] a problem known as Weighted Constraint Satisfaction. The total cost is a sum of simpler terms, also called cost functions, that each depend on just a few variables. The problem of finding the GMEC can be reduced to this problem, as shown by the form of the CPD energy, Eq. ( 1). The CFN variables represent side chain types and rotamers. The CFN cost functions represent energy contributions, and the total cost function represents the folding energy. Formally, using the terminology of protein design, a CFN can be defined as follows:

Definition 1. A CFN (R, W
) is defined by:

• a set R of rotamer variables r i , where i is an amino acid position; i ∈ {1, . . . , n} def = I. Each r i takes its values in a finite domain D i (a list of side chain types and rotamers).

• a set W of cost functions w S , each involving a set of rotamer variables {r i |i ∈ S}, S ⊆ I.

We denote D S = i∈S D i the Cartesian product of the domains of all the variables indexed in S ⊆ I. A vector of D I represents an assignment of all the variables (types and conformation), in other words a sequence-conformation. It will be denoted r I or simply r. A vector of D S will be denoted r S .

The total cost w(r) of a sequence-conformation r is defined as the sum over all cost functions:

w(r) = W w S (r S ) (2)
Here, r S is the projection of r onto S, while w(r) could also be written w I (r I ). The total cost can be identified with the CPD folding energy. A cost function that depends on a single variable, say w {i} , is called a unary cost function; it corresponds to a term E(r i ) in the folding energy. A cost function that depends on two variables r i , r j is called a binary function; it corresponds to a term E(r i , r j ). The backbone energy E bb can be viewed as a cost function that has no variabledependence or equivalently depends on an empty set of variables: E bb ≡ w ∅ . The optimization problem is to identify a sequence-conformation r that gives the minimum total cost.

Guaranteed CFN solving

The Weighted Constraint Satisfaction problem can be solved using various algorithms, the most usual being a guaranteed branch-and-bound algorithm [START_REF] Allouche | Anytime hybrid best-first search with tree decomposition for weighted CSP[END_REF][START_REF] Hurley | Multi-Language Evaluation of Exact Solvers in Graphical Model Discrete Optimization[END_REF] that combines best and depth-first searches with tight incremental lower bounds. These algorithms rely on the representation of the space D I of sequence-conformations as a binary tree. The root of the tree is associated with the original CFN. Any node in the tree has two children obtained by choosing both a variable r i and a value ρ in its domain. The left child is a new CFN with r i set to ρ (D i becomes {ρ}). The right child is a CFN where ρ has been removed (D i becomes D i \{ρ}). The leaves of this tree correspond to CFNs with all variables assigned to a distinct state r I ∈ D I . Together, they represent the complete set D I and their energies. Because this tree has guaranteed exponential size, it is never explored fully. Instead, the algorithm relies on the energy of the best solution (leaf) found so far as a global upper bound on the energy of the GMEC. It computes a lower bound on the energy of the CFN associated with the current node using so-called "local consistency-enforcing algorithms". [START_REF] Cooper | Soft arc consistency revisited[END_REF] These algorithms reformulate a CFN in an equivalent CFN that gives direct access to a non-naive lower bound on the optimum energy. If at a given node, this local lower bound is equal to or larger than the global upper bound, the exploration of the branch is stopped. Depth-first search will then backtrack and explore the most recent unexplored branch, with the guarantee that the optimal solution cannot be missed. Figure 1 shows a complete binary tree of depth 4.

Partial tree search with limited discrepancies

The amino acid position i and the rotamer value ρ chosen for testing at each node determine the precise tree structure. In practice, the branch-and-bound algorithm spends only a minor fraction of its time to find the optimal leaf and the rest to prove optimality. Because we are not focused on guaranteed optimization, we used a partial-tree search mechanism called "Limited Discrepancy

Search" (LDS) [START_REF] Harvey | Limited discrepancy search[END_REF] which we now rapidly describe: LDS assumes that for every position i, a rotamer value ρ 0 ∈ D i that most likely participates in a low energy solution can be determined heuristically.

A branch in the tree that sets r i to a rotamer different from ρ 0 is called a discrepancy. LDS is parameterized by a bound δ on the maximum number of discrepancies that are allowed in any explored branch. With δ = 0, a single branch defined by the heuristically chosen rotamers is explored. With larger δ bounds, an increasingly large fraction of the tree can be explored. When δ = n (the number of positions in the protein sequence), all branches can potentially be explored.

To obtain good solutions, Fontaine et al. [START_REF] Fontaine | Exploiting Tree Decomposition for Guiding Neighborhoods Exploration for VNS[END_REF] proposed setting the number δ of allowed discrepancies to three or less. This gave a good balance between computational cost and quality on several large benchmark problems. The size of the partial tree that can be explored scales as [START_REF] Larrosa | Limited Discrepancy AND/OR Search and Its Application to Optimization Tasks in Graphical Models[END_REF] where k is the number of mutators. In practice, because of the tight lower bounds provided by CFN local consistencies, only a tiny unpredictable fraction of this is explored. Beyond the lower bound, local consistencies also provide, for free, the preferred rotamer ρ 0 ∈ D i that is needed by LDS: this is the so-called "support" as defined by Cooper et al., [START_REF] Cooper | Soft arc consistency revisited[END_REF] in Definition 4.5. It is at least guaranteed to have a minimum one-body energy.

(k + δ) δ+1 ,

Variable Neighborhood Search overview

The fundamental mechanism driving VNS is a local search within the neighborhood of a given state, or sequence-conformation r ∈ D I . Neighboring states are ones obtained by changing one or more types or rotamers of r. The amino acid positions allowed to mutate are the mutators. If the best state r within the neighborhood has a better energy than r, it is adopted and the process is repeated. Such a local search will eventually be trapped in a local minimum. VNS then restarts, using a new, larger set of mutators and thus a larger neighborhood.

VNS iteratively performs three steps. In the first step, a subset K of k(i) amino acid positions are selected randomly to be mutators. The number k(i) is drawn from a sequence of integers, in-Figure 1: A binary search tree with four amino acid positions with two possible rotamer values each. At the root node, an unassigned position (r 1 ) is chosen and either assigned a chosen rotamer (left branch) or this rotamer removed from the domain (right branch). If one uses δ = 2 in LDS, the dashed branches on the rightmost parts of the tree and subtrees can never be explored. For δ = 0, the only branch explored always follows left branches, defining an initial "greedy" candidate solution, which energy already provides an upper bound on the GMEC energy.

dexed by i, which depends on the current iteration number. For example, k(i) could be drawn from a Luby sequence, defined below. In the second step, a search in the corresponding neighborhood yields a new solution r . In the third step, if r has a lower cost, it replaces r and the number of mutators is reset to its initial value k(0); otherwise, i is incremented and the next mutator number in the sequence will be used. The algorithm is shown in Fig. 2. For neighborhood exploration, we use a partial LDS, depth-first, branch-and-bound algorithm that exploits local consistencies for both energy bounds and the heuristic choice of rotamers.

Neighborhood growth

At the end of each VNS iteration, the number of mutators is modified according to the sequence k(i). If the latest iteration was successful (energy improvement), i is reset to 0 and the algorithm shifts back to the beginning of the mutator number sequence. If not, i is increased. The simplest mutator number sequence is defined by k(i + 1) = k(i) + 1. This is called the k ++ sequence.

Because LDS complexity can rapidly increase with the mutator number, it seems more suitable Otherwise, the next mutator number in the sequence will be used.

Algorithm 2:

The VNS schema: an initial sequence-conformation r is generated by a greedy assignment of variables (line 1). VNS iterates until the mutator number k reaches n or the runtime exceeds a limit (line 3). k(i) mutators are selected (line 4). The corresponding rotamer variables are optimized (the other rotamers being fixed) using limited discrepancy search. If the energy was improved, we shift back to the start of the mutator number sequence; otherwise, the next mutator number will be used. If i + 1 is a power of 2, and i = 2 t -1, then luby(i) = 2 t-1

Function VNS/LDS+CFN(R, W, δ) 1 r ← genInitSol() // initialize sequence-conformation 2 i ← 0 // i: progress in the mutator sequence 3 while (k ≤ n) ∧ (not T imeOut) do 4 select a subset K ⊂ I of k(i) mutators //
If not, and

2 t-1 ≤ i < 2 t -1, then luby(i) = luby(i -2 t-1 + 1)
where t is the solution of the equation in the left part of either of the conditions above. This expression defines the following sequence for luby(i): {1,1,2,1,1,2,4,1,1,2,1,1,2,4,8,1,. . .} In this series, powers of 2 are spaced increasingly and irregularly. We then defined a mutator sequence k m (i), parameterized by its initial value k m (0) and a constant integer m. If luby(i) = 2 m , we increase the mutator number by one; otherwise, we keep the previous mutator number:

If luby(i) = 2 m , k(i) = k(i -1) + 1 Otherwise, k(i) = k(i -1)
With the smallest m value (m = 0), the increase of the mutator number is maximum, but it is already slower than the k ++ rule, with a sequence of increments {1,1,0,1,1,0,0,1,1,0,1,1,0,0,0,1,. . .}. With a larger m, the mutator number increases even more slowly. With a slower increase, the probability is higher that any amino acid position will be considered as a mutator before larger neighborhoods are considered. Thus, local energy improvements are more likely to be found earlier.

Structure-based versus random mutator sets

The default VNS approach draws the k mutator positions from a uniform random distribution. In the CPD context, we expect that the protein 3D structure can guide a more efficient heuristic. The first mutator was always taken from a uniform distribution. The remaining mutators were chosen using a deterministic or a probabilistic scheme. The deterministic scheme used the nearest neighbors of the first mutator, based on the distances between side chain C β atoms. In the probabilistic schemes, the remaining mutators were sampled either uniformly at random or with a probability p i that decreased with distance in the 3D protein structure:

p i ∝ 1 d i (3)
where d i is the distance of position i from the first mutator. This last heuristic will be denoted d-prob.

VNS heuristics

By combining a specific mutator number sequence with neighborhood growth heuristics, we obtained six different VNS heuristics. The mutator number sequence was given either by the k ++ or the Luby-based sequences. The neighborhood growth was done either with a uniform random method, a nearest neighbor method, or a d-prob method. The various heuristics are represented as leaves of the tree in Figure 3.

k ++ random d-prob k-nearest k m (i) (Luby) random d-prob k-nearest

Numerical methods

Models and test cases

To evaluate VNS in the CPD context, we used test problems from earlier work. [START_REF] Mignon | Comparing three stochastic search algorithms for computational protein design: Monte Carlo, Replica Exchange Monte Carlo, and a multistart, steepestdescent heuristic[END_REF] They were based on nine protein domains, from the SH3, SH2 and PDZ structural families in Table 1. For each domain, the energy matrix was computed with the Proteus CPD package, [START_REF] Simonson | Computational protein design: The proteus software and selected applications[END_REF] using the Amber ff99SB molecular mechanics energy function, a Coulomb Accessible Surface Area (CASA) implicit solvent model [START_REF] Mignon | Comparing three stochastic search algorithms for computational protein design: Monte Carlo, Replica Exchange Monte Carlo, and a multistart, steepestdescent heuristic[END_REF][START_REF] Lopes | Computational sidechain placement and protein mutagenesis with implicit solvent models[END_REF] and the Tuffery rotamer library. [START_REF] Tuffery | A New Approach to the Rapid Determination of Protein Side Chain Conformations[END_REF] The model represents a good balance between accuracy and computational effort. [START_REF] Gaillard | Full Protein Sequence Redesign with an MMGBSA Energy Function[END_REF] The first nine test problems (Full tests) corresponded to the full redesign of each protein, except for GLY and PRO residues, which were not mutated. Additional test problems corresponded to the redesign of 20 or 30 amino acid positions in each protein (N20 and N30 tests, respectively). Five choices of mutating positions were considered in each case.

The positions were close together in the 3D structure. amino acids that were not allowed to mutate kept their wild type amino acid types, but were allowed to explore all possible rotamers. The test set contained 9 × (5 + 5 + 1) = 99 problems in total. The search space size for each problem is the product of the total number of types and rotamers allowed at each amino acid position. These sizes ranged from 10 61 to 10

Numerical protocols

The six VNS heuristics were implemented in the toulbar2 CFN solver. The C++ code of toulbar2 is available under an MIT license at https://github.com/toulbar2/toulbar2. For each test problem, we applied VNS with the six heuristics described above (Fig. 3). The initial number of mutators was k(0) = 4 with a fixed discrepancy value of 3, a value set according to results previously obtained on a large collection of problems (not including these) in The Probabilistic Inference Challenge (PIC2011)(http://www.cs.huji.ac.il/project/PASCAL/). For the Luby based variants, the parameter m was set to 3, based on short test runs on 20 instances randomly picked in the FULL, N20 and N30 set of the 1CKA and 1ABO proteins (over 990 overall). We imposed a CPU time limit of two hours. To estimate variability, ten runs were done for each test problem, with different initial random seeds.

The results were compared to those obtained by three other search methods: [START_REF] Mignon | Comparing three stochastic search algorithms for computational protein design: Monte Carlo, Replica Exchange Monte Carlo, and a multistart, steepestdescent heuristic[END_REF] Monte Carlo (MC), Replica Exchange Monte Carlo (REMC), and a multi-start, steepest-descent minimization heuristic (SDH). The results for the MC runs and the N20 and N30 REMC runs were presented in previous work. The SDH and the Full REMC runs are reported here. The SDH method uses a large number of "minimization cycles". In each cycle, a random sequence-conformation is chosen, then the energy is minimized, one residue at a time, with the others fixed in their current state. The entire sequence is scanned repeatedly this way until the energy no longer improves, ending the cycle. [START_REF] Mignon | Comparing three stochastic search algorithms for computational protein design: Monte Carlo, Replica Exchange Monte Carlo, and a multistart, steepestdescent heuristic[END_REF] The new REMC runs used an improved selection of two-position moves (which provides better results than the selection strategy used earlier [START_REF] Mignon | Comparing three stochastic search algorithms for computational protein design: Monte Carlo, Replica Exchange Monte Carlo, and a multistart, steepestdescent heuristic[END_REF] ). The new SDH runs used a slightly modified treatment of histidine residues, allowing multiple protonation states. All these results were obtained with the Proteus software, [START_REF] Simonson | Computational protein design: The proteus software and selected applications[END_REF] using the following settings:

• MC: The simulation temperature was constant and corresponded to a thermal energy of kT = 0.2 kcal/mol. The simulations lasted 10 9 steps and never exceeded 9 hours.

• Steepest descent heuristic (SDH): Each run included 110,000 minimization cycles, except for a few cases that were extended to 330,000 or 990,000 cycles; run times did not exceed 24 hours.

• REMC: The REMC runs used eight replicas and 0.75 x 10 9 MC steps per replica. The thermal energies of the replicas were spaced in a geometric progression between 0.125 and 2 or 3 kcal/mol. The total CPU time per core was never more than 3 hours, for a total cpu-time on an 8-core machine of 24 hours at most. [START_REF] Mignon | Comparing three stochastic search algorithms for computational protein design: Monte Carlo, Replica Exchange Monte Carlo, and a multistart, steepestdescent heuristic[END_REF] For each test we computed the difference ∆E = E Proteus -E VNS between the best energy obtained with the Proteus protocols and the best VNS energy. A positive value indicates an improvement with VNS. For the nine Full tests, new REMC runs were performed here. The new runs used an improved selection of two-position moves (which provides better results than the selection strategy implemented in previous work [START_REF] Mignon | Comparing three stochastic search algorithms for computational protein design: Monte Carlo, Replica Exchange Monte Carlo, and a multistart, steepestdescent heuristic[END_REF] ).

Results

Qualitative analysis

For a given VNS heuristic and category of problems (Full, N30, or N20), Table 2 reports the frequency of improvement (∆E > 0) and stagnation (∆E = 0) over all runs and test problems.

All six VNS variants provided energies of the same or better quality than MC+REMC+SDH for most of the tests, especially for the larger tests. On the easiest category (N20), the energy was unchanged for about 2/3 of the tests and improved for about 30%, with a mean improvement of about 0.1 kcal/mol. Indeed, for most of the N20 tests, MC+REMC+SDH found the GMEC, [START_REF] Mignon | Comparing three stochastic search algorithms for computational protein design: Monte Carlo, Replica Exchange Monte Carlo, and a multistart, steepestdescent heuristic[END_REF] which cannot be improved. On the most difficult, Full category, the best VNS methods improved over MC+REMC+SDH in 94-99% of the cases, with a mean energy improvement of about 3±1 kcal/mol. The best VNS energies, overall, were given by the Luby/d-prob heuristic. For the N20 set, we used toulbar2 to identify proven GMECs using options -dee: -d:

-l=3 -m=2. [START_REF] Traoré | A new framework for computational protein design through cost function network optimization[END_REF] Using the best solution provided by the VNS heuristics as an initial upper bound, toulbar2 was able to prove the optimality of 36 out of the 45 instances within a 72 hours timeout but never improved over the best VNS solution. Over this set of 36 instances with a proven GMEC, VNS with Luby/d-prob found the GMEC in 357 of 360 runs (99.1%). For the 3 remaining designs, the best solution could not be improved and no optimality proof reached within the 72 hours timeout. SDH reached the GMEC in 25 of the 36 cases. When all N20 instances are considered, VNS improved over (SDH+MC+REMC) in 14 cases, with improvements ranging from 0.04 to 1.39 kcal/mol. While the average improvement was very small over N20, there were no runs with ∆E < 0.

VNS Luby/d-prob also gave the best results on N30, offering improved energies for almost 2/3 of the N30 tests, including nine improvements by more than 1.5 kcal/mol. The best energy improvement obtained over all runs and problems was 5.4 kcal/mol.

To see how the results of the best VNS (Luby/d-prob) distribute over the different test problems, we show a box-plot of the distribution of ∆E over the nine Full tests in Figure 5 (boxplots for the N20 and N30 datasets, on a per instance basis, are available as Figures S1 andS2). If we consider the best VNS runs (among the 10 runs per test), the energy improvements ranged from 1.3 to 11.2 kcal/mol. In 8 cases out of 9, the median energy obtained with VNS Luby/d-prob also improved over the best energy obtained by MC+SDH+REMC with a maximum of 11.2 kcal/mol for 1G90.

The worst median ∆E was obtained on 1A81, with a value -0.30 kcal/mol. In this case, both the average and maximum improvements remained positive (respectively, 0.2 and 2.5 kcal/mol).

Optimality proof for side chain positioning

VNS does not provide an optimality proof for its best solution, which may not be the GMEC. As a partial optimality test, we considered the Full tests and asked whether for the obtained sequence, the obtained side chain rotamers were indeed optimal. Using Luby/d-prob, for every test problem, an optimal side chain positioning was found among the 10 runs. In the MC+REMC+SDH tests, [START_REF] Mignon | Comparing three stochastic search algorithms for computational protein design: Monte Carlo, Replica Exchange Monte Carlo, and a multistart, steepestdescent heuristic[END_REF] the optimal side chain rotamers were obtained only for five of the nine Full instances.

Over all problem categories and all runs, Luby/d-prob appears to improve conformational sampling, since 91% of all lowest-energy sequences corresponded to optimal side chain positions. As another test of reliability, we restarted sequence optimization from the best VNS solutions using the Hybrid best first search algorithm 42 (toulbar2 command line with no options) and with a previously published Depth First Search protocol [START_REF] Traoré | A new framework for computational protein design through cost function network optimization[END_REF] (both available in toulbar2). After 24 hours timeout with each protocol, none of these did provide any energy improvement, suggesting that the VNS energies were already very good. No optimality proof could be obtained either.

Neighborhood growth speed

We compared the k ++ and Luby neighborhood growth sequences during the VNS searches. We considered the nine Full problems. The histogram of mutator numbers used in the VNS searches is shown in Figure 6. With k ++ , the mutator number increased rapidly to sizes that define very hard subproblems. Using Luby, with the parameter m = 3, neighborhood growth was slower and energy exploration more effective. The largest mutator sets included about 40 amino acids, which is fewer than the average number of residues within 14 Å of a given amino acid residue. With

Luby, seven times more subproblems (or nodes) were explored on average. The largest mutator numbers were 49 and 40, respectively, with k ++ and Luby. The median mutator number was 8 with both methods.

Computational efficiency

As shown in Figure 7, the average time by run is similar for all the VNS heuristics and is around 10, 15 and 50 minutes for respectively the N20, N30 and FULL cases. The use of the Lubybased heuristic has a limited effect on the run time. However, it does increase the probability of identifying the GMEC (without proving optimality). An indicative rate of best solution recovery is provided in Table S1.

With the Luby/d-prob strategy, and considering only "FULL" instances, Table 4 shows that the average time to obtain the best solution for the considered proteins ranges between 14 and 88 minutes. The observed variability is clearly related to the different sizes of the proteins and is very likely also influenced by the interaction terms defined by the energy functions for each of the considered proteins.

In comparison, [START_REF] Mignon | Comparing three stochastic search algorithms for computational protein design: Monte Carlo, Replica Exchange Monte Carlo, and a multistart, steepestdescent heuristic[END_REF] SDH and REMC methods maximum run times reached almost 24 hours of sequential cpu-time (REMC using at most 3 hours on 8 cores). MC maximum run time was close to nine hours.

Concluding discussion

Cost function networks and their algorithms, implemented in Toulbar2, define a powerful framework for large CPD problems, able to provably obtain the GMEC for large problems in short times. For other problems that are too large for the GMEC to be found, we would nevertheless like to exploit the power of CFNs. This is done here by applying CFN optimization to a series of subproblems, whose size is gradually increased until the energy does not improve further. This iterative approach is referred to as Variable Neighborhood Search, and has been extensively used for other discrete optimization problems. [START_REF] Hansen | Variable neighbourhood search: methods and applications[END_REF][START_REF] Jarboui | 3rd International Conference on Variable Neighborhood Search (VNS'14)[END_REF] The key parameters of the method are the sequence of subproblem sizes and the detailed method for choosing the subproblems. Here, we defined a subset of mutator positions and we increased their number with the help of a Luby integer sequence.

The best strategy for choosing the mutators was a random one (d-prob) that choses with a high probability positions close together in the 3D protein structure. This strategy appears to balance the need to distribute mutators throughout the structure and the need for concerted rearrangements of groups of nearby amino acids in order to optimize packing and lower the folding energy. In this respect, our VNS strategy contrasts with heuristic schemes like SDH, where individual positions are optimized, or Monte Carlo where one or two amino acids are changed at a time. The largest energy improvements in our tests occurred when there were 4-8 mutators. Thanks to the CFN optimization, these positions are rigorously optimized, and all relevant concerted rearrangements are taken into account.

The largest CPD problems where Toulbar2 identified the GMEC 38 used the Talaris force field distributed with Rosetta, where different distance cutoffs are applied to different energy terms, between 6 and 10 Å. The problems considered here used a larger energy cutoff, which presumably makes the energy surface more rugged and complex, hence the need for a complex, multi-position exploration. It may be that a simplified energy function with an aggressive cutoff could be used to guide the VNS method, for example to choose the detailed structure of the search tree. The full energy function would then be optimized with the help of the resulting tree. For non-globular proteins that have an elongated structure, a smaller energy cutoff could make the problem decomposable into subproblems that interact only through few amino acids. The VNS heuristics developed in this paper could be extended to exploit such situations.

Because of its increased efficiency, our method can be highly useful in various CPD contexts.

The benchmarked problems used Tuffery's rotamer library [START_REF] Tuffery | A New Approach to the Rapid Determination of Protein Side Chain Conformations[END_REF] with 248 rotamers for a fully mutable residue. Since intensive CFN algorithms are only applied on relatively small subproblems in VNS, there is no visible reason that would preclude the use of VNS and the proposed heuristic with denser rotamer libraries [START_REF] Shapovalov | A Smoothed Backbone-Dependent Rotamer Library for Proteins Derived from Adaptive Kernel Density Estimates and Regressions[END_REF] or energy-based conformer libraries. [START_REF] Subramaniam | Backbone dependency further improves side chain prediction efficiency in the Energy-based Conformer Library (bEBL)[END_REF] More generally, whether it is for a better modeling of the protein flexibility or to account for multiple desirable or undesirable states of the protein being designed, [START_REF] Davey | Multistate approaches in computational protein design[END_REF][START_REF] St-Jacques | Modern Biocatalysis: Advances Towards Synthetic Biological Systems[END_REF] computational design will require to process increasingly harder and larger problems. In this context, our approach could be highly useful to handle the greatly increased search space of such refined models. This is also supported by the fact that VNS was recently parallelized. [START_REF] Ouali | Replicated parallel strategies for decomposition guided VNS[END_REF] The VNS algorithm is not designed for the extensive enumeration of low-energy states, unlike Monte Carlo. Nevertheless, the visited sequences give information about local minima and the energy landscape of the system, through the collection of suboptimal solutions explored during the search. Another limitation of the VNS method is that the folding energy is optimized without considering alternative objectives, such as protein solubility or specific binding of other proteins.

This defines multi-objective optimization problems, which are computationally difficult and often addressed with simple, heuristic methods. For example, to design ligand binding while maintaining protein stability, it is common to optimize the energy of the protein-ligand complex. A more rigorous approach would optimize the ligand binding free energy while imposing a limit on the folding energy. Such an approach was discovered and applied recently, with the help of adaptive importance sampling Monte Carlo. [START_REF] Villa | Comparing pairwise-additive and manybody generalized Born models for acid/base calculations and protein design[END_REF] Another example of a multi-objective optimization would be to optimize folding energy while also optimizing the folding of the associated messenger RNA, or while enforcing a DNA coding sequence that overlaps with another protein coding sequence. [START_REF] Opuu | Computational design of fully overlapping coding schemes for protein pairs and triplets[END_REF] The VNS framework offers possibilities along these lines, since it can take advantage of "constraint programming expressivity" 41 in order to enrich the CPD model with new types of constraints.

In conclusion, we have adapted the Variable Neighborhood Search algorithm to large protein design problems. We have shown that a slower growth of the mutator set and search space and a probabilistic mutator selection significantly improve the performance. Complex, concerted structure rearrangements are explored for each mutator thanks to the CFN framework. As CPD models become more realistic and complex, finding the GMEC and proving its nature will become increasingly intractable. In the context of very large problems and/or more complex physical models, the search approach described in this paper should be highly useful.

Figure 2 :

 2 Figure2: Variable Neighborhood Search procedure with mutator sequence k(i). The protein backbone is shown on the left as a chain of gray beads, representing amino acid residues. In step (1), we choose a subset of k(i) mutators (blue positions in the middle). In step (2), their amino acid types and rotamers are optimized, with the other positions fixed. The new types/rotamers are colored on the right. In step (3), if the energy has improved, the next iteration will use a mutator number k(0). Otherwise, the next mutator number in the sequence will be used.

Figure 3 :

 3 Figure 3: Representation of six different VNS heuristics defined by the mutator number sequence used and the neighborhood growth method used.

Figure 4 :

 4 Figure 4: Average and standard deviation of the average energy improvement using different VNS heuristics and test categories. The average and standard deviation for a VNS method H and a test category C (Full, N30, N20) were computed over all the tests in C. Blue and green colors are respectively associated with k ++ and Luby mutator growth methods. Squares, circles and triangles are respectively associated with k-nearest, random, and d-prob position selection. Each category is highlighted by an ellipsoid.

Figure 5 :

 5 Figure 5: Box plot of ∆E values with VNS Luby/d-prob from ten runs on each Full test.

Figure 6 :

 6 Figure 6: Histogram of the mutator numbers explored during the Full design tests using the k ++ (blue) or Luby (green) sequences.

Figure 7 :

 7 Figure 7: Average run times in minutes on FULL (blue circle), N30 (orange cross) and N20 (green triangle) for each considered heuristic.

( 11 )

 11 Gaillard, T.; Panel, N.; Simonson, T. Protein side chain conformation predictions with an MMGBSA energy function. Proteins: Struct., Funct., Bioinf. 2016, 84, 803-819. (12) Henrion, M. Search-based methods to bound diagnostic probabilities in very large belief nets. Seventh conference on Uncertainty in Artificial Intelligence. 1991; pp 142-150. (13) Georgiev, I.; Lilien, R. H.; Donald, B. R. The minimized dead-end elimination criterion and its application to protein redesign in a hybrid scoring and search algorithm for computing partition functions over molecular ensembles. J. Comput. Chem. 29, 1527-1542. (14) Viricel, C.; SIMONCINI, D.; Allouche, D.; De Givry, S.; Barbe, S.; Schiex, T. Approximate counting with deterministic guarantees for affinity computation. International Conference on Modelling, Computation and Optimization in Information Systems 2015. 2015; pp 165-176, Proceedings of the 3rd International Conference on Modelling, Computation and Optimization in Information Systems and Management Sciences.(15) Viricel, C.; Simoncini, D.; Barbe, S.; Schiex, T. Guaranteed weighted counting for affinity computation: Beyond determinism and structure. International Conference on Principles and Practice of Constraint Programming. 2016; pp 733-750. (16) Viricel, C.; de Givry, S.; Schiex, T.; Barbe, S. Cost function network-based design of protein-protein interactions: predicting changes in binding affinity. Bioinformatics 2018, 733-750. (17) Ojewole, A. A.; Jou, J. D.; Fowler, V. G.; Donald, B. R. BBK* (Branch and Bound Over K*): A Provable and Efficient Ensemble-Based Protein Design Algorithm to Optimize Stability and Binding Affinity Over Large Sequence Spaces. J. Comput. Biol. 2018, 25, 726-739, PMID: 29641249. (18) Druart, K.; Bigot, J.; Audit, E.; Simonson, T. A hybrid Monte Carlo method for multibackbone protein design. J. Chem. Theory Comput. 2017, 12, 6035-6048.

  their rotamers are r K 5 use LDS+CFN.Search to select new rotamers r K // The sequence-conformation is now r = r K + r I\K to only increase it if smaller neighborhoods have been sufficiently explored. The simplest way to achieve this would be to increase the mutator number less frequently. With this aim, we have explored a new update rule based on the Luby series,[START_REF] Luby | Optimal speedup of Las Vegas algorithms[END_REF] defined by:

	6	if w(r ) < w(r) then
	7	r ← r
	8	i ← 0	// Shift back to initial mutator number

9 else

i ← i + 1 // next mutator number 10 return r

Table 1 :

 1 229 . Protein set for testing.

	Type	Name	PDB Length Vars. a
	PDZ NHERF 1G9O	91	76
	PDZ Syntenin 1R6J	82	72
	PDZ	DLG2	2BYG	97	82
	SH3	Abl	1ABO	56	48
	SH3	Crk	1CKA	57	49
	SH2	Syk	1A81	108	92
	SH2	Grb2	1BM2	98	88
	SH2	Zap70	1M61	109	88
	SH2	Src	1O4C	104	96
	a The number of amino acids, excluding Gly and Pro.

Table 2 :

 2 Frequency of energy improvement with VNS: each row represents a VNS heuristic and reports the percentage of runs in which the energy improved (+) or was unchanged (0) compared to SDH+MC+REMC.For each VNS heuristic H and test problem P , we averaged the energy gain over the ten runs; the result is denoted ∆E(H, P ). The standard deviation over the runs is denoted σ(∆E)(H, P ):For each heuristic H and problem category C, we computed the average of ∆E(H, P ) and σ(∆E)(H, P ) over the corresponding tests, denoted ∆E C and σ(∆E) C . The results are given in Table3and shown as a scatter plot in Figure4.

		Full		N30		N20
	VNS	+	0	+	0	+	0
	Luby/d-prob 98% 0% 65% 31% 31% 68%
	Luby/random 99% 0% 63% 31% 31% 69%
	Luby/nearest 94% 0% 60% 31% 28% 66%
	k ++ /d-prob	80% 0% 56% 27% 27% 66%
	k ++ /nearest	80% 0% 54% 21% 21% 62%
	k ++ /random 76% 0% 56% 25% 25% 67%

Table 3 :

 3 Averages ∆E C of the mean energy improvement, compared to MC+REMC+SDH, over the test category C, and average σ(∆E) C of its standard deviation. -prob is the best method, offering both quality and reliability in most test categories and especially in the most challenging "Full" category. In this category, the results of k ++ and Luby are clearly separated. Luby outperform k ++ , providing better energies and smaller standard deviations, especially with the structure-based neighborhood growth heuristic d-prob.

			Full		N30		N20
	VNS	∆E C σ(∆E) C	∆E C σ(∆E) C	∆E C σ(∆E) C
	Luby/d-prob	3.17	1.03	0.77	0.35	0.15	0.01
	Luby/random	2.95	1.06	0.76	0.32	0.15	0.01
	Luby/nearest	2.78	1.36	0.69	0.33	0.13	0.07
	k ++ /d-prob	0.77	2.40	0.59	0.47	0.13	0.11
	k ++ /nearest	1.03	2.21	0.43	0.59	0.10	0.20
	k ++ /random	0.55	2.47	0.58	0.48	0.14	0.05
	Luby/d

Table 4 :

 4 Family, PDB identified, problem size (number of variables) and CPU times (in minutes) for the full redesign of each protein, averaged over ten runs using Luby d-prob.

	Family PDZ	SH2	SH2	SH2	SH2	PDZ	PDZ	SH3	SH3
	PDB	1A81 2BYG 1O4C 1M61 1BM2 1R6J	1G9O 1ABO 1CKA
	#. vars	92	97	96	88	88	72	76	48	49
	cpu	74±24 70±34 55±14 50±18 49±18 45±23 37±20 15±7 11±3

Acknowledgements

This work was granted access to the HPC resources on the TGCC-Curie supercomputer, the Computing mesocenter of Région Midi-Pyrénées (CALMIP, Toulouse, France) and the Geno-Toul (Toulouse) Bioinformatic platform. This work was supported by the EMERGENCE program of IDEX Toulouse (E-CODE project), the French National Research Agency (ANR-12-MONU-0015-03 for JC, TSi, SB, ANR-16-C40-0028 for TSc and DA) and the French National Institute for Agricultural Research (INRA).

Associated Contents

Supplementary materials are provided that describe more detailed information on the per-instance energy gap as well as the frequency of recovery of the best energy on each problem category. This information is available free of charge via the Internet at http://pubs.acs.org.