Vladimir P Zhdanov
email: zhdanov@catalysis.ru.

Bengt Kasemo

Protein adsorption and desorption on lipid bilayers

Keywords: Proteins, Lipids, Adsorption and desorption kinetics, Lipid redistribution, Electrostatic interactions -------

The protein surface usually exhibits one or a few charged spots. If a lipid bilayer contains a significant amount of lipids with oppositely charged head groups, protein adsorption on a bilayer may be energetically favourable due to the protein-lipid electrostatic interaction. The specifics of this case is that the lipids are highly mobile and the protein adsorption is accompanied by the redistribution of lipids between the areas covered and not covered by protein. We present a kinetic model illustrating that this effect is especially interesting if the fraction of the surface covered by charged lipids is relatively low. In this situation, with increasing protein coverage, the protein desorption rate constant rapidly increases while the adsorption rate constant drops, so that there is critical fraction of the area covered by protein.

Adsorption above this fraction is hindered both kinetically and thermodynamically.

A C C E P T E D M A N U S C R I P T 1 Introduction

Protein adsorption on solid surfaces has long attracted the attention of researchers working in the natural sciences, medicine and industry [START_REF]Proteins at Interfaces II[END_REF][START_REF] Hlady | Protein adsorption on solid surfaces[END_REF][START_REF]Biopolymers at Interfaces[END_REF][START_REF] Gray | The interaction of proteins with solid surfaces[END_REF][START_REF] Tsapikouni | Protein-material interactions: From micro-tonano scale[END_REF], because it plays an important role in biology and in numerous applications, e.g., for biosensors and protein chips, medical implants and in the food industry. The kinetics of protein adsorption are often apparently simple, but in reality very complex, due to multiple binding sites on the protein surface, many nearly iso-energetic conformational states of proteins, adsorption-induced conformational changes, diffusion limitations in the bulk phase, surface diffusion, ordering or aggregation of adsorbed proteins, surface roughness, and a multitude of related microscopic factors. Protein adsorption on soft matter in general and lipid bilayers in particular has attracted much less attention although it is also of interest from various perspectives. Maybe the most important one is protein interaction with cell membranes. In this case, specific binding sites are usually involved, like membrane proteins or peptides, but the interaction (nonspecific binding) with the plane membrane is also of interest, e.g., as a possible precursor for the specific binding. Another case where protein binding to lipid bilayers is of interest is biosensors based on supported lipid membranes. A unique new ingredient in adsorption on lipid membranes, compared to solid surfaces, is the fluidity of the bilayer which opens up for rapid adjustment of the lipids to optimize the interaction energy (see Ref. [START_REF] Flenner | Subdiffusion and lateral diffusion coefficient of lipid atoms and molecules in phospholipid bilayers[END_REF], focused on diffusion of lipids in bilayers, and more general related reviews [START_REF] Heimburg | Monte Carlo simulations of lipid bilayers and lipid protein interactions in the light of recent experiments[END_REF][START_REF] Shillcock | The computational route from bilayer membranes to vesicle fusion[END_REF]).

In the past decades, supported lipid membranes have received rapidly increasing attention [START_REF] Kasemo | Biological surface science[END_REF][START_REF] Parikh | Membrane-substrate interface: Phospholipid bilayers at chemically and topographically structured surfaces[END_REF], and are today established model systems to mimic various aspects of cell membranes both for diagnostics (sensors, biochips [START_REF] Boireau | Unique supramolecular assembly of a redox protein with nucleic acids onto hybrid bilayer: towards a dynamic DNA chip[END_REF]) and, e.g., as interesting surfaces for (stem) cell cultures [START_REF] Thid | Supported phospholipid bilayers as a platform for neural progenitor cell culture[END_REF] and also as protein-resistant surfaces [START_REF] Gulcev | Factors affecting the behavior and effectiveness of phospholipid bilayer coatings for capillary electrophoretic separations of basic proteins[END_REF]. This situation motivates to look into the protein adsorption on lipid bilayers in more detail and especially on the adsorption and desorption kinetics. Such adsorption is expected to be negligible if the lipids are neutral, because in this case the protein-

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

lipid interaction is weak (for peptids or, more specifically, for antimicrobial peptides, the situation is often different [START_REF] Bechinger | Rationalizing the membrane interactions of cationic amphipathic antimicrobial peptides by their molecular shape[END_REF], because their size is relatively small and they can easily be incorporated into the bilayers near the hydrophilic lipid head groups).

In vivo, i.e., in real cells, the bilayers contain however both neutral and charged lipids, and such mixtures are also of interest for basic studies. The fraction of the charged lipids is usually relatively low but not negligible. As model systems, charged lipid bilayers may be formed on solid surfaces by adsorption and rupture of vesicles composed of mixtures of zwitterionic, negatively charged, and/or positively charged lipids [START_REF] Richter | Pathways of lipid vesicle deposition on solid surfaces: A combined QCM-D and AFM study[END_REF][START_REF] Richter | Following the formation of supported lipid bilayers on Mica: A study combining AFM, QCM-D, and ellipsometry[END_REF]. All proteins are to some extent charged as well except at the isoelectric point but even there a protein is likely to have local charges. If the charges of the bilayer and protein are opposite, protein adsorption on a bilayer may therefore be favourable even in the absence of specific binding sites. As already noted, this type of adsorption is peculiar, because the lipids forming a bilayer are very mobile [START_REF] Flenner | Subdiffusion and lateral diffusion coefficient of lipid atoms and molecules in phospholipid bilayers[END_REF] and can easily change their spatial distribution in a bilayer in order to increase the interaction with a protein.

The amount of protein adsorbed is determined by the balance of the rates of adsorption and desorption. The activation energy for protein adsorption is usually relatively low. To be adsorbed, a protein should however come in contact with charged lipids. For this reason, the protein adsorption rate should strongly depend on the concentration of charged lipids (globally, it can be limited by diffusion). The activation energy for protein desorption should be appreciable in order to observe measurable amount of adsorption. In the case under consideration, the latter activation energy is determined primarily by the Coulomb interaction between charged amino-acid residues of a protein, charged lipids, and electrolyte. Using the conventional Poisson-Boltzmann theory (the relevant studies are numerous [START_REF] Evans | The Colloidal Domain, where Physics, Chemistry, and Biology Meet[END_REF][START_REF] Asthagiri | Influence of structural details in modeling electrostatically driven protein adsorption[END_REF][START_REF] Koehl | Electrostatics calculations: latest methodological advances[END_REF][START_REF] Mulgrew-Nesbitt | The role ofelectrostatics in proteinmembrane interactions[END_REF][START_REF] Hinderliter | Cooperative adsorption of proteins onto lipid membranes[END_REF][START_REF] Wagner | Electrostatic interactions across a charged lipid bilayer[END_REF][START_REF] Mengistu | Debye-Hückel theory of mixed charged-zwitterionic lipid layers[END_REF]), one can estimate that the contribution of electrostatic interactions to the protein adsorption energy is about 4-6 kcal/mol per charged aminoacid-lipid pair ("charged" means here that the charges of an amino-acid residue and

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

lipid are opposite while globally a pair is neutral), and accordingly a few (4-5) pairs are sufficient for nearly irreversible adsorption at ambient or body temperature. The statistics of the formation of these binding pairs is expected bo be fairly sensitive to the fraction of charged lipids. Our goal here is to illustrate how the rate constants of protein adsorption and desorption can depend on this fraction.

Concerning the subject under consideration, we may notice that various aspects of adsorption of flexible unfolded charged macromolecules on mixed lipid membranes were earlier analyzed in Refs. [START_REF] Dias | Polyion Adsorption onto catanionic surfaces. A Monte Carlo study[END_REF][START_REF] Tzlil | Flexible charged macromolecules on mixed fluid lipid membranes: Theory and Monte Carlo simulations[END_REF] (for the experiment, see e.g. Ref. [START_REF] Bechinger | Rationalizing the membrane interactions of cationic amphipathic antimicrobial peptides by their molecular shape[END_REF]). The applicability of the results obtained there to proteins is however limited, because, despite partial denaturation, adsorbed proteins usually remain in a folded state.

A few available treatments [START_REF] May | Lipid Demixing and Protein-Protein Interactions in the Adsorption of Charged Proteins on Mixed Membranes[END_REF][START_REF] Harries | Adsorption of charged macromolecules on mixed fluid membranes[END_REF][START_REF] Dias | Colloid adsorption onto responsive membranes[END_REF][START_REF] Trusova | Electrostatically-controlled protein adsorption onto lipid bilayer: Modeling adsorbate aggregation behavior[END_REF] of adsorption of charged proteins on charged lipid bilayers (with emphasis on the electrostatic interaction [START_REF] Harries | Adsorption of charged macromolecules on mixed fluid membranes[END_REF][START_REF] Dias | Colloid adsorption onto responsive membranes[END_REF],

adsorption isotherms [START_REF] Harries | Adsorption of charged macromolecules on mixed fluid membranes[END_REF][START_REF] Trusova | Electrostatically-controlled protein adsorption onto lipid bilayer: Modeling adsorbate aggregation behavior[END_REF] and aggregation [START_REF] Trusova | Electrostatically-controlled protein adsorption onto lipid bilayer: Modeling adsorbate aggregation behavior[END_REF]) implicitly assume that a protein is spherical and does not change its shape during adsorption.

