Mondher Bouteraa

Chérif Nouar
email: cherif.nouar@univ-lorraine.fr

Rayleigh-Bénard convection in non-Newtonian Carreau fluids with arbitrary conducting boundaries

Keywords: Rayleigh-Bénard convection, non-Newtonian fluid, shear-thinning fluid, pattern selection

The objective of the present work is to investigate the Rayleigh-Bénard convection in non-Newtonian fluids with arbitrary conducting boundaries. A linear and weakly nonlinear analysis is performed. The rheological behavior of the fluid is described by the Carreau model. As a first step, the critical Rayleigh number and wavenumber for the onset of convection are computed as a function of the ratios ξ b and ξ t of the thermal conductivities of the bottom and top slabs to that of the fluid. In the second step, the preferred convection pattern is determined using an amplitude equation approach. The stability of rolls and squares is investigated as a function of (ξ b , ξ t ) and the rheological parameters.

. For a significant deviation from the critical conditions, the nonlinear convection terms and nonlinear viscous terms become stronger, reducing overall the stability domain of squares. The largest Nusselt number, N u, is obtained for perfectly conducting boundaries. For a given ξ b , ξ t , the stable solution yields the largest Nusselt number. The enhancement of heat transfer due to shear-thinning effects is significantly reduced for poorly heat conducting plates.

Introduction

The problem of Rayleigh-Bénard convection (RBC) in Newtonian and non-Newtonian fluids layer heated from below and cooled from above remains one of the classical problems of fluid dynamics and heat transfer. In spite of intensive studies made in the past and extensive research work undertaken so far to understand the competition between convective structures (rolls, squares and hexagons) which are often influenced by the boundary conditions (see Holmedal et al. [START_REF] Holmedal | Planform selection in Rayleigh-Bénard convection between finite slabs[END_REF]; Clever and Busse [START_REF] Clever | Convection rolls and their instabilities in the presence of a nearly insulating upper boundary[END_REF]) and other parameters such as temperaturedependence of viscosity (see White [START_REF] White | The planforms and onset of convection with a temperaturedependent viscosity[END_REF]; Palm [START_REF] Palm | On the tendency towards hexagonal cells in steady convection[END_REF]; Richter [START_REF] Richter | Experiments on the stability of convection rolls in fluids whose viscosity depends on temperature[END_REF]; Olivier and Booker [START_REF] Oliver | Planform of convection with strongly temperaturedependent viscosity[END_REF]; Busse and Frick [START_REF] Busse | Square-pattern convection in fluids with strongly temperature-dependent viscosity[END_REF]; Jenkins [START_REF] Jenkins | Rolls versus squares in thermal convection of fluids with temperature-dependent viscosity[END_REF]), there are still many outstanding issues that need to be answered. May be one of the most important question to be addressed is the effect of conductive horizontal plates on the heat transfer and the convection patterns. For instance, in geophysical problems and particularly in the context of the Earth's mantle convection, continents and oceans impose different thermal boundary conditions at the top of the mantle: continents act as insulators while a fixed temperature is imposed by oceans. These different thermal boundary conditions affect the convective flow and the heat transport in the Earth's mantle [START_REF] Guillou | On the effect of continents on mantle convection[END_REF]. Actually, in most numerical investigations of RBC, the plates are assumed to be infinitely heat conducting, and a fixed temperature at the boundaries is imposed, while in engineering and geophysical problems as well as in laboratory experiments the boundaries have a finite conductivity. This may lead to a discrepancy between the experimental and the numerical/theoretical results. The ratios ξ b and ξ t between the thermal conductivities of the bottom and top slabs and that of the fluid may have a significant effect and must be taken into account as additional parameters [START_REF] Riahi | Nonlinear convection in a porous layer with finite conducting boundaries[END_REF].

The influence of the thermal conductivity of the boundaries on Rayleigh-Bénard convection was first investigated in the Newtonian case by Busse and Riahi [START_REF] Busse | Nonlinear convection in a layer with nearly insulating boundaries[END_REF] using a weakly nonlinear analysis. They considered the situation where ξ b = ξ t = ξ << 1 and found that the wavelength of convection flow becomes very large in comparison with the height of the layer and only square patterns are stable. This result was confirmed and extended to the fully nonlinear problem by Proctor [START_REF] Proctor | Planform selection by finite amplitude thermal convection between poorly conducting slabs[END_REF] using a 'shallow water theory'. Afterwards, Jenkins and Proctor [START_REF] Jenkins | The transition from roll to square-cell solutions in Rayleigh-Bénard convection[END_REF] determined the critical value of the thermal conductivities ratios ξ b = ξ t = ξ c at which the preferred planform changes from square cell to roll.

For P r > 10, they found that the preferred planform is rolls when ξ > 1, and squares when ξ < 1. Le Gal et al. [START_REF] Gal | Square versus roll pattern at convective threshold[END_REF] carried out experiments to study Rayleigh-Bénard convection in silicone oil confined between two glass plates. So that ξ b = ξ t = ξ = 7. Near the threshold of convective instability, at ǫ < 0.021, where ǫ is the relative distance from the onset of instability, they observed cells of square planforms. But when 0.024 < ǫ < 0.057, the amplitude of two mutually perpendicular roll sets underwent periodic oscillations in antiphase with another; as ǫ was increased and convection became more intense, one set became predominant and then only roll still stable. This experiment was subsequently modified by Le Gal and Croquette [START_REF] Legal | Appearance of a square pattern in a Rayleigh-Bénard experiment[END_REF] : glass was replaced by plexiglass and water was used as the working fluid, so that ξ = 0.4. In contrast to the preceding experiment, squares were observed in a wide range of ǫ values without any signs of destabilization. The authors think that in the first case, the silicone oil behaves as a mixture and the observed features were governed by the thermophoresis .

Although extensive studies have been devoted to understand the influence of the thermal boundary conditions on the Rayleigh-Bénard convection in Newtonian fluids, only a limited number of works have dealt with complex fluids.

In comparison with the Newtonian system, the nonlinearity of the rheological law introduces an additional coupling in the velocity component. Recently, Bouteraa and Nouar [START_REF] Bouteraa | Weakly nonlinear analysis of Rayleigh-Bénard convection in a non-Newtonian fluid between plates of finite conductivity: Influence of shear-thinning effects[END_REF] have investigated the influence of shear-thinning effects on the convection in a horizontal layer of a shear-thinning fluid between two horizontal symmetric plates of finite thermal conductivity. The rheological behavior of the fluid is described by the Carreau model. The authors found that: (i) the characteristic time of instability τ 0 increases significantly when ξ < 1, (ii) the critical value of the shear-thinning degree α c above which the bifurcation becomes subcritical increases with decreasing ξ, and (iii) the critical value ξ c at which the planform changes from square-cell solution (ξ < ξ c ) to two-dimensional roll solution (ξ > ξ c ) decreases with increasing shear-thinning effects.

In some experimental situations, ξ t = ξ b . For Newtonian fluids, Riahi [START_REF] Riahi | Nonlinear convection in a porous layer with finite conducting boundaries[END_REF] [START_REF] Riahi | Nonlinear thermal convection with finite conducting boundaries[END_REF] has studied this problem and demonstrated, using a linear stability analysis of stationary flows the enormous influence of thermal boundary conditions (when ξ b = ξ t ) on the competition between the convection patterns. He found that squares are stable when rolls are unstable and vice versa, and always hexagonal patterns are unstable. No hysteresis effect is found. In addition, Riahi [START_REF] Riahi | Nonlinear thermal convection with finite conducting boundaries[END_REF] has also shown that square planforms are preferred in a bounded region Ω in the (ξ b , ξ t )-space coordinate system and rolls are favored only outside Ω. When P r < 0.025, the region Ω is quite small and disappears as P r = 0. However, for P r > 7, Ω is largest and nearly independent of P r. Using nonlinear developments, Clever and Busse [START_REF] Clever | Convection in a layer heated from below with a nearly insulating boundary[END_REF] [3] demonstrated in the case of stressfree nearly insulating top plate and highly conducting no-slip lower plate, that two-dimensional rolls are stable near the onset, but become unstable at higher Rayleigh number and are replaced by which is called hexharoll convection.