In reality, proteins are usually not fully spherical. In addition, proteins are fairly flexible, have a multiplicity of almost isoenergetic states, and can easily change their shape after adsorption. The corresponding conformational changes may range from local reconfiguration near the protein-surface contact to complete denaturation (see, e.g., Refs. [START_REF] Zhdanov | Monte Carlo simulations of the kinetics of protein adsorption[END_REF][START_REF] Castells | Surface-induced conformational changes in lattice model proteins by Monte Carlo simulation[END_REF][START_REF] Raffaini | Protein adsorption on the hydrophilic surface of a glassy polymer: a computer simulation study[END_REF][START_REF] Skepö | Model simulations of the adsorption of statherin to solid surfaces: Effects of surface charge and hydrophobicity[END_REF][START_REF] Kubiak | Molecular dynamics simulations of hen egg white lysozyme adsorption at a charged solid surface[END_REF]). For these reasons, a protein is expected to form a flat contact with a lipid bilayer. This case was treated in Refs. [START_REF] Heimburg | Binding of peripheral proteins to mixed lipid membranes: Effect of lipid demixing upon binding[END_REF][START_REF] Mbamala | Domain Formation Induced by the Adsorption of Charged Proteins on Mixed Lipid Membranes[END_REF][START_REF] Loew | Stability of protein-decorated mixed lipid membranes: The interplay of lipid-lipid, lipid-protein, and protein-protein interactions[END_REF] with emphasis on the adsorption isotherms.

Our analysis below also implies a flat contact between a protein and bilayer.

The main novel ingredient of our work compared to Refs. [START_REF] Heimburg | Binding of peripheral proteins to mixed lipid membranes: Effect of lipid demixing upon binding[END_REF][START_REF] Mbamala | Domain Formation Induced by the Adsorption of Charged Proteins on Mixed Lipid Membranes[END_REF][START_REF] Loew | Stability of protein-decorated mixed lipid membranes: The interplay of lipid-lipid, lipid-protein, and protein-protein interactions[END_REF]] is that we focus on the protein adsorption and desorption kinetics or, more specifically, on the dependence of the rate constants of protein adsorption and desorption on the fraction of charged lipids. There are also more specific differences. In Refs. [START_REF] Heimburg | Binding of peripheral proteins to mixed lipid membranes: Effect of lipid demixing upon binding[END_REF][START_REF] Mbamala | Domain Formation Induced by the Adsorption of Charged Proteins on Mixed Lipid Membranes[END_REF][START_REF] Loew | Stability of protein-decorated mixed lipid membranes: The interplay of lipid-lipid, lipid-protein, and protein-protein interactions[END_REF], for example, the charged amino-acid-lipid pairs forming a protein-bilayer contact are considered to be independent so that the contact energy is linearly proportional to

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

the number of pairs. In our model, the number of protein charges in the contact area is considered to be constant, and the total binding energy is assumed to be maximum when the screening of these charges by charged lipids located in that area is perfect. With decreasing the number of charged lipids there (i.e., with decreasing the number of charged amino-acid-lipid pairs), the total binding energy is considered to decrease rapidly (due to the lack of screening and the corresponding Coulomb repulsion) so that the adsorption becomes energetically unfavourable before breaking the last bond.

Model

In our model (Fig. 1), the charged amino-acid residues are assumed to form one or a few spots near or on the surface of a protein (the location of charges near the surface is usually energetically favourable). After (in reality during) adsorption, a protein slightly changes its shape so that one of the spots, containing m charged amino-acid residues, forms a flat contact with n oppositely charged lipid heads on the surface of the lipid bilayer. The corresponding protein-lipid-bilayer interaction energy is designated to be E n . Strictly speaking, E n is the free energy difference, because it contains entropic contributions. By definition, E n < 0 for the bound states.

To specify E n , we notice that the protein and lipid-bilayer dielectric constants are low, the protein-lipid-bilayer interface is expected to narrow (a few Å), and the corresponding dielectric constant is low as well. In this case,

|E

ACCEPTED MANUSCRIPT

sight, one could postulate that l = 1. In reality, however, as already noted in the Introduction, the protein energy is expected (due to the lack of screening) to become higher than that in the bulk (i.e., E n > 0) at n > 1, and accordingly we should have l > 1. Thus, the desorption process can be viewed as a set of reversible steps of association of an adsorbed protein, P a , with charged lipids, L,

P a L n + L ⇀ ↽ P a L n+1 , (1)
combined with the desorption steps,

P a L n → P u + nL, (2)
where P u is an unbound protein, and n is in the range from l to m -1 in (1) and from l to m in (2).

According to the scheme above, the apparent rate constant of protein desorption is represented as

k d = n=m n=l k n p n , (3)
where k n is the rate constant of step (2), and p n is the probability that a protein has n bonds with charged lipids.

The probability p n depends on the details of the charge distribution on the protein surface and may depend on the interplay of steps ( 1) and ( 2). In our model, in analogy with the transition-state theory, we assume, that steps (1) are rapid [compared to steps (2)] and close to equilibrium, and use the simplest mean-field approximation in order to calculate p n . Specifically, we consider that for a given n the protein-lipid interaction energy, E n , is the same irrespective of the mutual arrangement of the lipids, and the number of the corresponding configurations of

the charged lipids is C m n (C m n = m!/[n!(m -n)!]
is the binomial coefficient). In this case, the grand canonical distribution yields

p n = C m n exp nµ -E n k B T m n=l C m n exp nµ -E n k B T , (4)

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

where µ is the chemical potential of charged lipids, defined under the condition that in the absence of the protein-lipid interaction the energy of a single charged lipid is zero.

The rate constant k n depends on E n and also on the protein-lipid interaction in the activated state for desorption. The latter interaction contributes to the activation energy for protein adsorption. As already noted, the activation energy for protein adsorption is usually relatively low, and accordingly we neglect the proteinlipid interaction in the activated state. In this approximation, we have

k n = κ d exp(E n /k B T ), (5)
where κ d is a constant independent of n.

Substituting ( 4) and ( 5) into (3), we obtain

k d = κ d m n=l C m n exp nµ k B T m n=l C m n exp nµ -E n k B T . ( 6
)
With decreasing n, |E n | rapidly decreases, and if the fraction of charged lipids is not too low, we can take only two terms (with n = m and m -1) in the denominator of expression (6), i.e.,

m n=l C m n exp nµ -E n k B T ≃ exp mµ -E m k B T 1 + m exp - µ + ∆E k B T ,
where ∆E = E m-1 -E m . In this case, expression ( 6) is represented as

k d ≃ k m m n=l C m n exp (n -m)µ k B T 1 + m exp - µ + ∆E k B T . ( 7
)
To describe the process of protein adsorption, we consider that it happens provided that a protein contacts m lipids and at least l of them are charged. Specifically, the protein adsorption rate constant, k a , is assumed to depend only on the corresponding probabilities to contact at least l charged lipids, because as already noted we neglect the protein-lipid interaction in the activated state. In this case, we have

k a = κ a m n=l C m n exp nµ k B T m n=0 C m n exp nµ k B T , (8)

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

where κ a is a constant independent of n. Taking into account that

m n=0 C m n exp nµ k B T = exp µ k B T + 1 m ,
we rewrite expression [START_REF] Shillcock | The computational route from bilayer membranes to vesicle fusion[END_REF] as

k a = κ a m n=l C m n exp nµ k B T exp µ k B T + 1 m . (9)
Eq. ( 9) describes the effect of the protein-lipid interaction on the rate constant of protein adsorption. If the bilayer is partly covered by protein, this rate constant depends also on the hardcore interaction between proteins. The latter effect can be included into the rate constant κ a in analogy with numerous earlier studies of protein adsorption on solid surfaces (see reviews [START_REF]Proteins at Interfaces II[END_REF][START_REF] Hlady | Protein adsorption on solid surfaces[END_REF][START_REF]Biopolymers at Interfaces[END_REF][START_REF] Gray | The interaction of proteins with solid surfaces[END_REF][START_REF] Tsapikouni | Protein-material interactions: From micro-tonano scale[END_REF] and, e.g., Refs. [START_REF] Chatelier | Adsorption of Globular Proteins on Locally Planar Surfaces: Models for the Effect of Excluded Surface Area and Aggregation of Adsorbed Protein on Adsorption Equilibria[END_REF][START_REF] Minton | Effects of excluded surface area and adsorbate clustering on surface adsorption of proteins. II. Kinetic models[END_REF][START_REF] Zhdanov | Kinetics of irreversible adsorption of deformable proteins[END_REF][START_REF] Fang | Kinetics and thermodynamics of protein adsorption: A generalized molecular theoretical approach[END_REF][START_REF] Almeidaa | A simple theory of peptide interactions on a membrane surface: Excluded volume and entropic order[END_REF]). This conventional aspect is beyond the scope of our present study and we do not specify κ a in more detail.

To use the equations presented, we need the relation between µ and θ. In general, this relation is complex. In our examples below, we consider that θ is relatively low (this is the case in vivo and expected often to be the case in vitro), neglect lateral lipid-lipid interaction, and employ the standard mean-field expression,

µ = k B T ln[θ/(1 -θ)], (10)
derived in analogy with that corresponding to the ideal lattice gas.

In our equations above, θ represents the fraction of charged lipids in the area not covered by protein. If the amount of adsorbed protein is low, θ is the same as in the absence of protein, and the equations can be directly used to describe experiments.

With increasing protein coverage, the proteins may trap appreciable part of charged lipids, θ may be lower than that in the absence of protein, and accordingly the protein coverage, θ and µ should be calculated consistently. In particular, θ and µ remain to be related as [cf. Eq. ( 10)]

θ = exp(µ/k B T )/[1 + exp(µ/k B T )]. (11)

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

If for example we consider that the adsorbed proteins are primarily in the states with n = m and m -1 [cf. Eq. ( 7)], the average number of charged lipids per protein is given by n = mp m + (m -1)p m-1 or, more explicitly [with the use of Eq. ( 4)], by

n = m + (m -1) exp[-(µ + ∆E)/k B T ] 1 + m exp[-(µ + ∆E)/k B T ] . (12)
The balance of charged lipids is described as

θ(1 -f ) a + n f b = θ 0 a , (13)
where θ 0 is the fraction of charged lipids in the absence of protein, f is the fraction of the bilayer covered by protein, a is the area per lipid, and b is the protein-bilayer contact area.

Eqs. ( 6), ( 7) and ( 9) alone or in combination with Eqs. ( 11)-( 13) allow one to calculate the protein desorption and adsorption rate constants under various conditions.

Results of calculations

To illustrate the dependence of the protein desorption rate constant on θ (at T = 300 K), we use expression [START_REF] Flenner | Subdiffusion and lateral diffusion coefficient of lipid atoms and molecules in phospholipid bilayers[END_REF]. It shows that the normalized rate constant, k d /k m , depends on three parameters, m, l, and ∆E. To be specific, we employ ∆E = 6 kcal/mol, m = 5, and l = 2, 3 and 4. With these parameters, the desorption rate constant rapidly increases with decreasing θ (Fig. 2), because the protein state with n = m becomes less probable. This effect is especially strong for l = 2. With increasing l, it becomes weaker. Concerning ∆E, we may notice that for physically reasonable values of the latter parameters the results are fairly insensitive to its variation. For example, the results for ∆E = 4 kcal/mol (not shown) nearly coincide with those presented in Fig. 2.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

According to expression [START_REF] Kasemo | Biological surface science[END_REF], the normalized adsorption rate constant, k a /κ a , depends only on two parameters, m and l. With decreasing θ, this rate constant rapidly decreases (Fig. 3), because the contacts with n ≥ l becomes less probable.

During adsorption, θ may decrease. If θ 0 is high (about 0.5 or higher), Eqs.

(11)-( 13) predict that θ remains to be high with increasing f up to saturation of the surface by protein. This means that the protein desorption rate constant remains to be close to minimum, i.e., k d ≃ k m , while the adsorption rate constant is close to maximum, i.e., k a ≃ κ a . This situation is close to the case of adsorption on a surface with fixed uniformly distributed charges.

If θ 0 is relatively low, Eqs. ( 11)-( 13) predict (see, e.g., Fig. 4 for θ 0 = 0.1) that θ ≃ θ 0 only if f is low. With increasing f , θ decreases nearly linearly down to negligible values at f ≃ bθ 0 /(ma). This means that with increasing f the protein desorption rate constant increases, while the adsorption rate constant drops. At f > bθ 0 /(ma), the decrease and increase of these rate constants become dramatic. For these reasons, we can introduce the critical parameter, f c = bθ 0 /(ma). Kinetically, the state with f > f c can hardly be reached even if at f < f c the adsorption is easily feasible.

Conclusion

We have proposed a model describing the likely dependences of the rate constants of protein adsorption and desorption on the fraction of charged lipids in a lipid bilayer. Our analysis shows that the redistribution of charged lipids in the bilayer during protein adsorption results in interesting features of adsorption and desorption quite different from those on a surface with immobile charges and binding sites. In particular, our model indicates that there is critical fraction of the area covered by protein. Adsorption above this fraction is hindered.

Concerning the limitations of our model, it is appropriate to articulate that we

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

have assumed that although a protein is deformable its adsorption is not accompanied by denaturation. In reality, the protein-lipid electrostatic interaction may of course result in protein denaturation and it may change the kinetics of adsorption and especially desorption.

Finally, we may notice that the protein adsorption on a lipid bilayer represents the first step of the incorporation of a protein into a bilayer. For this reason, our results are of interest from the latter perspective as well. f

  n | is maximum for perfect charge screening, i.e., for n = m. This state is the most probable if the fraction of charged lipids, θ, is not too low. With decreasing θ, the states with n < m may become more probable than that with n = m. With decreasing n, |E n | rapidly decreases and eventually the protein adsorption becomes energetically unfavourable. It happens at n < l, where l (l < m) is the number of bonds in the bound state with the minimum binding energy. At first A C C E P T E D M A N U S C R I P T

Fig. 1 .

 1 Figure captionsFig.1. Scheme of protein adsorption on a bilayer composed of lipids with neutral and charged heads. The charged amino-acid residues form a spot on the surface of a protein. After adsorption, a protein changes its shape so that the charged spot forms a flat contact with the oppositely charged lipid heads of the lipid bilayer.

Fig. 2 .

 2 Fig. 2. Rate constant of protein desorption as a function of the fraction of charged lipids for ∆E = 6 kcal/mol, m = 5, and l = 2, 3 and 4.

Fig. 3 .

 3 Fig. 3. Rate constant of protein adsorption as a function of the fraction of charged lipids for m = 5 and l = 2, 3 and 4.

Fig. 4 .

 4 Fig. 4. Fraction of charged lipids in the area not covered by protein as a function of the fraction of the area covered by protein (according to Eqs. (11)-(13)) for θ 0 = 0.1, m = 5, ∆E = 6 kcal/mol, and b = 10a.

Figures 1- 4 A

 4 Figures1-4

A C C E P T E D

M A N U S C R I P T


Giulia Morra

Chiara Baragli

Giorgio Colombo
email: g.colombo@icrm.cnr.it.

Selecting Sequences that Fold into a Defined 3D Structure: A New Approach for Protein Design Based on Molecular Dynamics and Energetics

Keywords: Protein Folding, Protein Design, Molecular Dynamics, Self-organization

The problem of finding amino acid sequences able to fold into a defined three-dimensional (3D) structure is at the basis of successful protein design efforts.

Herein, we present the results of the application of a novel, all-atom molecular dynamics based, energy decomposition approach to the selection of sequences able to fold into a given 3D conformation. First, the energy decomposition approach is applied to natural sequences associated to a well-defined structure to identify the principal energetic coupling interactions necessary to stabilize it, defining the specific energetic signature for the fold. Then, several different sequences are threaded on the defined 3D structure and only those sequences whose energetic signature (pattern) is close to that of the natural sequence, according to a similarity criterion, are selected as able to populate the specific fold. Furthermore, it is possible to evaluate the fitness of a certain sequence for a fold by combining the information provided by the energetic signature to that contained in the contact map, which recapitulates the fold topology. The results show that the better fit between the energetic properties of a sequence and the topology corresponds to a better stabilization of the protein fold by that sequence.We applied this approach to a library of natural and artificial WW domain sequences, previously developed by the Ranganathan group, containing sequences that are experimentally known to be able and unable to fold into native structures. The results show that our approach can correctly identify 70% of the sequences known to populate the typical WW domain fold.

A C C E P T E D M A N U S C R I P T

Introduction.

Successful protein design relies on the correct identification of sequences that fold into defined threedimensional (3D) structures. This problem, also known as the "inverse protein folding" problem, can be tackled effectively by specifying what information in the sequence is necessary and sufficient to determine a certain fold.

Since the seminal work of the Eisenberg group [START_REF] Bowie | A method to identify protein sequences that fold into a known three-dimensional structure[END_REF], considerable progress has been made the development of computational methods for identifying amino acid sequences compatible with a target structure [START_REF] Watters | The highly cooperative folding of small naturally occurring proteins is likely the result of natural selection[END_REF][START_REF] Kuhlman | Design of a novel globular protein fold with atomic-level accuracy[END_REF][START_REF] Lee | Molecular dynamics in the endgame of protein structure prediction[END_REF][START_REF] Dahiyat | De novo protein design: Fully automated sequence selection[END_REF][START_REF] Desjarlais | De novo design of the hydrophobic cores of proteins[END_REF]. Mayo and coworkers reported one of the most notable examples, namely the complete redesign of a zinc finger protein [START_REF] Dahiyat | De novo protein design: Fully automated sequence selection[END_REF]. Using different atomistic energy functions that mimic the physical interactions between aminoacids, several other groups have achieved outstanding successes in redesigning natural folds, in the de novo construction of novel folds, or in the re-design of enzymes [START_REF] Kortemme | Computational redesign of protein-protein interaction specificity[END_REF][START_REF] Kraemer-Pecore | A de novo redesign of the ww domain[END_REF][START_REF] Harbury | High resolution protein design with backbone freedom[END_REF].