From experimental point of view, Darbouli et al. [START_REF] Darbouli | Rayleigh-Bénard convection for viscoplastic fluids[END_REF] have investigated Rayleigh-Bénard convection for viscoplastic fluids confined in a cylindrical cell. They used two different horizontal plates of finite thermal conductivity. The bottom and upper walls are made respectively of copper alloy and glass. They used distilled water as Newtonian fluid to validate their experimental setup and an aqueous solution of Carbopol 940 as viscoplastic fluid. In these situations, the ratios ξ t and ξ b are estimated to ξ t = 2 and ξ b = 201.6 for both fluids (authors estimated that the solution of Carbopol 940 has the same thermal conductivity than water). Hence, it is no longer possible to rely on the assumption that the plates are held at fixed and uniform temperatures, which corresponds to plates with infinite thermal conductivity.

The purpose of the present work is to study the influence of arbitrary thermal-conducting top and bottom boundaries on nonlinear processes of Rayleigh Bénard convection, and to see the influence of the shear-thinning effect on the preferred flow pattern. The finite conductivity of the slabs remains one explanation for differences between results obtained in experiments and numerical investigations. We hope that our findings will shed new light on the interpretation of the results obtained by Darbouli et al. [START_REF] Darbouli | Rayleigh-Bénard convection for viscoplastic fluids[END_REF] although the fluid used is not only shear-thinning but has also a yield stress.

Physical and mathematical model

General equations and parameters

We consider a horizontal layer of a shear-thinning fluid of height d confined between two horizontal plates that are infinite in extent and which have In the fluid, 0 < ẑ < d, the hydrostatic solution and the temperature profile are:

d P dẑ = -ρ0 ĝ 1 -β T -T0 and Tcond = T0 + ∆ T 1 + Λ/ξ (b) + Λ/ξ (t) 1 2 - ẑd , (1)
where, ĝ is the acceleration due to gravity. Here, the z-axis is directed upwards, with the origin located at the bottom plate. The reference temperature T0 is the temperature in the middle of the fluid layer and ρ0 is the fluid density at T0 .

Here, T0 = T1 -1/2 + Λ/ξ b ∆ Tf , where T1 is the temperature on the outer surface of the bottom plate and ∆ Tf the temperature difference between the top and the bottom of the fluid layer:

∆ Tf = ∆ T /(1 + Λ/ξ t + Λ/ξ (b) ).
The temperature profile in the top and bottom plates are:

Tcond = T0 + ∆ T ξ t + Λ (1 + ξ t /ξ b ) 1 + Λ 2 1 - ξ t ξ b - 1 2 ξ t - ẑd , (2)
d ≤ ẑ ≤ (1 + Λ) d and Tcond = T0 + ∆ T ξ b + Λ (1 + ξ b /ξ t ) 1 2 ξ b - Λ 2 1 -ξ b /ξ t - ẑd , -Λ d ≤ ẑ ≤ 0. ( 3
)
When the bottom and top plates are poor thermal conductors, a large part of ∆ T occurs across the plates, and remains only a small part ∆ Tf of ∆ T , acting as the driving force for the convection. When ∆ Tf , exceeds a critical value, the convection sets in and a so-called convective patterns emerge. The stability of the hydrostatic solution is considered by introducing temperature and pressure perturbation as well as a fluid motion. Boussinesq approximation is adopted, i.e., the temperature dependence of the fluid properties can be neglected except for the temperature induced density difference in the buoyant force that drives the flow. The heat production due to viscosity is neglected. Distances are scaled with d, velocity with κ/ d, where κ is the thermal diffusivity of the fluid, time with d2 /κ, temperature with ∆ Tf , pressure and stresses with κμ 0 / d2 . Using these scales, the dimensionless perturbation equations read:

∂∇ 2 w ∂t -e z • [∇ × ∇ × [(u.∇) u]] = P r ∆ 2 w + Ra P r ∇ 2 H θ - (4) P r [∇ × ∇ × [∇ • (µ -1) γ]] • e z , ∂θ ∂t + (u • ∇) θ = w + ∇ 2 θ , (5)
in the fluid, and

∂ θ ∂t = κ(b,t) p κ ∇ 2 θ , (6)
in the bounding slabs. Here, e z denotes the unit vector in the vertical direction,

w = u(x, t
).e z is the vertical velocity field and θ(x, t) represents the temperature deviations from their values in the conduction state. The temperature perturbation in the slabs is denoted θ(x, t). The position vector x has components x, y, z. The pressure was eliminated by applying twice the rotation to the momentum balance equation. The state of the fluid is controlled by two dimensionless variables. These are the Rayleigh number Ra which measures the ratio between the driving buoyancy and the damping forces, and the Prandtl number P r which represents the ratio between thermal and viscous diffusion times:

Ra = ρ0 ĝ β∆ Tf d3 κ μ0 ; P r = μ0 ρ0 κ . (7)

Rheological model and parameters

In the present study the behavior of the shear-thinning fluid is described by the Carreau model given by:

μ -μ∞ μ0 -μ∞ = 1 + λ2 Γ nc-1 2 with Γ = 1 2 γij γij . (8)
Here, μ0 and μ∞ are the viscosities at low and high shear rate, (n c < 1) the power-law exponent characterizing the shear-thinning regime and λ the characteristic time of the fluid. The location of the transition from the Newtonian plateau to the shear-thinning regime is determined by λ, since 1/ λ defines the characteristic shear rate marking the onset of shear-thinning. The infinite shear viscosity, μ∞ , is generally associated with a breakdown of the fluid, and is frequently significantly smaller (10 -3 to 10 -4 times smaller) than μ0 , see Bird et al.

[21] and Tanner [START_REF] Tanner | Engineering rheology[END_REF]. The ratio μ∞ /μ 0 will be thus neglected in the following.

The dimensionless effective viscosity is then

µ = μ μ0 = 1 + λ 2 Γ nc-1 2 with λ = λ d2 /κ . (9)
The Newtonian behavior, μ = μ0 , is obtained by setting n c = 1 or λ = 0.

For a small amplitude disturbance, the viscosity can be expanded about the hydrostatic solution,

µ = 1 + n c -1 2 λ 2 Γ + 1 2 n c -1 2 n c -3 2 λ 4 Γ 2 + ... (10)
At lowest nonlinear order, a relevant rheological parameter is the 'degree of shear-thinning'

α = dµ dΓ Γ=0 = 1 -n c 2 λ 2 . ( 11
)

Boundary conditions

For the perturbation velocity w, we imposed a realistic boundary conditions on the top and bottom plates: The slabs are rigid, thus enforcing the no-slip boundary conditions (NSBC), which implies

w = Dw = 0 at z = 0, 1. (12)
For thermal boundary conditions, a constant and uniform temperature is assumed on the outer surface of each plate

θ = 0 at z = -Λ, 1 + Λ ( 13
)
The continuity condition of temperature and heat flux at the interface liquid/solid reads

θ = θ at z = 0, 1, (14)
Dθ = ξ (b,t) D θ at z = 0, 1. ( 15
)
Where D ≡ ∂ ∂z and ξ (b,t) = k(b,t) p k .

Remark: Following Holmedal et al. [START_REF] Holmedal | Planform selection in Rayleigh-Bénard convection between finite slabs[END_REF], applying the transformation

u, v, w, p, θ, θb,t , z -→ u, v, -w, p, -θ, -θb,t , 1 -z (16)
to Eqs. ( 4) -( 6), ( 14) and ( 15), it can be shown that the problem is invariant if ξ b , ξ t -→ ξ t , ξ b . Therefore, the preferred convection pattern, its stability and the heat transfer is not modified if the bottom slab is replaced by the top slab and vice versa.

3. Linear stability analysis

Critical conditions and critical modes