The basic principle in atomistic design is the optimization of a target potential function providing sequences with defined thermodynamic minima corresponding to the native configurations, well separated from alternative conformational states. Having a deep free-energy minimum for the native state conformation ensures the production of a sequence with high thermal stability. However, the native states of natural proteins reside in shallow free energy minima corresponding to partially stable and dynamic folds, characterized by a multiplicity of (similar) conformations.

Starting from these concepts, Ranganathan et al. proposed a different strategy based on applying Statistical Coupling Analysis (SCA) to multiple sequence alignments of a protein family, to identify the mutual inter-residue dependencies evidenced by conserved statistical correlations between amino acid distributions at specific sites [START_REF] Lockless | Evolutionarily conserved pathways of energetic connectivity in protein families[END_REF][START_REF] Socolich | Evolutionary information for specifying a protein fold[END_REF][START_REF] Suel | Evolutionary conserved networks of residues mediate allosteric communication in proteins[END_REF][START_REF] Shulman | Structural determinants of allosteric ligand activation in rxr heterodimers[END_REF][START_REF] Russ | Natural-like function in artificial ww domains[END_REF]. The application of this approach showed that a small set of residues at specific positions (in a certain protein family) coevolves among a majority which are largely uncoupled, and that the strongly coevolving residues are organized into spatially connected networks stabilizing their respective structures through packing interactions. If used in the design and selection of new sequences folding to a certain target-structure, this purely statistical, mechanism-free method should in principle produce sequences with the same marginal stabilities and biological functions as those of natural proteins. Ranganathan and coworkers were actually able to design artificial WW domains showing thermodynamic and structural properties in excellent agreement with the ones of their natural counterparts [START_REF] Socolich | Evolutionary information for specifying a protein fold[END_REF][START_REF] Russ | Natural-like function in artificial ww domains[END_REF]. Moreover, the authors showed that the artificial sequences could

A C C E P T E D M A N U S C R I P T

perform the same function as the native ones, showing class-specific recognition of proline-containing target peptides [START_REF] Socolich | Evolutionary information for specifying a protein fold[END_REF][START_REF] Russ | Natural-like function in artificial ww domains[END_REF].

In this paper, we aim to analyze the energetic determinants of the sparse architecture of residueresidue interactions necessary to stabilize a certain fold, based on the analysis of the conformational dynamics and interactions in the native state of a natural protein, and to use this information for the selection of other, non-natural sequences able to fold to the same target 3D structure.

This approach is based on a recently introduced Energy Decomposition Method, aimed at identifying the key residues (interaction hot spots) for the stabilization and folding of the protein to a defined 3D structure [START_REF] Tiana | Understanding the determinants of stability and folding of small globular proteins from their energetics[END_REF][START_REF] Ragona | Determinants of protein stability and folding: Comparative analysis of beta-lactoglobulins and liver basic fatty acid binding protein[END_REF][START_REF] Colacino | The determinants of stability in the human prion protein: Insights into the folding and misfolding from the analysis of the change in the stabilization energy distribution in different condition[END_REF][START_REF] Colacino | Similar folds with different stabilization mechanisms: The cases of prion and doppel proteins[END_REF][START_REF] Morra | Relationship between energy distribution and fold stability: Insights from molecular dynamics simulations of native and mutant proteins[END_REF]. Previously, we showed the ability of this method to capture the essential changes occurring in the energetics of a protein upon single amino acid mutation [START_REF] Tiana | Understanding the determinants of stability and folding of small globular proteins from their energetics[END_REF][START_REF] Ragona | Determinants of protein stability and folding: Comparative analysis of beta-lactoglobulins and liver basic fatty acid binding protein[END_REF][START_REF] Colacino | The determinants of stability in the human prion protein: Insights into the folding and misfolding from the analysis of the change in the stabilization energy distribution in different condition[END_REF][START_REF] Colacino | Similar folds with different stabilization mechanisms: The cases of prion and doppel proteins[END_REF][START_REF] Morra | Relationship between energy distribution and fold stability: Insights from molecular dynamics simulations of native and mutant proteins[END_REF]. These changes are mainly related to stability variations in an ensemble of single-point mutants of a certain protein [START_REF] Morra | Relationship between energy distribution and fold stability: Insights from molecular dynamics simulations of native and mutant proteins[END_REF].

The main obstacle in trying to define what properties of a sequence are necessary to define a certain fold and how structural constraints impact on the selection of a certain sequence, using atomic level resolution for the study of interactions, is represented by the vast complexity of the energetic interactions between amino acids. The Energy Decomposition Method alleviates this problem by providing a simplified view of stabilizing interactions, extracting the major contributions to energetic stability of the native structure from all-atom molecular dynamics (MD) simulations. In this method, for a protein of N residues, the matrix of average non-bonded interactions between pairs of residues is built from an MD trajectory. The energy map is then simplified through eigenvalue decomposition (Principal Component Analysis) [START_REF] Tiana | Understanding the determinants of stability and folding of small globular proteins from their energetics[END_REF][START_REF] Ragona | Determinants of protein stability and folding: Comparative analysis of beta-lactoglobulins and liver basic fatty acid binding protein[END_REF][START_REF] Colacino | The determinants of stability in the human prion protein: Insights into the folding and misfolding from the analysis of the change in the stabilization energy distribution in different condition[END_REF][START_REF] Colacino | Similar folds with different stabilization mechanisms: The cases of prion and doppel proteins[END_REF][START_REF] Morra | Relationship between energy distribution and fold stability: Insights from molecular dynamics simulations of native and mutant proteins[END_REF].

The eigenvector associated with the lowest eigenvalue is made of N components, each one describing the contribution to stabilization energy provided by the corresponding protein residue. Each of the components describes the contribution of the respective aminoacid to the stabilization energy of the protein. Analysis of the N components of the eigenvector associated with the lowest eigenvalue was shown to single out those residues (hot sites) behaving as strongly interacting and possible stabilizing centers. In general, these residues constitute a network of strongly coupled interactions typical for a certain fold. This vectorial representation of the sequence (sequence eigenvector, SE) may be thought of as the "energetic" signature of that fold [START_REF] Morra | Relationship between energy distribution and fold stability: Insights from molecular dynamics simulations of native and mutant proteins[END_REF].

The lowest eigenvalue represents an effective coupling parameter: a variation in the first eigenvalue due to mutations or structural changes can be interpreted as a change (rescaling) in the strength (intensity) of all stabilizing interactions introduced by the mutation. A more detailed and quantitative description of the method is given in Materials and Methods.

A C C E P T E D M A N U S C R I P T

A similar reasoning could be applied to the analysis of the structural properties. The native state structure, or designed target geometry, can be described in terms of the matrix of its native contacts (the contact matrix). This provides the essential geometrical definition of the topology of the native structure. It is known that the native state topology is a major determinant of the folding free-energy landscape of many (small) proteins. The vectorial representation of the topology of the native state is defined by the principal eigenvector of the native contact matrix (contact eigenvector, CE), which depends on the desired 3D structure [START_REF] Morra | Relationship between energy distribution and fold stability: Insights from molecular dynamics simulations of native and mutant proteins[END_REF].

The validity of the vectorial representation of stabilization energy was previously checked in the context of the calculation of the relative stability of single mutants of several proteins [START_REF] Morra | Relationship between energy distribution and fold stability: Insights from molecular dynamics simulations of native and mutant proteins[END_REF], showing good correlations between theoretical and experimental data. The components of the first eigenvector define the main attractive couplings that stabilize a certain folded state. In related protein mutants that can still fold properly to the native structure, mutations can either modulate the coupling intensity of these specific interactions (reflected in the value of the eigenvalue), or modulate the height of some peaks in the first eigenvector, without disrupting the overall signature of the profile [START_REF] Morra | Relationship between energy distribution and fold stability: Insights from molecular dynamics simulations of native and mutant proteins[END_REF]. We could also show that the similarity, defined in terms of the Pearson's correlation coefficient, between the sequence eigenvector (SE) of a certain sequence and the contact eigenvector (CE) of the native structure correlated reasonably well with the relative stability of the corresponding protein [START_REF] Morra | Relationship between energy distribution and fold stability: Insights from molecular dynamics simulations of native and mutant proteins[END_REF].

Building on these considerations, we set out to test the possibility of this approach to discriminate, in a large ensemble of sequences, those that are able to fold to a desired structure vs. those that are not. We selected a subset of the same natural and artificial WW-domain sequences, tested by Ranganathan in his seminal paper on SCA, as a means for sequence selection [START_REF] Socolich | Evolutionary information for specifying a protein fold[END_REF]. The conformational dynamics and energetics of each sequence were probed by all-atom Molecular Dynamics (MD) simulations in explicit water at 300K. The discrimination between folding vs. non-folding sequences was based on the calculation of the similarity between the SE for each simulated sequence and the equivalent vector in the natural WW domain protein, used for reference and not included in the set of Native sequences.

Moreover, based on previous results demonstrating that higher correlations between SE and CE successfully identify proteins endowed with higher structural stability, we investigated the Pearson's coefficient similarity between SE's and CE's as a measure for discriminating productive folders from non-folders, independently of the previous knowledge of the SE of the native protein. This aspect may be relevant in the design or modification of novel folds where the main information available may be the geometry of the desired target, in the absence of statistically relevant information on the SE's of known structural homologues. Ideally, the knowledge of the topology of the target structure can be No specific new free-energy function was built for the selection. Analysis of the results showed that the combination of energy decomposition and topological analysis is able to correctly identify 70% of the sequences folding to the natural WW structure.

A C C E P T E D M A N U S C

Results

The analysis of the foldability of different sequences presented here is based on a description of the complex non-bonded energy of a protein through the approximated stabilization energy E nb app [START_REF] Morra | Relationship between energy distribution and fold stability: Insights from molecular dynamics simulations of native and mutant proteins[END_REF].

The simplification is achieved by means of the first (most negative) eigenvalue and first eigenvector obtained by decomposition and diagonalization of the energy matrix from all-atom MD simulations (see Methods).

The eigenvector describes for each residue the amount of energy coupling it shares with all other residues in the native state of the protein. Upon mapping these couplings on the 3D structure of the protein, a connected network of strong interactions is revealed, involving distant residues in the sequence. These residues correspond to the most intense peaks. Most importantly, the first eigenvector reports on an organization of the energetics of the native state that is typical for a certain fold, defining an energetic signature for that fold. The information contained in the main energetic eigenvector (Sequence Eigenvector, SE) of a related set of sequences known to fold to the desired 3D structure can be used as a template to search for other sequences able to stabilize the same fold in native conditions, without limiting to single or double mutations.We applied the Energy decomposition method to the native YAP65 WW domain sequence whose X-ray structure is available in the Protein Data Bank (PDB id: 1k9r), in order to determine the energetic signature of the protein. Then, we considered four

A C C E P T E D M A N U S C R I P T
sets of 8 sequences each, randomly extracted from the natural and artificial, folding and non folding, sequences from the four groups present in the study of Socolich et al . [START_REF] Socolich | Evolutionary information for specifying a protein fold[END_REF]. These are:

-Native sequences (N): they occur naturally and do fold into the WW domain structure.

-CC (coupled conservation) sequences: artificial sequences, created on the basis of the SCA and on the premise that conservation of the pattern of coupled interactions seen in natural sequences is sufficient to favor the folding to the desired structure [START_REF] Socolich | Evolutionary information for specifying a protein fold[END_REF]. Experimentally, a significant percentage of them proved to fold to the native structure by circular dichroism (CD) and NMR analysis. Our selected subset of CC sequences comprises only sequences that fold to a native WW structure under the considered experimental conditions in [START_REF] Socolich | Evolutionary information for specifying a protein fold[END_REF][START_REF] Russ | Natural-like function in artificial ww domains[END_REF].

-IC: artificial sequences, created from multiple sequence alignments based on the hypothesis that conservation is a property of one single site, independently of others. Experimentally, no folding to the native WW domain geometry was observed [START_REF] Socolich | Evolutionary information for specifying a protein fold[END_REF].

-Random sequences (R): they were built by randomly mutating native sequences. As expected, they do not fold [START_REF] Socolich | Evolutionary information for specifying a protein fold[END_REF].

As shown by Socolich et al. [START_REF] Socolich | Evolutionary information for specifying a protein fold[END_REF] conservation of amino acid composition as inferred by a multiple sequence alignment is not sufficient to discriminate productive folders (CC) from non-folding sequences (IC).

All sequences and their labels are reported in the Supplementary Material

Correlation between sequence, energetics and folding to a given 3D structure.

The energy decomposition method was first applied to the wild type sequence of the YAP65 WW domain (PDB id: 1k9r). The resulting SE profile recapitulates the information on which residues at which positions are most important in the stabilization of the 3D native structure. Peaks in the SE correspond to those residues whose pair interactions with the rest of the amino acids contribute to the protein stability (Figure 1). In particular, pair interactions between two peak-related residues provide a significant stabilizing energy to the protein. As a natural consequence, the SE reports therefore also on the residue-residue couplings defining the network of interactions that contribute to the folding core of the protein, whose participating residues are indicated by the SE regions above a defined threshold [START_REF] Tiana | Understanding the determinants of stability and folding of small globular proteins from their energetics[END_REF].The SE of the YAP65 WW domain (Figure 2) indicates as energetically relevant residues a subset comprising the segments E8-S13 and Q17-D25 (numbering as in Socolich et al. [START_REF] Socolich | Evolutionary information for specifying a protein fold[END_REF]). Relevant residues correspond to positions in the principal eigenvector characterized by a high value of the

A C C E P T E D M A N U S C R I P T component.
Graphically, most relevant residues (and the values of their components) are identified by peaks in the eigenvector profile. These segments are known to contain strongly coevolving residues (such as E8, Y21, H23) responsible for ligand binding and function in WW domain. Hence, according to our analysis these residues also provide a stabilizing contribution to the protein.Following this first step, the SE's for a set of eight natural WW domain sequences were calculated. The peaks in the profile of the SE's, indicating residues participating in stabilizing pair interactions, are rather conserved and the distribution between values over and under the threshold reflects the wild type case.The Pearson's correlation coefficient was used to define the similarity between the SE of 1k9r with each of the SE's calculated for the other native sequences (N set, Figure 2)

The results showed generally high values, in the range between 0.7 and 1 for the Pearson's correlation (see Table 1), and an average value of 0.80, indicating that the specific pattern of interactions defined by the Sequence Eigenvector (SE) profile may actually be considered important in the discovery of the energetic determinants of the folding features of the protein. Next, attention was focused on the artificial sequences proposed and tested by Ranganathan [START_REF] Socolich | Evolutionary information for specifying a protein fold[END_REF].

First, the CC group was analyzed, following the same procedure as for native sequences: for each of the eight CC sequences SE was calculated and its similarity to SE of 1k9r was measured by Pearson's coefficient. For this set, the SE profiles still retain the modulation of the wild type (Figure 3), however they show some increased variability in the peak intensities. Pearson's coefficients' values are lower than for the Native set, but still in general good agreement with the values computed for native sequences (between 0,6 and 1; Table 1, with an average of 0.72). Finally the same analysis was performed on IC and Random sequences. The IC set comprises sequences that are not folding to the native state in spite of the native-like amino acid composition site by site, resulting from simple sequence alignment analysis. The SE profiles are as noisy as in the CC case, hence these sequences do not seem to be distinguishable from the folding ones. However, when looking at the Pearson's coefficients, they turn out to be consistently lower, with an average of 0.65. For random sequences, which are expected to produce totally uncorrelated SE's, Pearson's coefficient values of for IC and Random sequences are very low values (Table 1) with an average of 0.45.

These results show that the SE is actually capable of capturing the significant energetic features for the fold. In general, higher correlations are found between the energy couplings of the native sequence and those of proteins known to fold to the typical WW domain structure. In contrast, random sequences and also the sequences missing important coupling interactions in spite of high sequence similarities (the IC set) are characterized by a much lower correlation suggesting that the energetic signature captured

A C C E P T E D M A N U S C R I P T
by the SE is able to discriminate between the folding propensities of different sequences. The method based on the similarity with the reference WT protein correctly identifies native sequences and to a significant extent also the sequences belonging to the CC set, thus indicating that information on paircorrelations is included in the SE. The different behavior of CC and IC sequences suggests that the foldability of a sequence appears to be encoded in more subtle sequence properties than the residue composition and single site distribution. In order to check whether the differences in the correlations between SE's of CC and IC and the SE of the Wild Type (WT) does not trivially result from a sequence similarity, we re-evaluated specifically for each sequence in our dataset the alignment score with the WT sequence [START_REF] Larkin | Clustalw and clustalx version 2[END_REF] (Table 1). While Native sequences generally show high alignment scores with the WT, the CC and IC sequences have comparable similarities and cannot be distinguished from one another based on sequence features only. Moreover, the correlation between the results of the two methods was evaluated: the Pearson's coefficients between the SE's of IC and CC peptides with WT-SE were plotted against the respective sequence alignment scores, yielding a linear correlation coefficient of 0.17, excluding random sequences from the calculation (two-sided p-value=0.406).