In the linear theory, u and θ are assumed infinitesimal. The nonlinear terms in ( 4)-( 5) can be neglected. We obtain :

1 P r ∂∆w ∂t = ∆ 2 w + Ra∆ H θ , (17)
∂θ ∂t = w + ∆θ , (18)
∂ θ ∂t = κ(b,t) p κ ∆ θ. ( 19
)
At this stage, no non-Newtonian effects enter the problem. For equations ( 17)-185 [START_REF] Clever | Convection in a layer heated from below with a nearly insulating boundary[END_REF], we seek a normal mode solution

w(x, y, z, t) = F 11 (z)f (x, y) exp (s t) (20)
θ(x, y, z, t) = G 11 (z)f (x, y) exp (s t) (21)
θ(x, y, z, t) = G11 (z)f (x, y) exp (s t) ,

where, s = s r +is i is a complex number and f (x, y) satisfies the two-dimensional Helmoltz equation ∆ H f = -k 2 f . Here k is the norm of the horizontal wavenumber k. Substituting (20)-( 22) into ( 17)-( 19), leads to the differential equations

s P r -1 D 2 -k 2 F 11 = -k 2 RaG 11 + D 2 -k 2 2 F 11 , (23)
s G 11 = F 11 + (D 2 -k 2 )G 11 , (24)
s G11 = κ(b,t) p κ (D 2 -k 2 ) G11 . ( 25
)
The boundary conditions are:

F 11 = DF 11 = 0 at z = 0, 1 (26)
G11 = 0 at z = -Λ, 1 + Λ, (27)
G 11 = G11 at z = 0, 1 (28)
DG 11 = ξ (b,t) D G11 at z = 0, 1 . (29)
Denoting by the superscript * the complex conjugate, multiplying (23) by F * 11 , (24) by -k 2 RaG * 11 , (25) by -k 2 Ra κ/κ b,t ξ (b,t) G * 11 , summing and integrating over the interval [-Λ, 1 + Λ], we have