Interestingly, the discrimination between folders and non-folders cannot be obtained by looking at standard global properties of the structure like RMSD (data not shown). Within 5ns of MD trajectory, no significant unfolding or global structural rearrangements are possible even for a small protein like the WW domain. Interestingly, neither the total interaction energy, that is provided for each structure by the force field parameters, nor the approximate stabilization energy calculated by the Energy Decomposition Method (as already proved in [START_REF] Morra | Relationship between energy distribution and fold stability: Insights from molecular dynamics simulations of native and mutant proteins[END_REF]), are able to distinguish between folders and nonfolders.

Still, the distribution of stabilizing interactions as it is described in the SE turns out to be a well defined measure for specific sequence properties such as the ability to select favorable native contacts providing stability and cooperativity to structure formation, hence determining the sequence foldability.

Correlation between energetics and topological properties of the target 3D structure. In order to add the topological information on the fold to our analysis, the contact matrix for the WW domain was calculated, and subsequently subjected to eigenvalue and eigenvector analysis. The eigenvector associated to the highest positive eigenvalue (principal eigenvalue) was considered to be representative of fold topology, and is generally referred to as the Contact Eigenvector (CE). In analogy to SE, CE

A C C E P T E D M A N U S C R I P T

indicates which residues constitute the essential determinants of the domain architecture [START_REF] Bastolla | Principal eigenvector of contact matrices and hydrophobicity profiles in proteins[END_REF]: in general, CE singles out residues that have a high number of contacts with other residues (Figure 4).

In a previous study [START_REF] Morra | Relationship between energy distribution and fold stability: Insights from molecular dynamics simulations of native and mutant proteins[END_REF], we demonstrated that the degree of correlation between the SE and CE is a measure of the fitness of a certain sequence to a certain fold. A suitable sequence for a certain fold places the strongly interacting residues (hot spots) where they can stabilize the structure. An example is the buried core of the protein where several residues must be tightly packed to develop energetic interactions responsible for correct folding. As a consequence, the SE of a folding sequence should be similar to the CE of the fold. Importantly, by applying this concept, it was possible to rank the differences in stability of a diverse set of mutants of a series of proteins with remarkably different folds. The comparison between SE and CE allows to shed light on the degree of compatibility between a specific sequence, which provides characteristic energetic interactions, and the 3D structure that is being evaluated. Here, we attempt to evaluate to what extent the similarity between CE and SE correlates with folding properties of the sequence.

In this context, CE was calculated for the crystal structure of 1k9r as a representative of the WW fold.

Subsequently the Pearson's coefficients between the SE profiles of the native sequences from the pdb and the CE were measured: as expected, the results show great similarity (with an average of 0.87; Table 1). Slightly lower values were obtained for CC sequences (average 0.82). However, also

Pearson's correlations between IC sequences SE's and WW domain CE gave similar results to those of the N sequences (0.88). Finally, the similarities between the SE's and the native CE for random sequences display in general minimal values (0.65).

Taken together, this body of results suggests that the measure of similarity of the SE of a certain sequence to that of the native sequence may effectively discriminate between sequences that are either able or unable to fold to a target structure. Moreover, they suggest that while a high degree of similarity between CE of the target structure and the SE's of different sequences is necessary for a sequence to fold, it is not a sufficient criterion to determine whether a sequence can actually populate that specific fold. According to this result, the large majority of non-random sequences in our dataset satisfy this criterion, hence they are compatible with the native structure in terms of "topological" requirements, such as the number of contacts formed (residue size and hydrophobicity) at each single site. This is reflected in the linear correlation coefficient between SE-CE similarities and sequence similarities, which results equal to 0.50, in contrast to the 0.17 calculate above (two sided p-value 0.009).

A C C E P T E D M A N U S C R I P T

Combining energetic and topological information to select viable sequences for a defined 3D fold.

Finally, the information on the correlations between the SE's of different sequences and that of the natural sequence folding to the WW domain 3-stranded geometry (1k9r) was combined with the information on the correlations between the SE's of different sequences and the CE recapitulating the properties of the WW-domain fold. The two quantities were plotted in a graph (Figure 5). The folding sequences (native and CC) define an ensemble mainly located in the right-upper part of the graph, separated from nonfolding sequences (IC and random). The limit of 0.7 for Pearson's coefficient on both axes, defines an area that contains mostly folding sequences: 11 out of 16 folding sequences (69%) are in this part of the graph, and 13 out of 16 nonfolding sequences (81%) are out of it. Hence, by setting the acceptance threshold to 0.7 we obtain 5 false negatives and only 3 false positives.

Among the false negatives, only one native sequence is not recognized (N7), possibly because of a very high peak (due to the strong relative contribution at position 23) which slightly alters the overall relative distribution of peaks (as an effect of normalization). By introducing a more restrictive threshold, such as 0.8, the number of false positives drops to zero, whereas the number of false negatives increases by two units.

The combined analysis of considering both the energetics of the sequences and the topological features of the fold proves the possibility to discriminate between folding and non-folding sequences, with nearly 70% accuracy and a limited computational cost. In this procedure, we used the SE of 1k9r as a reference, representative for the determinant residue-residue coupling interactions necessary to fold into the WW domain structure. Although the chosen WWdomain was selected randomly among all structures present in the Protein Data Bank, it might introduce some bias in the classification of folding/nonfolding sequences, which is based on the similarity to its SE.

Therefore, in order to be independent of previous knowledge, a final test was performed without any previous assumption on the similarities with a certain reference sequence. In this line of thought, for each sequence, the correlation between its SE and SE's of all other sequences were measured by

Pearson's coefficient. The resulting values were clustered by means of a cluster analysis. This procedure highlights the presence of a distinct, dominant cluster containing sequences that are similar to one another in terms of SE's. Strikingly, most of the sequences in this cluster were proved experimentally to fold to the required WW-domain geometry. In detail: 13 out of 17 folding sequences, including the wild type, (76%) belong to this cluster and 10 out of 16 non-folding sequences (63%) do not. However, the number of false positives also increases to 7, and includes only IC sequences.These

A C C E P T E D M A N U S C R I P T

results confirm that folding sequences show similar energetic features, so they can be clustered together and distinguished from non-folding sequences; moreover, SE captures with good approximation the most significant energetic characteristics of each sequence. Therefore, the energy decomposition method can be a useful tool to investigate whether a sequence is likely to fold on a required structure.

In this context, the CE-SE similarity criterion should be considered only as a necessary, but not sufficient, condition for a sequence to fold into a required geometry.

Discussion and Conclusion.

In this paper we applied a Molecular Dynamics based strategy to predict, given a protein fold and a set of sequences with very similar chemical features (composition, single site distribution etc), which ones are able to fold to the given structure and which ones are not, based on the hypothesis that foldability where site correlations are neglected, is the information used the to build the artificial IC sequences, which hence satisfy the required constraints in terms of local burial as the N and CC sequences do. This observation correlates with the experimental finding pointed out by Ranganathan [START_REF] Socolich | Evolutionary information for specifying a protein fold[END_REF], that IC sequences are generally soluble, which is not the case for the majority of R sequences. This fact could indicate that the hydrophobic collapse resulting in protein solubility is accounted for by the correct hydrophobic-polar pattern along the chain, while the unique shape of the native fold requires a network of specific cooperative interactions, also involving distant sites, somehow not present in the IC sequences and whose absence is not detected when comparing their SE to the protein's CE. This is confirmed by the statistically significant correlation we found between sequence similarity and SE-CE similarity. The SE-CE similarity evaluates the basic agreement of a candidate sequence with a target structure. In light of the results shown here, the SE-CE similarity should be considered only as a necessary, but not sufficient, condition for a sequence to fit to a certain topology. Equivalent information could be obtained by comparing the sequence to a multiple alignment of natural sequences folding to the given structure. In this respect, our approach offers the advantage of requiring in principle only the structure of the natural template protein, or information on the topology of a designtarget structure recapitulated by the CE. In a typical design application the exact native structure is not available. In this context, we speculate that at the present stage the SE-CE correlation of our method might be used as mainly as a preliminary screen for sequences that can be efficiently threaded on the structure.

While the CE-SE comparison can identify sequences having a single site distribution that fit the requirements of the structure, the specific energy-related details of the cooperative interactions seem to be captured by the comparison of each sequence's SE with the template WT-SE, which in fact proves able to discriminate between CC and IC sequences. Moreover, the set of native sequences obtains the highest average similarity score with the template. The natural sequences (N) are likely to have evolved a more complex cooperativity not entirely described by the correlation data used to design the CC set, which nevertheless were shown to be essential and also sufficient for the domain to fold.

SE reflects the energetic properties of the native state. In a previous paper, we could show that the approximate stabilization energy recapitulated by the SE can be used as an effective approximation of the enthalpic part of the folding free energy. This was applied for estimating the stability difference between two sequences differing for only one mutation on a total of about 60 residues. In that case, in light of the small perturbation, it was reasonable to assume that the energetic properties of the unfolded

A C C E P T E D M A N U S C R I P T

state are essentially identical for the two sequences, allowing a good approximation of the folding enthalpy [START_REF] Morra | Relationship between energy distribution and fold stability: Insights from molecular dynamics simulations of native and mutant proteins[END_REF]. This hypothesis may not hold here, considering sequences differing at many positions.