s 1+Λ -Λ |DF 11 | 2 + k 2 |F 11 | 2 + k 2 Ra |G 11 | 2 + k 2 Ra ξ (b,t) κ κ(b,t) p G11 2 dz = k 2 Ra 1+Λ -Λ (F 11 G * 11 + F * 11 G 11 ) dz -k 2 Ra 1+Λ -Λ |DG 11 | 2 + k 2 |G 11 | 2 dz (30) -k 2 Raξ (b,t) 1+Λ -Λ |D G11 | 2 + k 2 | G11 | 2 dz, with, F 11 = 0 as z ∈ [-Λ, 0] ∪ [1, 1 + Λ, G 11 = 0 as z ∈ [-Λ, 0[∪]1, 1 + Λ] and
G11 = 0 when z ∈]0, 1[. The right hand side of equation ( 30) is real, therefore s i = 0: the linearized equations ( 23) -( 25) satisfy the principle of exchange of stability. The conduction state loses its stability to a convective stationary state in a smooth way.

The set of differential equations ( 23)-( 25) is an eigenvalue problem where s is the eigenvalue and X 11 = (F 11 , G 11 , G11 ) the eigenvector. It can be written

s M • X 11 = L • X 11 . (31)
Actually, Eq. 25 can be solved analytically:

G11 (z) = G 11 (z = 1) sinh k(t) (1 + Λ -z) sinh k(t) Λ ; 1 ≤ z ≤ 1 + Λ , (32)
G11 (z) = G 11 (z = 0) sinh k(b) (Λ + z) sinh k(b) Λ ; -Λ ≤ z ≤ 0 , (33)
with k(b,t) = k 2 + sκ/κ (b,t) p

. Hence, the eigenvalue problem ( 23)-( 25) can be restrained to the fluid domain, i.e. Eqs ( 23) -( 24), with the boundary conditions

DG 11 = ±ξ (b,t) k(b,t) G 11 coth k(b,t) Λ ; z = 0, 1. (34)
The eigenvalue problem ( 23)-( 24) with the boundary conditions (34) is solved using a Chebyshev collocation method. By setting s = 0, one obtains the marginal stability curve Ra(k). The minimum of this curve gives the critical Rayleigh number Ra c and the corresponding wavenumber k c . Ra c allows us to determine the critical ∆ Tf at which the system changes from the hydrostatic state to the state of cellular motion and k c provides information about the horizontal periodicity of the patterns at the onset of convection. We recall that the transition from the rest state to the convective state is independent of the Prandtl number. Mathematically, the decrease of Ra c is caused by the weakening of the thermal boundary conditions ( 28), [START_REF] Heutmaker | Convection patterns: time evolution of the wave-vector field[END_REF] for the temperature fluctuations θ as ξ b or ξ t or both decrease from 10 3 to 10 -3 . These results are in very good quantitative agreement with those obtained by Sparrow et al. [START_REF] Sparrow | Thermal instability in a horizontal fluid layer: effect of boundary conditions and a non-linear temperature profile[END_REF], Proctor [START_REF] Proctor | Planform selection by finite amplitude thermal convection between poorly conducting slabs[END_REF], Jenkins and

Proctor [START_REF] Jenkins | The transition from roll to square-cell solutions in Rayleigh-Bénard convection[END_REF], Carriere et al. [START_REF] Carriere | Wavelength selection in Rayleigh-Bénard convection between horizontal boundaries of finite conductivity[END_REF] and Cerisier et al. [START_REF] Cerisier | Thermal influence of boundaries on the onset of Rayleigh-Bénard convection[END_REF]. From a physical point of view, it means that when a thermal fluctuations occurs in the fluid close to infinitely conducting boundaries, it easily relaxes. However, when either ξ b or ξ t or both have a low conductivity, the wall temperature fluctuations can persist and become a supplementary source of instability for the bulk, and then decreasing Ra c . The experimental values of Ra c and k c found by Kebiche et al.

[26] and Darbouli et al. [START_REF] Darbouli | Rayleigh-Bénard convection for viscoplastic fluids[END_REF] are much larger than those predicted by the linear theory when the finite thermal conductivity of the plates is taken into account.

A possible explanation of this discrepancy is given in §5.

In Figure 2, the eigenfunctions, F 11 , G 11 and G11 are presented for different combinations of the conductivities of the upper and lower plates. Note that the symmetry of F 11 is little influenced by the thermal boundary conditions.

As indicated above, the temperature fluctuation doesn't relax at the wall when ξ t,b < 1. With decreasing ξ b,t , the fluctuation of the temperature gradient decreases. When, ξ t,b << 1, ∂θ/∂z → 0 and the situation of fixed heat flux at the boundary is recovered. Contours of the temperature perturbation are represented in Fig. 3, for three situations. The first one (Fig. 3a) represents the classical case where the boundaries are "good heat conductor". The second situation (Fig. 3b) corresponds to the experimental device used by Kebiche et al. [START_REF] Kebiche | Experimental investigation of the Rayleigh-Bénard convection in a yield stress fluids[END_REF]. The horizontal plates are made of Polycarbonate with ξ t = ξ b = 0.25.

The third situation (Fig. 3c) corresponds to the experimental device used by Darbouli et al. [START_REF] Darbouli | Rayleigh-Bénard convection for viscoplastic fluids[END_REF] where ξ t = 2 and ξ b = 201.6. As we can see, the coupling between the temperature perturbation θ in the fluid and that in the horizontal plates θ is specific for each situation. With decreasing ξ b,t , the horizontal scale of the convection motion increases and the vertical diffusion rate becomes greater than the horizontal one.

Characteristic time of the instability

Near the onset of convection, the growth rate Re(s) of the perturbation may be approximated using Taylor expansion,

s = ǫ τ 0 + O ǫ 2 with ǫ = Ra -Ra c Ra c . ( 35
)
The characteristic time τ 0 of instability to grow is obtained by evaluating (ds/dǫ) ǫ=0 ; s is calculated for different values of ǫ near ǫ = 0. Figure 4 shows the variation of τ 0 as a function of ξ (b,t) . As it can be observed, the characteristic time of instability increases significantly as the thermal conductivities ratios decrease. This can be explained by the fact that the thermal disturbance evolves on a very long time scale in the solid compared to that of the fluid when the slabs are a very poor heat conductors compared to the fluid. The following asymptotic limits are found: (i) τ 0 = 0.055 for perfectly heat conducting boundaries (ξ b > 100 and ξ t > 100); (ii) τ 0 = 0.1 for one insulating boundary and the other one perfectly conducting; (iii) τ 0 → ∞ for insulating boundaries. (a)

ξ t = 1 ξ t = 10 3 ξ t = 10 -1 ξ t = ξ b ξ t = 10 - 3 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 0 0.5 1 1.5 2 2.5 3 3.5 k c (b)
ξ b ξ t = 10 -3 ξ t = 10 3 ξ t = 1 ξ t = 10 -1 ξ t = ξ b
ξ t = ξ b ξ t = 10 -1 ξ t = 10 -3 ξ t = 1 ξ t = 10 3

Weakly nonlinear stability analysis: Pattern selection

As it is known, linear stability analysis yields the onset of instability and the critical wave number. It is unable to give unique predictions for the form or evolution of the resulting patterns at a finite distance above the onset of convection [START_REF] Heutmaker | Convection patterns: time evolution of the wave-vector field[END_REF]; this is so because the eigenvalue problem is degenerate, since to one eigenvalue Ra there corresponds an infinite number of possible patterns with the same wave number k [START_REF] Lebon | Understanding Non-equilibrium Thermodynamics[END_REF]. The pattern selection is determined by the nonlinear terms. A weakly nonlinear analysis is adopted as a first approach to investigate nonlinear effects in the competition between rolls and squares [START_REF] Bouteraa | Weakly nonlinear analysis of Rayleigh-Bénard convection in shear-thinning fluids: nature of the bifurcation and pattern selection[END_REF], [START_REF] Bouteraa | Weakly nonlinear analysis of Rayleigh-Bénard convection in a non-Newtonian fluid between plates of finite conductivity: Influence of shear-thinning effects[END_REF]. Hexagons configuration is not considered here. Further calculations show that this three-dimensional planform is unstable when the temperature dependence of the viscosity is not taken into account [1][11][18] [START_REF] Busse | Nonlinear convection in a layer with nearly insulating boundaries[END_REF]. The wide region, Φ, corresponds to the Newtonian case and can be bounded approximately by the lines ξ b + ξ t ≤ 2 and ξ t = ξ b + 0.7. The ascertainment that squares are stable in this domain can be rationalized as follows. When the walls are poor heat conductors compared to the fluid, the convection intensity is low and the vertical velocity of the liquid w is small [START_REF] Verzicco | A comparison of turbulent thermal convection between conditions of constant temperature and constant heat flux[END_REF]. So, this reflects a weak interaction between sets of rolls that constitute squares knowing that no two rolls predominates the other. As a consequence, only squares planforms are stable.

By introducing the shear-thinning character, the convection intensity increases [START_REF] Pierre | Experimental investigation of natural convection heat transfer in confined space for non-Newtonian fluid[END_REF] [34] [35][36], thereby increasing the interaction between the two sets of rolls that constitute squares. As a consequence, one of the two sets of rolls nonlinearly damps out the other and the squares become unstable contrary to rolls which become the preferred form of convection [START_REF] Bouteraa | Weakly nonlinear analysis of Rayleigh-Bénard convection in shear-thinning fluids: nature of the bifurcation and pattern selection[END_REF]. As illustrated in Fig. 5, the region Φ decreases with increasing shear-thinning effects.

With increasing the reduced Rayleigh number ǫ = (Ra -Ra c )/Ra c , the nonlinear convection terms and non linear viscous terms become stronger, reducing the stability region of squares as illustrated in Figs. 6 and7.

These results are consistent with the maximum heat-transfer principle, "the only stable solution is the one of maximum heat transport" [START_REF] Malkus | Finite amplitude cellular convection[END_REF]. The Nusselt number is defined by

N u = 1 - ∂ θ ∂z z=0 . (36)
The overbar denotes the horizontal average.

Rematk : Deviating from the critical conditions, s > 0, the temperature profile in the plates, and the boundary conditions at the interface fluid/solid, involve a parameter denoted k which depends on s κ(b,t) p /κ , see for instance Eqs.

(32)- [START_REF] Liang | Experiments on buoyancy-driven convection in non-Newtonian fluids[END_REF]. For the computations at ǫ > 0, we have assumed as Chapman and

Proctor [START_REF] Sparrow | Thermal instability in a horizontal fluid layer: effect of boundary conditions and a non-linear temperature profile[END_REF] and Carriere et al. [START_REF] Carriere | Wavelength selection in Rayleigh-Bénard convection between horizontal boundaries of finite conductivity[END_REF] that κ(b,t) Therefore, the ratio of the thermal diffusivity of the plate, κp , to that of the fluid, κ, is mainly governed by the thermal conductivities ratio ξ.

Conclusion

In the present study we have investigated a linear and a weakly nonlinear stability analysis of a horizontal Carreau fluid layer bounded by two horizontal plates of arbitrary conductivity. Symmetry properties show that the the results are unchanged if the upper slab is replaced by the lower one and vice-versa.

The critical Rayleigh number and wavenumber associated to the onset of convection are determined as a function of the conductivities ratios ξ t and ξ b . When either ξ b or ξ t or both are low, the wavelength of the convection becomes larger and the critical Rayleigh number Ra c is reduced.

Stability analysis of the stationary solutions of the amplitude equations indicate that near the onset, the convection pattern is either a two-dimensional roll or a square. The influence of shear-thinning effects on the stability-domain of squares in the plane (ξ t , ξ b ) is highlighted. The results are in agreement with the maximum heat transfer principle.

The heat transfer and the convection intensity are strongly reduced when the plates are poor heat conductors. It is therefore more difficult to detect experimentally the onset of convection than in the case of 'good heat conductors'. This could explain the fact that Ra c found by Kebiche et al. [START_REF] Kebiche | Etude expérimentale de l'instabilité de Rayleigh-Bénard dans les fluides non-Newtoniens[END_REF] and Darbouli et al. [START_REF] Darbouli | Rayleigh-Bénard convection for viscoplastic fluids[END_REF] is larger that that predicted by the linear theory.

Finally, we think that the present findings are useful to experimental observations since the finite conductivity of the walls plays a primordial role in the Rayleigh-Bénard convection and it can be a source of divergence between experimental observations and theoretical predictions.

Subsequent developments of the present work will concern the secondary instabilities of the primary pattern described here. Details on the secondary instabilities can be found in Ref. [START_REF] Hoyle | Pattern formation. An introduction to methods[END_REF].

  a thickness Λ d, where Λ is of order unity. The outer surface of the bottom and top plates are kept at constant temperatures respectively T0 + ∆ T /2 and T0 -∆ T /2, with ∆ T > 0. The fluid has density ρ, thermal conductivity k, thermal coefficient expansion (at constant pressure) β and viscosity μ0 at zero shear rate. The top slab has a thermal conductivity kt p and a thermal diffusivity κt p . The corresponding quantities for the bottom slab are denoted kb p and κb p . Here and in what follows, (t) and (b) refer to the top and bottom and the quantities with hat (.) are dimensional. Because of the thermal expansion, the temperature difference between the two plates, induces a vertical density stratification. Heavy cold fluid is above a light warm fluid. For small ∆ T , the fluid remains motionless and the heat is transferred by conduction, with a linear temperature profile across the fluid layer.

Figure 1

 1 displays the variation of Ra c and k c as a function of the ratio ξ b for different values of ξ t . The dimensionless thickness of the plates is fixed at Λ = 1. As expected, the critical conditions are the same if the upper plate is exchanged with the lower one. We recover the following asymptotic limits: (i) Ra c = 1707.7 and k c = 3.12 for two perfectly conducting boundaries, (ii) Ra c = 720 and k c → 0 for two insulating boundaries and (iii) Ra c = 1296 and k c = 2.55 for one insulating boundary and the other one perfectly conducting. When either ξ b or ξ t or both are low, the wavelength of the convection becomes larger and the critical Rayleigh number Ra c is reduced.

Figure 1 :

 1 Figure 1: Critical Rayleigh number (a) and critical wavenumber (b) as function of ξ b for different values of ξ t . The thick black line corresponds to the case ξ t = ξ b . The symbols are the theoretical values of Rac and kc corresponding to experimental devices used in the literature: ( ) Kebiche et al. [27], (•) Darbouli et al. [20] and ( * ) Hassan et al [28].

Figure 2 :

 2 Figure 2: Eigenfunctions (a) F 11 (vertical velocity) and (b) G 11 , G11 (temperature perturbation) for different values of ξ b , ξ t : (1) ξ b = ξ t = 0.1; (2) ξ b = ξ t = 0.25; (3) ξ b = 201.6, ξ t = 2; (4) ξ b = ξ t = 10 3 .

Figure 3 :

 3 Figure 3: Temperature perturbation contours for three Rayleigh-Bénard situations : (a) ξ b = ξ t = 1000 "good heat conductor", (b) ξ b = ξ t = 0.25 corresponding to Kebiche et al.[27] experimental device and (c) ξ b = 201.6 ; ξ t = 0.25 corresponding to Darbouli et al. [20] experimental device.

Figure 4 :

 4 Figure 4: Characteristic time of instability vs thermal conductivities ratio ξ b at Pr = 10 in the case of ξ t = ξ b ; ξ t = 10 -3 ; ξ t = 10 -1 ; ξ t = 1 and ξ t = 10 3 .

Figure 5

 5 Figure 5 shows, at critical conditions, the stability regions of square planforms in the (ξ b , ξ t )-space coordinate system at P r = 10 for a Newtonian fluid and Carreau fluids with n = 0.5 and different values of the dimensionless constant time λ. Squares are stable in the region bounded by the axis and the stability curve. Rolls are stable outside this region. In agreement with Riahi [18], stability curves are symmetric about the line ξ b = ξ t .

pFigure 5 :

 5 Figure 5: Stability boundary of squares in the (ξ b , ξ t )-space coordinate system at P r = 10 for a Newtonian fluid (λ = 0) and Carreau fluids with n = 0.5 two values of λ: domain Ψ, λ = 0.01 and domain Γ, λ = 0.02 .

Figure 8

 8 Figure8shows the variation of the Nusselt number as a function of ǫ for different thermal boundary considered in laboratory experiments. The largest Nusselt number is obtained for perfectly heat conducting slabs. For ξ t << 1 and ξ b << 1, Nu is weakly larger for square-cell solution than for two-dimensional roll solution, but the difference is systematic. Shear-thinning effects reduce the viscous friction at the wall, which leads to an increase of the convection intensity and therefore of the Nusselt number as it is illustrated by Fig.9. However, the enhancement of the heat transfer by shear-thinning effects is significantly reduced for poorly heat conducting plates.

Figure 6 :

 6 Figure 6: Stability boundary of squares in the (ξ b , ξ t )-space coordinate system for a Newtonian fluid at P r = 10 and three different values of the reduced Rayleigh number ǫ. (thick line) ǫ = 0; (thin line) ǫ = 0.05; (dotted line) ǫ = 0.1.

Figure 7 :

 7 Figure 7: Stability boundary of squares in the (ξ b , ξ t )-space coordinate system for a Carreau fluid with n = 0.5 and λ = 0.01 at P r = 10 and two different values of ǫ.

Figure 8 :

 8 Figure 8: Nusselt number versus reduced the reduced Rayleigh number, ǫ, for a Newtonian fluid with different thermal boundary conditions used in laboratory experiments: (1) ξ t = ξ b = 0.25 (Kebiche et al. [26]); (2) ξ t = ξ b = 0.4 (LeGal and Croquette [15]; (3) ξ t = 2 and ξ b = 201.6 (Darbouli et al. [20]); (4) ξ t = ξ b = 10 3 . (dashed line): the planform convection is a square; (continuous line): the planform convection is a two dimensional roll.

Figure 9 :

 9 Figure 9: N u versus ǫ. Influence of shear-thinning effects for two different thermal boundary conditions: Curves (1) and (2) correspond to ξ t = 2 and ξ b = 201.6; curves (3) and (4) correspond to ξ t = ξ b = 1000. Newtonian fluid: curves (1) and (2); Carreau fluid with n = 0.5, λ = 0.01, curves (2) and (4).

Acknowledgments

The authors gratefully acknowledge the financial support of the French National Research Agency, under the grant ANR-10-BLAN-925-01.


Francesca Troilo

Daniela Bonetti

Carlo Camilloni

Angelo Toto

Sonia Longhi

Maurizio Brunori
email: maurizio.brunori@uniroma1.it

Stefano Gianni
email: stefano.gianni@uniroma1.it

The folding mechanism of the SH3 domain from Grb2

SH3 domains are small protein modules involved in the regulation of important cellular pathways. These domains mediate protein-protein interactions recognizing motifs rich in proline on the target protein. The SH3 domain from Grb2 (Grb2-SH3) presents the typical structure of an SH3 domain composed of two-three stranded antiparallel β-sheets orthogonally packed onto each other, to form a single hydrophobic core. Grb2 interacts, via SH3 domain, with Gab2, a scaffolding disordered protein, triggering some key metabolic pathways involved in cell death and differentiation. In this work we report a mutational analysis (Φ-value analysis) of the folding pathway of Grb2-SH3 that, coupled with molecular dynamic simulations, allows us to asses the structure of the transition state and the mechanism of folding of this domain. Data suggest that Grb2-SH3 folds via a native-like, diffused transition state with a concurrent formation of native-like secondary and tertiary structure (nucleation-condensation mechanism) and without the accumulation of folding intermediates. The comparison between our data and previous folding studies on SH3 domains belonging to other proteins, highlights that proteins of this class may fold via alternative pathways, stabilized by different nuclei leading or not to accumulation of folding intermediates. This comparative analysis suggests that the alternative folding pathways for this class of SH3 domains can be selectively regulated by the specific aminoacid sequences.

INTRODUCTION

One of the most informative approaches to address the folding mechanism of globular proteins is to compare experiments performed on homologous proteins [START_REF] Calosci | Comparison of successive transition states for folding reveals alternative early folding pathways of two homologous proteins[END_REF][START_REF] Capaldi | Ultrarapid mixing experiments reveal that Im7 folds via an on-pathway intermediate[END_REF][START_REF] Clarke | Folding studies of Ig-like beta-sandwich proteins suggest they share a common folding pathway[END_REF][START_REF] Friel | Structural analysis of the ratelimiting transition states in the folding of Im7 and Im9: similarities and differences in the folding of homologous proteins[END_REF][START_REF] Gianni | Unifying features in proteinfolding mechanisms[END_REF][START_REF] Zarrine-Afsar | The family feud: do proteins with similar structures fold via the same pathway?[END_REF] . In fact, by describing the folding of proteins sharing the same topology while displaying a different sequence, it is theoretically possible to draw some general rules on the basic principles governing folding. Comparative folding studies have been previously reported for example on the colicin immunity proteins Im7 and Im9 [START_REF] Capaldi | Ultrarapid mixing experiments reveal that Im7 folds via an on-pathway intermediate[END_REF][START_REF] Friel | Structural analysis of the ratelimiting transition states in the folding of Im7 and Im9: similarities and differences in the folding of homologous proteins[END_REF][START_REF] Capaldi | Im7 folding mechanism: misfolding on a path to the native state[END_REF][START_REF] Ferguson | Rapid folding with and without populated intermediates in the homologous four-helix proteins Im7 and Im9[END_REF] , on the immunoglobulin domains [START_REF] Clarke | Folding studies of Ig-like beta-sandwich proteins suggest they share a common folding pathway[END_REF] , on c-type cytochromes [START_REF] Gianni | Parallel pathways in cytochrome c(551) folding[END_REF][10][11] , on homeodomain-like proteins [START_REF] Gianni | Unifying features in proteinfolding mechanisms[END_REF][START_REF] White | Simulation and experiment conspire to reveal cryptic intermediates and a slide from the nucleation-condensation to framework mechanism of folding[END_REF] , on PDZ domains [START_REF] Calosci | Comparison of successive transition states for folding reveals alternative early folding pathways of two homologous proteins[END_REF][START_REF] Chi | A conserved folding mechanism for PDZ domains[END_REF][14][15] and others. Whilst all these studies suggest that the overall general features of folding are by-and-large defined by protein topology, it appears that a closer look at the folding pathway of the different homologues appears to highlight some features specific for each globular protein.

In the context of comparative folding studies, the SH3 domain represents a debated system. In fact, whilst earlier comparison between the src and the spectrin SH3 domains suggested this class of proteins to fold via a robust two-state mechanism characterized by a polarized and highly conserved transition state [16][17][18] , studies on the Sso7d domain revealed an additional complexity 19 . Indeed, while displaying less than 10% sequence identity, the Sso7d protein shares the typical topology of SH3 domains, except for the last β-strand that is a small α-helix in Sso7d. Interestingly, experimental and computational comparison of the folding of Sso7d with SH3 domains, revealed a substantial shift in the folding nucleus from the third to the second β-hairpin 19 . This finding highlighted that, even if protein topology plays a major role in the selection of the folding pathways, the specific nature of the interactions stabilizing the protein is still critical to describe folding mechanisms. Furthermore, it is of interest to note that recent studies highlighted how the folding of SH3 may also occur via a multi-state scenario, with accumulation of intermediates characterized at equilibrium [20][21][START_REF] Dasgupta | Four-state folding of a SH3 domain: saltinduced modulation of the stabilities of the intermediates and native state[END_REF] . It appears therefore that even for a deeply investigated protein system, such as the SH3 domain, folding demands a careful study to be fully understood.

The SH3 domain from Grb2 (Grb2-SH3) corresponds to the typical structure of an SH3 domain composed of two three-stranded antiparallel β-sheets orthogonally packed onto each other, to form a single hydrophobic core 23,24 . Physiologically, the domain is involved in binding a proline rich stretch of amino acids of Gab2 (encompassing residues 503 to 524), with this interaction triggering some key metabolic pathways involved in cell death and differentiation. From the perspective of its primary structure, it is interesting to note that Grb2-SH3 displays a similar degree of sequence identity towards Sso7d and spectrin SH3 (Figure 1), posing this system as an interesting candidate to understand further the folding mechanism of this highly studied class of proteins.

Here we present the characterization of the folding of the Grb2-SH3 domain. By carrying out kinetic experiments on 23 site-directed variants in combination with restrained molecular dynamics simulations, we present the structure of the main folding transition state. The transition state is stabilized by contacts involving both the first β-hairpin and the N-and C-termini of the protein, a finding which appears different from what previously observed for src, spectrin and fyn SH3 16,18,[START_REF] Northeyk | Hydrophobic core packing in the SH3 domain folding transition state[END_REF]26 . Furthermore, we present evidence that this protein folds via a nucleation-condensation mechanism, with a diffused, rather than structurally polarized transition state. The data are discussed in the context of previous work on other SH3 domains.

EXPERIMENTAL AND THEORETICAL METHODS

Site-Directed Mutagenesis

C-SH3 domain of Grb2 was subcloned in a pET28b+ plasmid vector. The constructs encoding the site directed variants of SH3 were obtained using the gene encoding Grb2-SH3 wt as a template to perform site-directed mutagenesis using the QuickChange Lightning Site-Directed Mutagenesis kit (Agilent technologies) according to the manufactorer's instructions. All mutations are conservative All mutations were confirmed by DNA sequencing.

Protein expression and purification

The C-SH3 domain of Grb2 wt and all the site directed variants were expressed in E. coli cells BL21 (DE3). Bacterial cells were grown in LB medium, containing 30 µg/ml of kanamycin, at 37°C until OD 600 = 0.7 -0.8 and then protein expression was induced with 1mM IPTG. After induction cells were grown at 37°C over night and then collected by centrifugation.

To purify the protein, the bacterial pellet was resuspended in buffer 50 mM TrisHCl, 0.5 M NaCl, pH 7.5 with the addition of antiprotease tablet (Complete EDTA-free, Roche), then sonicated and centrifuged. The soluble fraction from bacterial lysate was loaded onto a nickel-charged HisTrap Chelating HP (GE Healthcare) column equilibrated with 50 mM TrisHCl, 0.5 M NaCl, pH 7.5. The protein was then eluted with a gradient from 0 to 1 M imidazole by using an AKTA-prime system. Fractions containing the protein were collected and the buffer was exchanged to 25 mM Hepes pH 7.5 100 mM potassium acetate by using a HiTrap Desalting column (GE Healthcare). The purity of the protein was analyzed through SDS-page.

Protein concentration was estimated by measuring the absorbance of tryptophan residue at 280nm and calculated through the Lambert-Beer equation.

Equilibrium experiments

Equilibrium unfolding experiments were performed on a Fluoromax single photon counting spectrofluorometer (Jobin-Yvon, NJ, USA). C-SH3 protein and all the site directed variants, at a constant concentration of 3 µM, was excited at 280 nm and emission spectra were recorded between 300 and 400 nm, at increasing denaturant (urea) concentration. Experiments were performed at 25°C, using a quartz cuvette with a path length of 1 cm, in buffer 50 mM sodium phosphate buffer at pH 7.2.

Stopped-flow folding experiments

Unfolding and refolding kinetics experiments were carried out on a single-mixing SX-18 stopped-flow instrument (Applied Photophysics), monitoring the change of fluorescence emission. The experiments were performed at 25°C in buffer 50 mM sodium phosphate pH 7.2, by using urea as the denaturant. The excitation wavelength used was 280 nm and the fluorescence emission light was recorded by using a 320 nm cut-off glass filter. For each denaturant concentration usually 5 individual traces were averaged. The final concentration of Grb2-SH3 and all the variants was typically 1 µM. In all cases the fluorescence time courses obtained was satisfactorily fitted by using a single exponential equation.

Molecular Dynamics Simulations

Molecular dynamics simulations of SH3 were performed using the CHARMM22* force field 27 with the TIP3P water model 28 . All the simulations were run using GROMACS 29 and PLUMED2 30 . A time step of 2 fs was used together with LINCS constraints [START_REF] Hess | P-lincs: A parallel linear constraint solver for molecular simulation[END_REF] . Van der Waals and Coulomb interactions were implemented with a cutoff at 0.9 nm, and long-range electrostatic effects were treated with the particle mesh Ewald method on a grid with a mesh of 0.1 nm.

A standard 200 ns molecular dynamics simulation at 300 K was performed as a reference for the native state ensemble. The starting conformation was taken from an available X-Ray structure (PDB code 2VWF 23 ) and solvated with 4531 water molecules and 4 sodium ions.

The transition state ensemble was determined following a standard procedure based on the interpretation of Φ value analysis in terms of fraction of native contacts.

Briefly, given a set of experimental Φ values, a pseudo energy term has been added to the force field as the squared difference between experimental and simulated Φ values in order to maximize the agreement with the experimental value while keeping the simulation stable. Given two residues that are not nearest neighbours, the native contacts between them are defined as the number of heavy side-chain atoms located within 0.65 nm in the native structure. The Φ value for a residue i is calculated from the fraction of native contacts that it makes in a given conformation. With this approach only Φ values between 0 and 1 can be incorporated as structural restraints.

The transition state ensemble was generated using 1000 cycles of simulated annealing. Each cycle is 200 ps long, in which the temperature is varied between 300 K and 400 K. Only the structures sampled at the reference temperature are retained for further analysis, resulting in TSE of ~6000 conformations.

RESULTS

The kinetic folding mechanism of wild type Grb2-SH3

In order to characterize the folding mechanism of Grb2-SH3 we initially conducted experiments on the wild type protein. Urea-induced equilibrium denaturation of Grb2-SH3 measured at 25°C, pH 7.2 in 50 mM sodium phosphate buffer by decrease in Trp emission is reported in Figure 2. The observed transition is consistent with a simple two-state behavior, suggesting the absence of stable equilibrium intermediate(s) 32 .

The unfolding free energy in water derived from two-state analysis is 3.1 kcal mol -1 displaying an m D-N value of 0.73 kcal mol -1 M -1 . This value, which is proportional to the change in accessible surface area upon unfolding, is consistent with what expected from a protein of 56 amino acids [START_REF] Myers | Denaturant m values and heat capacity changes: relation to changes in accessible surface areas of protein unfolding[END_REF] .

The folding and unfolding kinetics of Grb2-SH3 were measured by stopped-flow fluorimetry. As expected for a two-state folder, under all investigated conditions, folding and unfolding time courses were consistent with a single exponential decay. Furthermore, in analogy to what previously observed on other SH3 domains, the urea dependence of the observed rate constant (k obs ) on urea concentration conforms to a V-shaped chevron plot (Figure 3), a typical signature of two state folding 32 . Since Grb2-SH3 unfolding displays a low cooperativity, with an m D-N of 0.73 kcal kcal mol -1 M -1 an accurate determination of the folding parameters from each independent experiment is complicated. Therefore, to decrease the fitting error and, at the same time, to test the robustness of two-state folding of Grb2-SH3, equilibrium

ACS Paragon Plus Environment

The Journal of Physical Chemistry and kinetic experiments were fitted globally to the following equations:

Equilibrium: Y obs = Y N + Y D e (m D-N ([urea]-[urea ] 1/2 ) 1+ e (m D-N ([urea]-[urea] 1/2 )
Kinetics:

k obs = k F e (-m F [urea ]) + k U e (-m U [urea]) ; m D-N = m F + m U
with shared m D-N values. The fitting parameters calculated from the global analysis are reported in Table 1.

The structure of the folding transition state of Grb2-SH3

In order to characterize the transition state of folding of Grb2-SH3, we carried out a Φ value analysis 34,35 , by producing 23 site directed variants. The Φ value is then calculated by dividing the effect of the substitution on the activation free energy by that of the stability of the native structure. The conservative variants were designed and the analysis carried out using the standard rules of Φ value analysis, as formalized previously 36 .

Unfolding and folding of all the variants were measured both at equilibrium, by urea induced denaturation, and by kinetics, using the stopped-flow fluorimetry. In all cases, in analogy to what observed for wild type Grb2-SH3, folding and unfolding kinetics were consistent with a single exponential decay. Figure 4 shows the equilibrium and kinetic experiments carried out on each site directed variant. In all cases, data were consistent with a two-state scenario, indicating that Grb2-SH3 folds via a robust mechanism, which does not involve any transient folding intermediates.

To determine the structure of the folding transition state of Grb2-SH3, we used the experimentally measured Φ values as restrains in molecular dynamics simulations.

This method, which has been previously used and validated on several different protein systems [START_REF] Calosci | Comparison of successive transition states for folding reveals alternative early folding pathways of two homologous proteins[END_REF][37][START_REF] Gianni | Structural characterization of a misfolded intermediate populated during the folding process of a PDZ domain[END_REF][39][40][41][START_REF] Gsponer | Determination of an ensemble of structures representing the intermediate state of the bacterial immunity protein Im7[END_REF] , is based on the incorporation of the Φ values as biases on the fraction of formed native contacts in a molecular dynamics simulation trajectory (cf.

Methods).

The structure of the folding transition state of Grb2-SH3, together with the associated contact map, is reported in Figure 5A. It is evident that the protein seems to fold via a native-like transition state that is characterized by the formation of the first β-hairpin, together with a consolidation of the interaction between the N-and C-termini of the protein. Structure gradually tapers off, with the region encompassing the β2-β3

interaction being the most disordered of the ensemble. The structural features of the transition state of folding of Grb2-SH3, in comparison to those previously depicted for other SH3 domains are analysed in the discussion section.

On the basis of the Φ value analysis of src and spectrin SH3 16,18 , it has been previously suggested the structure of the transition state of SH3 domains to be highly polarized. To test this hypothesis for Grb2-SH3, we analysed the Bronsted plot of this protein 43 . In fact, whilst a diffused native like structure is expected to return linear Bronsted plots, a polarized transition state is more likely to yield a scatter in the Bronsted plot, with only some positions playing a key role in stabilizing its structure Bronsted plot, suggesting this protein to fold via a native-like diffused, rather than polarized, transition state. This finding appears consistent with a nucleationcondensation mechanism [START_REF] Abkevich | Specific nucleus as the transition state for protein folding: evidence from the lattice model[END_REF]46 for this SH3, in agreement with what proposed earlier by Shakhnovich and co-workers 26 .

Robustness of two-state folding in Grb2-SH3

It has been proposed that some SH3 domains may retain some residual structure in their denatured state 47,48 and/or populate folding intermediates [20][21][START_REF] Dasgupta | Four-state folding of a SH3 domain: saltinduced modulation of the stabilities of the intermediates and native state[END_REF] . In order to test the robustness of the two-state folding in Grb2-SH3, we resorted to analyse the dependence of the folding parameters as a function of protein stability. In fact, comparing the parameters measured on different site-directed variants represents an efficient test to address the overall folding characteristics of transition and denatured state 49,50 . More specifically, since the dependence of activation and ground states free energies on the denaturant concentration (measured by the m U , m F and m D-N values) are dependent from the changes in accessible surface area between the pertinent state [START_REF] Myers | Denaturant m values and heat capacity changes: relation to changes in accessible surface areas of protein unfolding[END_REF] , an analysis of their dependence may be reveal signatures of shifts of the transition and denatured states along the reaction coordinate, as well as the accumulation of folding intermediates. The Journal of Physical Chemistry domain [20][21][START_REF] Dasgupta | Four-state folding of a SH3 domain: saltinduced modulation of the stabilities of the intermediates and native state[END_REF] , the folding mechanism of this protein is robust and consistent with twostate.

DISCUSSION

The first comparative Φ value analysis on globular proteins was presented in a two papers describing the folding of src and spectrin SH3 16,17 . These studies suggested this protein family to fold via a conserved mechanism characterized by a structurally robust transition state. Furthermore, it was pointed out that the structure of the transition state was primarily stabilized by interactions taking place in the third βhairpin of the protein, representing a polarized folding nucleus. Subsequently, also a Φ value analysis of fyn SH3 domain was reported, further supporting the robustness of the structure of the transition state [START_REF] Northeyk | Hydrophobic core packing in the SH3 domain folding transition state[END_REF] showing that even drastic non conservative mutations caused little structural rearrangements of the transition state 51 . A breakdown of such robustness could be observed in SSo7d, a protein sharing a similar topology with the other SH3 domains while displaying negligible sequence homology. In fact, in this case, a shift in the transition state nucleus from the third to the second β-hairpin was reported 19 .

In the context of previous work on SH3 domains, it is therefore interesting to note how the structure of the transition state of folding of Grb2-SH3 is different from that of src, spectrin and fyn SH3. In fact, Grb2-SH3 displays an extended folding nucleus, which involves the β-sheet comprising the N-and C-termini of the protein together with the first β-hairpin. Since the structural architecture of the folding nucleus of Grb2-SH3 appears to be distinct from that of Sso7d, it appears that this protein family may fold through a multitude of mechanisms comprising distinct regions the protein.

Such pathways may then be selectively stabilized over others by the amino acid sequence, indicating that, whilst the overall features of folding are defined by protein topology, the nature of the interactions stabilising the native state are still critical to influence protein folding mechanisms, In this context, alternative pathways may emerge when the sequence is changed extensively.

A number of studies have shown that proteins may fold with or without folding intermediates, depending on solvent conditions and changes in sequence composition [START_REF] Capaldi | Ultrarapid mixing experiments reveal that Im7 folds via an on-pathway intermediate[END_REF][START_REF] Capaldi | Im7 folding mechanism: misfolding on a path to the native state[END_REF][START_REF] Ferguson | Rapid folding with and without populated intermediates in the homologous four-helix proteins Im7 and Im9[END_REF][START_REF] White | Simulation and experiment conspire to reveal cryptic intermediates and a slide from the nucleation-condensation to framework mechanism of folding[END_REF]15,[START_REF] Gianni | Structural characterization of a misfolded intermediate populated during the folding process of a PDZ domain[END_REF] . Accordingly, whilst the folding of SH3 domains has been classically descripted with a two-state mechanism, Dokholyan and co-workers predicted 52 , by analysing different molecular dynamics simulations, that this class of protein may populate stable intermediates as a consequence of the local stabilization of individual structural elements. This finding was later supported by NMR and by pulse hydrogen exchange mass spectrometry, that revealed that presence of at least one folding intermediate in the case of Fyn 20,21 and PI3K SH3 22 respectively. In both cases, the stabilization of the intermediate appears to arise from the stabilization of non-native hydrophobic interactions, leading to a polarized structure formation upon folding. The analysis of the Bronsted plot of Grb2-SH3 suggests this protein to fold via a transition state with diffused native-like structure. In this case, therefore, the protein seems consistent with a nucleation-condensation mechanism, characterized by a concurrent formation of native secondary and tertiary structure [START_REF] Abkevich | Specific nucleus as the transition state for protein folding: evidence from the lattice model[END_REF]46 . On the light of this finding, it is not surprising to observe that, contrary to the SH3 domains of Fyn and PI3K, Grb2-SH3 seems to conform to two state folding, even when challenged with different site-directed variants, as illustrated by the robustness of the measured m F , m U and m D-N values, which are essentially independent of protein stability.

Taken together, our analysis of the folding pathway of Grb2-SH3 supports a view whereby this class of proteins may fold via alternative pathways, stabilized by different nuclei, that can be selectively balanced by sequence composition. In agreement with previous finding on other protein systems, local stabilization of such alternative nuclei may lead to the accumulation of intermediates, switching two-state to multi-state folding.

ACS Paragon Plus Environment

The Journal of Physical Chemistry Panel B: Bronsted plot. As explained in the Results, the linearity of the Bronsted plot suggests that this protein fold via a native-like diffused, transition state. This finding appears consistent with a nucleation-condensation mechanism [START_REF] Abkevich | Specific nucleus as the transition state for protein folding: evidence from the lattice model[END_REF]46 .

  high Φ value), with the others displaying low values of Φ 44 . As evident from Figure5B, the transition state of Grb2-SH3 clearly displays a linear

Figure 6

 6 depicts the correlation between the m D-N , m U , and m F and the ∆∆G D-N for the different site-directed variants. It is evident that, in the case of Grb2-SH3, no detectable change in m U , m F and m D-N values could be observed for the different variants, spanning a change in protein stability of about 3 kcal mol -1 . This observation suggests that, contrary to what observed in the case of fyn and PI3K SH3 Page 11 of 29 ACS Paragon Plus Environment

Figure 1 :

 1 Figure 1: Cartoon rapresentation of Grb2-SH3 (A), Spectrin-SH3 (B) and Sso7d-SH3

Figure 2 :

 2 Figure 2: Equilibrium denaturation experiment of the Grb2-SH3 domain carried out in

Figure 3 :

 3 Figure 3: Chevron plot of the Grb2-SH3 domain obtained in buffer 50 mM sodium

Figure 4 :

 4 Figure 4: Equilibrium denaturations and chevron plots of Grb2-SH3 and its site

Figure 5 :

 5 Figure 5: Structure of the folding transition state of Grb2-SH3, together with the

Figure 6 :

 6 Figure 6: Correlation between the m D-N (open circle), m U (rhombus), and m F (squares)

  Figure 1

Figure 4 Page

 4 Figure 4

Table 1 .

 1 Kinetic folding parameters of Grb2-SH3 and its site-directed variants.

		k F	k U	m F	m U	m D-N	[urea] 1/2
	MUTANT	(s -1 )	(s -1 )	(kcal/M mol)	(kcal/M mol)	(kcal/M mol)	(M)	φ φ φ φ
	WT	16.0±1.4	0.14±0.02 0.59±0.02 0.14±0.06 0.73±0.01	4.2±1.0
	T1S	13.0±1.3	0.21±0.02 0.59±0.06 0.12±0.12 0.71±0.10	3.1±0.2	0.28±0.26
	Y2A	8.6±0.9	0.69±0.07 0.59±0.10 0.07±0.10 0.73±0.04	1.8±0.1	0.26±0.07
	V3A	4.6±0.6	0.33±0.05 0.75±0.07 0.08±0.09 0.83±0.06	2.5±0.3	0.58±0.11
	A5G	1.7±0.7	0.36±0.04 0.63±0.05 0.12±0.08 0.73±0.06	0.2±0.1	0.70±0.17
	L6A	8.0±0.6	0.67±0.05 0.68±0.04 0.13±0.06 0.81±0.04	2.3±1.1	0.30±0.06
	F7A	6.4±0.64	0.21±0.02 0.66±0.07 0.11±0.10 0.77±0.07	3.2±0.1	0.68±0.15
	F19A	2.6±0.3	1.23±0.05 0.57±0.09 0.07±0.13 0.73±0.09	1.6±0.2	0.45±0.04
	F24A	3.7±0.3	0.70±0.03 0.70±0.04 0.12±0.06 0.81±0.04	1.8±0.2	0.47±0.05
	I25V	6.1±0.5	0.64±0.03 0.70±0.04 0.08±0.06 0.78±0.04	2.4±1.7	0.38±0.06
	H26A	33.0±8.0	0.23±0.05 0.62±0.05 0.15±0.06 0.77±0.03	3.6±1.5	*
	S31A	11.0±1.0	0.42±0.04 0.67±0.03 0.08±0.04 0.75±0.03	2.8±0.4	0.24±0.10
	A39G	16.0±2.0	0.92±0.07 0.71±0.04 0.06±0.06 0.77±0.04	2.4±1.3	-0.03±0.07
	H41A	13.0±1.3	0.56±0.06 0.65±0.07 0.07±0.07 0.72±0.02	3.0±0.1	0.11±0.09
	T44S	14.0±1.0	0.53±0.05 0.69±0.07 0.05±0.07 0.74±0.02	2.8±0.1	0.07±0.09
	Y51A	11.0±0.5	0.43±0.02 0.57±0.02 0.14±0.02 0.71±0.01	2.6±0.4	0.23±0.07
	T53S	11.0±3.0	0.23±0.05 0.53±0.06 0.12±0.09 0.65±0.07	2.7±0.3	0.41±0.37
	A54G	16.2±1.6	0.62±0.06 0.55±0.03 0.12±0.04 0.67±0.02	2.4±1.1	-0.05±0.09

Table 1 :

 1 The

					Page 22 of 29
	V55A	15.2±1.5	0.45±0.04 0.47±0.05 0.23±0.09 0.70±0.07	3.4±0.3	-0.01±0.11
		mutants F9A, L17A, V27A, F47A and V52A expressed poorly and
		could not be characterized.
		*This mutant shows ܩ߂߂ ିே < 0.4 kcal mol -1 , preventing reliable calculation of
		the Φ-value 36 .

ACKNOWLEDGEMENTS

Work partly supported by grants from the Italian Ministero dell'Istruzione dell'Università e della Ricerca (Progetto di Interesse 'Invecchiamento' to S.G.), Sapienza University of Rome (C26A155S48, B52F16003410005 and RP11715C34AEAC9B to S.G), the Associazione Italiana per la Ricerca sul Cancro (Individual Grant -MFAG 2016, 18701 to S.G.). FT is a recipient of a PhD fellowship from the Italo-French University.

ACS Paragon Plus Environment

The Journal of Physical Chemistry