The number of perturbations may actually alter the distribution of conformations and interactions defining the unfolded states for each of the sequences.

Given the impossibility to fully characterize the unfolded states for many different sequences, a more accurate calculation of the free energy differences between sequences could be obtained with free energy perturbation methods and thermodynamic cycles connecting different mutants. However, free energy perturbations require long equilibration times at intermediate values of the Hamiltonian coupling parameters and result in time-consuming efforts. Moreover, limitations due to the accuracy of the force field parameters and sampling issue may impact on the final outcome of the calculations.

These are out of the scope of this paper.

The present analysis and data set is based on a randomly and arbitrarily chosen template, namely the native YAP65 WW domain 1k9r, whose specific properties might introduce some bias in the sequence classification. With the aim of removing this bias, we attempted the clustering strategy outlined at the end of the Results section, considering the similarity between all sequence pairs and not only with respect to the template sequence and then clustering the data into two sets, in order to distinguish between folders and non folders. The performance of the classification method improves slightly, but the template structure still might have some influence on the results. In the future, one may think of optimizing the selection of the template structure, for instance by means of a clustering procedure on the PDB data entries. Such a procedure might be suitable for inverse protein folding applications, where the target structure is known and new sequences can be selected or optimized by means of the present classification scheme. The acceptance threshold of 0.80 set for the SE, WT-SE correlation proves able to discriminate between folders and non-folders in the case of WW domain, a small protein for which very similar sequences are analyzed. In spite of the significant similarities among the sequences in the data set, both the less and the more restrictive thresholds of 0.70 and 0.80 chosen, allow the discrimination of productive folders from non-folding sequences based only on a mathematical and numerical descriptor. This constitutes a valid prerequisite to the general applicability to different systems.

The method proves suitable for applications in the field of protein design, limited at the moment to the pre-selection of sequences that may fit to a certain fold, and consequently as a sequence classification tool. In this initial work, we made use of information regarding the template structure, while in a fully

A C C E P T E D M A N U S C R I P T

ab initio protein design application, the exact native structure of the protein may not be available and this may limit the performances of our method. It is important to recognize that further theoretical developments coupled with experimental characterization of designed sequences are required in this context. From the theoretical point of view, for instance, one may benefit from efforts to characterize the unfolded states that would give a better description of the free energy components involved in the stabilization of native states. Alternatively, one could use simulation approaches allowing wider sampling to properly weight the statistical relevance of different conformations. All these efforts require parallel advances in the performances of MD algorithms and in the development of hardware.

Recent reports on these subjects hold great promise with regards to the possibility to apply the methods described here on a much larger scale and with much better sampling. [START_REF] Klepeis | Long-timescale molecular dynamics simulations of protein structure and function[END_REF][START_REF] Shaw | Anton, a specialpurpose machine for molecular dynamics simulation[END_REF][START_REF] Van Der Spoel | Gromacs user manual version 3.2. Www.Gromacs.Org[END_REF][START_REF] Hess | Gromacs 4: Algorithms for highly efficient, load-balanced, and scalable molecular simulation[END_REF] Finally, the applications presented here require relatively short MD simulations. The choice of 5-ns long MD simulations was actually an arbitrary one, which we considered a good compromise between computational efficiency and sufficient relaxation of the structures, allowing the simulation of 32 different sequences at the all-atom level of resolution. In our experience, the Molecular Dynamics simulations spanning a short time range like 5ns, while not ensuring the complete convergence to the equilibrium structure, represent a viable compromise to reach a sufficient relaxation of the side chain and local backbone conformations, yielding a well-defined pattern of non-bonded pair interactions, with a limited computational effort. Therefore, the limited size of the required Molecular Dynamics trajectories (5 ns) as well as the speed up offered by the use of implicit solvent models in the calculation makes this approach compatible with a large scale application, such as multiple sequences screening, at least in the case of small to medium size molecules.

Materials and methods.

Socolich et al. [START_REF] Socolich | Evolutionary information for specifying a protein fold[END_REF] studied and experimentally tested the fold of many alternative sequences of the WW domain, using 1H-NMR tests and thermal denaturation. We have applied our method to a subset of these sequences in order to test its capability in discriminating the folding sequences. The sequences can be divided into 4 groups, they are labeled as in Socolich et al (N for Native sequences, CC for Coupled conservation, IC for Independent conservation and R for Random) [START_REF] Socolich | Evolutionary information for specifying a protein fold[END_REF], and reported in Table 1 and Supplementary Material.

We selected 8 sequences within each group, with the only requirements that the CC seqeuences would fold to a native WW domain under the experimental conditions, and that they did not contain gaps with respect to the template sequences, to avoid discrepancies when mounting the sequence onto the template structure.

These sequences were associated to the structure of a YAP65 WW domain, PDB id: 1k9r (NMR structure), extracted from the protein data bank. The template sequence corresponds to residues 8-40 of the original protein, in the following numbered from 1 to 33 in agreement with the numbering in Socolich et al [START_REF] Socolich | Evolutionary information for specifying a protein fold[END_REF].

We used Reconstruction program described in [START_REF] De Mori | Study of the villin headpiece folding dynamics by combining coarse-grained monte carlo evolution and all-atom molecular dynamics[END_REF] to thread the sequences on 1k9r fold: this program conserves only the alpha-carbon trace in 1k9r structure and adds the side chains of the residues of the sequence. Then, it determines the coordinates of the side chain atoms, using a library of rotamers. This library contains empirical data, extracted from 100 NMR structures (chosen from PDB select). Finally, the structures were minimized in order to avoid non-realistic contacts or positions of the residues.

For each structure, after a 1000 steps minimization via the Steepest Descent algorithm, 5-ns MDs NVT simulation in a octahedral water box with explicit solvent and periodic boundary conditions are run using the GROMACS package (version 3.2.1), [START_REF] Lindahl | Gromacs 3.0: A package for molecular simulation and trajectory analysis[END_REF] with the GROMOS96 43Al force field [START_REF] Van Gunsteren | Gromos force field[END_REF].

The simple point charge model SPC is applied, to model water molecules [START_REF] Berendsen | The missing term in effective pair potentials[END_REF]. All bond lengths are constrained by means of the LINCS algorithm [START_REF] Hess | A linear constraint solver for molecular simulations[END_REF]. Electrostatic interactions are treated via PME implementation of the Ewald summation method. Temperature is set to 300 K and controlled by

Berendsen thermostat [START_REF] Berendsen | Molecular dynamics with coupling to an external bath[END_REF]. The timestep is set to 2 fs.We neutralized the N-and C-terminals, since charged terminals form a saline bond that decreases the mobility of the protein, resulting in a less representative sampling of conformational space.

Energy Decomposition Method.

The energy decomposition method is based on the calculation of an interaction matrix M ij on a representative protein structure derived from an MD trajectory [START_REF] Tiana | Understanding the determinants of stability and folding of small globular proteins from their energetics[END_REF][START_REF] Ragona | Determinants of protein stability and folding: Comparative analysis of beta-lactoglobulins and liver basic fatty acid binding protein[END_REF][START_REF] Colacino | The determinants of stability in the human prion protein: Insights into the folding and misfolding from the analysis of the change in the stabilization energy distribution in different condition[END_REF][START_REF] Colacino | Similar folds with different stabilization mechanisms: The cases of prion and doppel proteins[END_REF][START_REF] Morra | Relationship between energy distribution and fold stability: Insights from molecular dynamics simulations of native and mutant proteins[END_REF]. The matrix contains the

A C C E P T E D M A N U S C R I P T

interaction energies between residue pairs, comprising all the nonbonded interresidue atomic energy components (namely, van der Waals and electrostatic couplings between all atoms of two residues).

Solvent effects such as electrostatic shielding and pair non polar terms are are implicitly taken into account by means of the generalized Born approximation (GBSA). In the following we provide a detailed explanation of the method, which was already presented previously in [START_REF] Tiana | Understanding the determinants of stability and folding of small globular proteins from their energetics[END_REF][START_REF] Morra | Relationship between energy distribution and fold stability: Insights from molecular dynamics simulations of native and mutant proteins[END_REF].

The square matrix M ij of non-bonded interactions between all residues i and j in the protein can be diagonalized and re-expressed in terms of eigenvalues and eigenvectors, in the form:

∑ = = N k k j k i k ij w w M 1 λ (1)
where N is the number of amino acids in the protein, k λ is an eigenvalue, and k i w is the i-th component of the associated normalized eigenvector. [START_REF] Bowie | A method to identify protein sequences that fold into a known three-dimensional structure[END_REF] λ is defined as the most negative and the following eigenvalues are labelled in increasing order. In the following we refer to the first eigenvector as the eigenvector corresponding to the eigenvalue 1 λ . The total non-bonded energy E nb is defined as:

∑ ∑ ∑ ∑ = = = = = = = N k k k N k j k k i k N j i N j i ij nb W w w M E 1 1 1 , 1 , λ λ (2)
where

∑ = = 1 , j i j k k i k w w W . If 1 1 W λ is larger than k k W λ for 1 ≠ k j i k w w 1 1 λ
, the sum over i,j of M ij is dominated by the contribution due to the first eigenvalue and eigenvector, such that the total non bonded energy can be approximated by:

1 1 1 , 1 1 1 W w w E E N j i j i nb app nb λ λ = = ≈ ∑ = (3)
The hot spots, or the residues giving a more relevant contribution to the stabilization energy, are defined as those sites whose component is higher than a threshold value t which is calculated as the value corresponding to a normalized vector whose components provide all the same contribution for each site (flat eigenvector). These hot spot values correspond to the "peaks" in the graphical representation of the eigenvector as a function of the sequence. This corresponds to a case in which each residue contributes with the same weight to structural stability. In this approximation the threshold value depends only on the number N of residues in the protein and is calculated as:

N w i 1 1 = for each i [15-19].

A C C E P T E D M A N U S C R I P T

We verified that the most frequented cluster was significantly more populated than the others, not to neglect significant structural deviations captured from other clusters. The Energy Decomposition method was applied on the representative structure of this cluster.

The energy pair decomposition is obtained with the MM-PBSA adapting the algorithm implemented in AMBER8 using the GB approximation (GB model of Onufriev et al. [START_REF] Onufriev | Modification of the generalized born model suitable for macromolecules[END_REF]). The pairwise energy contributions are calculated as sum of gas phase interactions, solvation free energy with GB and hydrophobic term obtained with the LCPO method [START_REF] Weiser | Approximate atomic surfaces from linear combinations of pairwise overlaps (lcpo)[END_REF]. The contact map C ij of a structure is a matrix that describes which residues are in contact in the starting conformation and is defined by looking at Cα atom pairs.

If the distance between any two Cα atoms is below a cutoff value, the corresponding matrix entry is set to l, otherwise it is set to 0. The distance cutoff is set to 0.65 nm. For the sake of homogeneity with the energy matrix, also contacts between nearest neighbors i,j+1 are included. Therefore:

   > ≤ = 5 . 6 0 5 . 6 1 ij ij ij r r C
Upon diagonalization of the matrix C ij , we obtain a set of eigenvalues and eigenvectors. The principal eigenvector, corresponding to the most positive eigenvalue, has all components c i of equal sign, which is also true for the first eigenvector of matrix M ij .

The similarities among energy patterns and between contact map and energy patterns are calculated using the Pearson's correlation coefficient:

To evaluate patterns similarities we also performed a cluster analysis with Matlab 7 release 14. In order to compute the distances among patterns, we used the function PDIST, method CORRELATION:

given an m-by-n data matrix X, which is treated as m (1-by-n) row vectors x 1 , x 2 , ..., x m , the various distances between the vector x r and x s are defined as follows:

r(c, w) = (c i -< c >)(w i -< w >) i ∑ (c i -< c >) 2 i ∑       (w i -< w >) 2 i ∑       d rs = 1 - x r -x r ( ) x s -x s ( ) ′ x r -x r ( ) x r -x r ( ) ′       1 2
x sx s ( ) x s -x s ( )

′       1 2
We then clustered the patterns with the LINKAGE function (default method), which creates a hierarchical cluster tree. showing the peaks of the SE, residues E7 M8 A9 R17 Y18 F19 L20 with VdW spheres

Caption for Table:

Similarity values between the Sequence Eigenvector of the Native structure of YAP65 WW domain.

PDB id: 1k9r and each of the sequences tested in this paper. The similarity measure is based on the calculation of the Pearson's Coefficient. The sequences are reported in the Supp. Mat.

Supplementary Material.

The Supp. Mat. File reports the list and the alignment of the sequences analyzed.

  requires cooperativity and therefore correlation among different sites. The analysis was based on purely physico-chemical and structural properties; hence no previous knowledge of the folding abilities of sequences was used. The procedure entails the evaluation of the Sequence Eigenvector (SE) of each sequence on the target fold by means of the Energy Decomposition Method developed in our group[START_REF] Tiana | Understanding the determinants of stability and folding of small globular proteins from their energetics[END_REF][START_REF] Ragona | Determinants of protein stability and folding: Comparative analysis of beta-lactoglobulins and liver basic fatty acid binding protein[END_REF][START_REF] Colacino | The determinants of stability in the human prion protein: Insights into the folding and misfolding from the analysis of the change in the stabilization energy distribution in different condition[END_REF][START_REF] Colacino | Similar folds with different stabilization mechanisms: The cases of prion and doppel proteins[END_REF][START_REF] Morra | Relationship between energy distribution and fold stability: Insights from molecular dynamics simulations of native and mutant proteins[END_REF]. The Energy Decomposition Method was applied to a set of 32 WW domain sequences, including 8 native sequences (N) and two groups of 8 designed sequences, one obeying both site conservation statistics and a set of selected residue-residue correlations (CC) and the other satisfying only single site conservation (IC). 8 random sequences were also considered for comparison. The similarity between each sequence's SE and the template SE from native WW domain protein 1k9r was shown to correlate with the folding capability of the sequence under examination. Around 70% of the effectively folding sequences has a Pearson's correlation coefficient of at least 0.70 with the template, whereas not folding sequences reach lower values. A threshold of 0.70 for the Pearson's correlation coefficient could therefore be used to classify between folders and non-folders to the WW domain. We also tested a classification criterion based on the comparison of a sequence's SE to the structure CE, recapitulating the structural properties of the fold. The similarity of the sequence's SE and the WW domain structure CE is not a sufficient criterion to distinguish between folding and not folding sequences, since it yields equivalent scores for N, CC and IC sequences. Apparently, the similarity to a structure's CE guarantees the correct placement of hydrophobic and polar residues in the highly buried and exposed sites, respectively. The natural occurring frequency of a given type of residue with specific chemical properties at a given site (such as polar or non polar amino acids) corresponds to the conservation statistics obtained in the multiple sequence alignment. The natural single site distribution,A C C E P T E DM A N U S C R I P T

Figure 1 .

 1 Figure captions: Figure 1. a) Structure of the YAP 65 WW domain (Pdb entry 1k9r). The terminal segment depicted in green was not considered in the calculation. The protein regions corresponding to significant contributions to stability (SE components over threshold) are depicted as blue lines. b) Same as a)

Figure 2 .

 2 Figure 2. Top, SE profile of the 8 native sequences mounted on the template structure. Bottom, SE profile of the WT template sequence.

Figure 3 .

 3 Figure 3. Top, SE profile of the 8 Random sequences, Middle: SE profile of the 8 IC sequences, Bottom: SE profile of the 8 CC sequences.

Figure 4 .

 4 Figure 4. CE profile calculated over residues 8-40 of the experimental structure 1k9r.

Figure 5 .

 5 Figure 5. Pearson correlation coefficients of each sequence's SE with respect to the template SE (x axis) plotted versus the Pearson correlation coefficients of each sequence's SE with respect to the CE of the X ray template structure.

A

  C C E P T E D M A N U S C R I P T A C C E P T E D M A N U S C R I P T A C C E P T E D M A N U S C R I P T

TABLE 1 .

 1 Similarity values between the Sequence Eigenvector of the Native structure of YAP65 WW domain. PDB id: 1k9r and each of the sequences tested in this paper. The similarity measure is based on the calculation of the Pearson's Coefficient. The sequences are reported in the Supp. Mat.

	Sequence	Pearson's

coeff SE- wtSE Clustalw alignment score Wt Pearson coeff. SE-CE

	N33 R19		0.834896 0.390216		42 6	0.914411 0.576323
	N40 Average R	0.708581 0.45		42 9.75	0.907512 0.65
	Average N CC8 CC13 CC14 CC18 CC22 CC24 CC43 CC45 Average CC	0.78 0.708989 0.794943 0.794792 0.62209 0.669008 0.645391 0.558861 0.934383 0.72	43.88 45 18 48 27 24 51 48 39 M A N U S C R I P T 37.5	0.87 0.901651 0.722429 0.878284 0.767779 0.700183 0.879911 0.837438 0.854095 0.82
	IC1 IC5 IC10 IC16 IC23 IC25 IC35 IC41	0.68257 0.613809 0.671559 0.59906 0.627277 0.597115 0.688846 A C C E P T E D 0.742172	42 36 36 27 30 45 48 41	0.931625 0.805095 0.909656 0.808594 0.884833 0.86078 0.925319 0.878261
	Average IC	0.65		39	0.88
	Wt		1		-	0.744945
	R2		0.308312		6	0.774204
	N6		0.815371		39	0.913731
	R5		0.405056		9	0.569796
	N7		0.553304		51	0.825885
	R6		0.706564		9	0.894776
	N8		0.90147		51	0.799291
	R8		0.595654		6	0.559788
	N11		0.821357		45	0.831807
	R10		0.096114		12	0.530359
	N22		0.724998		24	0.89256
	R13		0.71476		21	0.781025
	N28		0.883323		57	0.896717
	R16		0.421356		9	0.512525

A C C E P T E D M A N U S C R I P T A C C E P T E D M A N U S C R I P T