Peter Dürrschmidt
email: peter.duerrschmidt@wacker.com

Johanna Mansfeld

Renate Ulbrich-Hofmann
email: renate.ulbrich-hofmann@biochemtech.uni-halle.de

Refolding of the non-specific neutral protease from Bacillus stearothermophilus proceeds via an autoproteolytically sensitive intermediate

Keywords: Autoproteolysis, disulfide bond, folding intermediate, neutral protease, reactivation, refolding hydrochloride, k obs , observed rate constant, pWT, pseudo-wild type, TCA, trichloroacetic acid

A very thermostable variant of the thermolysin-like protease from Bacillus stearothermophilus (G8C/N60C) was previously created by introduction of a disulfide bond into the cysteine-free pseudo-wild type variant (pWT) and thus fixing the unfolding region 56-69. In the present paper, we show that G8C/N60C and pWT can be reactivated from the completely unfolded states, accessible at ≥ 7.5 M guanidine hydrochloride, and analyze the kinetics of folding, autoproteolytic degradation and aggregation. From changes in the fluorescence spectra with time of renaturation, it can be concluded that a folding intermediate with native-like structure, but which is still inactive and sensitive to autoproteolysis, is rapidly formed after renaturation has initiated. The critical region 56-69 of pWT is involved in the autoproteolytic sensitivity of the intermediate as we conclude from the differences in the chevron plots of the first order rate constants of reactivation and the fragmentation patterns in SDS-PAGE of pWT and G8C/N60C.

Please cite this article as: Peter Dürrschmidt, Johanna Mansfeld, Renate Ulbrich-Hofmann, Refolding of the non-specific neutral protease from Bacillus stearothermophilus proceeds via an autoproteolytically sensitive intermediate, Biophysical Chemistry (2010), doi: 10.1016/j.bpc.2010.01.001 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

A C C E P T E D

M A N U S C R I P T

Introduction

Proteolytic enzymes differ from other enzymes by their capability of autoproteolysis. Generally, they are not susceptible to autoproteolysis under native conditions but become susceptible to autoproteolysis when their native structure is impaired. For this reason proteases are irreversibly inactivated if they are exposed to denaturing influences such as high temperatures or chemical denaturants. During protein biosynthesis, proteases are liable to be degraded as long as their native folded structure has not been attained. In nature several strategies have evolved to protect proteases from autoproteolytic degradation. Many proteases are produced as inactive proenzymes, which attain their active conformation only after a specific proteolytic cleavage and release of the propeptide [START_REF] Baker | The role of pro regions in protein folding[END_REF][START_REF] Bryan | Prodomains and protein folding catalysts[END_REF]. Thus, most digestive proteases as well as bacterial extracellular proteases are secreted as inactive zymogens, e.g. elastase as proelastase or subtilisin as prosubtilisin.

Alternatively, many proteases have specific inhibitors that prevent activity. Thus, serpins are a well-known class of structurally similar compounds with particular importance as inhibitors of serine proteinases [START_REF] Silverman | The serpins are an expanding superfamily of structurally similar but functionally diverse proteins. Evolution, mechanism of inhibition, novel functions, and a revised nomenclature[END_REF].

Due to phenomenon of autoproteolysis, it is generally not possible to follow the unfolding of proteases without concomitant protein degradation [START_REF] Corbett | Domain unfolding and the stability of thermolysin in guanidine hydrochloride[END_REF][START_REF] Sättler | Temperature-gradient gel electrophoresis for analysis and screening of thermostable proteases[END_REF][START_REF] Szeltner | Conformational stability and catalytic activity of HIV-1 protease are both enhanced at high salt concentration[END_REF]. Similarly, folding (refolding) processes cannot be observed by biophysical methods. Therefore, the thermodynamic stability (ΔG°) of proteases is also mostly inaccessible to biochemical investigation.

Many proteases also need a propeptide for correct folding [START_REF] Eder | Pro-sequence-assisted protein folding[END_REF]. From extensive studies on the αlytic protease from Lysobacter enzymogenes, it was concluded that mature proteases are thermodynamically unstable but kinetically trapped by a high energetic barrier to unfolding, whereas folding is catalyzed by the propeptide [START_REF] Jaswal | Energetic landscape of alpha-lytic protease optimizes longevity through kinetic stability[END_REF]. These constraints also seemed to be valid for the unspecific proteases of the thermolysin type [START_REF] Marie-Claire | Intramolecular processing of prothermolysin[END_REF]. Recently, however, we could show that the

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

ACCEPTED MANUSCRIPT 4 neutral protease from Bacillus stearothermophilus can be expressed in E. coli with an N-terminal His tag without the help of the propeptide [START_REF] Mansfeld | The propeptide is not required to produce catalytically active neutral protease from Bacillus stearothermophilus[END_REF]. The enzyme accumulated as inclusion bodies and, unexpectedly, it could be renatured in high yields after solubilisation of the protein aggregates in guanidine hydrochloride (GdnHCl). Moreover, we could show in a study comparing this protease to its extremely stable disulfide-engineered variant G8C/N60C in GdnHCl-induced denaturation [START_REF] Dürrschmidt | Differentiation between conformational and autoproteolytic stability of the neutral protease from Bacillus stearothermophilus containing an engineered disulfide bond[END_REF][START_REF] Dürrschmidt | An engineered disulfide bridge mimics the effect of calcium to protect neutral protease against local unfolding[END_REF] that these enzymes can be unfolded without significant autoproteolysis in the presence of very high concentrations of GdnHCl. These findings motivated us to investigate the reactivation of this protease and its thermostable variant G8C/N60C, although there are only a few related studies in the literature, probably because of the low renaturation yields and the complex interplay of folding and irreversible processes, such as autoproteolysis and aggregation. For example, however, almost complete restoration of enzymatic activity was reported for subtilisin BPN´ in the presence of organic salts and a temporary inhibitor [START_REF] Matsubara | Achievement of renaturation of subtilisin BPN' by a novel procedure using organic salts and a digestible mutant of Streptomyces subtilisin inhibitor[END_REF]. Moreover, high yields of renaturation could also be achieved with immobilized thermolysin in the presence of salts of organic acids [START_REF] Nohara | Media selection for refolding of thermolysin by use of immobilized preparation[END_REF].

The neutral protease from B. stearothermophilus as one of the thermolysin-like members of the M4 family (with 85 % sequence identity to thermolysin) consists of 319 amino acid residues and belongs to the large family of the zinc-dependent metalloendopeptidases. Because of its phylogenetic relationship with the human matrix metalloproteinases and other physiologically important mammalian proteases such as endothelin-and angiotensin-converting enzymes, neprilysin or neurolysin, it has also importance as a model system for these proteases. Based on the resolved X-ray structures of its homologous bacterial relatives from B. thermoproteolyticus (known as thermolysin) and B. cereus, a 3D model of this enzyme is available [START_REF] Vriend | Prediction and analysis of structure, stability and unfolding of thermolysin-like proteases[END_REF]. Fig. 1 shows this model for the disulfide-containing variant G8C/N60C. The enzyme is organized in two

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT 5
domains with the active site between, containing the catalytic zinc ion and four binding sites for calcium ions with stabilizing function.

Figure 1

Extensive mutational studies had revealed that mutations in the loop region 56-69 (Fig. 1) are critical for the thermal stability of the protease whereas mutations in other regions have only marginal effects [START_REF] Eijsink | Structural determinants of the stability of thermolysin-like proteinases[END_REF]. In the engineered enzyme variant G8C/N60C [START_REF] Mansfeld | Extreme stabilization of a thermolysin-like protease by an engineered disulfide bond[END_REF] a disulfide bond connects this region with the underlying β-sheet structure (Fig. 1), which causes a dramatic stabilization against thermal inactivation. The variant was produced by site-directed mutagenesis starting from the cysteine-free variant of the protease (C288L), which does not significantly differ from the wild-type enzyme either in activity or in stability [START_REF] Dürrschmidt | Differentiation between conformational and autoproteolytic stability of the neutral protease from Bacillus stearothermophilus containing an engineered disulfide bond[END_REF][START_REF] Eijsink | The effect of cavity-filling mutations on the thermostability of Bacillus stearothermophilus neutral protease[END_REF]. This enzyme, called pseudowild type (pWT), is used as reference enzyme in studies on G8C/N60C. The introduction of the disulfide bond produced a shift in the half-life at 92.5°C from < 0.3 min to 35.9 min and increased the temperature after which half of the initial activity is lost within 30 min (T 50 ) by 16.7°C [START_REF] Mansfeld | Extreme stabilization of a thermolysin-like protease by an engineered disulfide bond[END_REF].

In GdnHCl-induced denaturation of pWT and G8C/N60C, unfolding and autoproteolysis could be separated because the proteolytic activity is strongly inhibited by GdnHCl even in the presence of low denaturant concentrations [START_REF] Dürrschmidt | Differentiation between conformational and autoproteolytic stability of the neutral protease from Bacillus stearothermophilus containing an engineered disulfide bond[END_REF][START_REF] Dürrschmidt | An engineered disulfide bridge mimics the effect of calcium to protect neutral protease against local unfolding[END_REF]. The results suggested that the disulfide bond stabilizes the protease against autoproteolysis more than against global unfolding [START_REF] Dürrschmidt | Differentiation between conformational and autoproteolytic stability of the neutral protease from Bacillus stearothermophilus containing an engineered disulfide bond[END_REF]. Obviously, unfolding of pWT proceeds via local unfolding of the critical region which is sensitive to autoproteolysis. The disulfide bond fixes this region and protects therefore against autoproteolysis. A similar effect is obtained at high Ca 2+ concentrations (100 mM) by saturation of the low affinity calcium binding site III being located in the critical loop region between amino acids 56-69 [START_REF] Dürrschmidt | An engineered disulfide bridge mimics the effect of calcium to protect neutral protease against local unfolding[END_REF].

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT 6
In the present paper we show that pWT and G8C/N60C can undergo substantial reactivation starting from the completely GdnHCl-unfolded states. Comparison of the reactivation kinetics with the kinetics of refolding by fluorescence spectroscopy and autoproteolytic degradation by SDS-PAGE reveals that refolding proceeds via a native-like intermediate which is susceptible to autoproteolysis. The combination of the kinetic refolding data with the results of unfolding [START_REF] Dürrschmidt | An engineered disulfide bridge mimics the effect of calcium to protect neutral protease against local unfolding[END_REF] allows an estimate of the standard Gibbs free energies of the reversible transition from the native state to the intermediate and also allows a comparison of their kinetic and thermodynamic stabilities, which has been difficult to achieve with proteases until now.

Materials and Methods

Chemicals

Casein was purchased from Merck (Darmstadt, Germany), 2-aminobenzoyl-Ala-Gly-Leu-Ala-4nitrobenzylamide (Abz-AGLA-Nba) from Bachem (Heidelberg, Germany), GdnHCl from ICN Biomedicals GmbH (Eschwege, Germany), sodium deoxycholate from Sigma-Aldrich Chemie GmbH (Deisenhofen, Germany), and tris(hydroxymethyl)aminomethane (Tris) from Amersham Biosciences (Uppsala, Sweden). All other reagents were of the purest quality available.

Enzymes

pWT and G8C/N60C were produced as described previously [START_REF] Mansfeld | Extreme stabilization of a thermolysin-like protease by an engineered disulfide bond[END_REF][START_REF] Eijsink | The effect of cavity-filling mutations on the thermostability of Bacillus stearothermophilus neutral protease[END_REF]. Immediately before use, the enzyme solution was dialyzed against 50 mM Tris/HCl buffer, pH 7.5, containing 5 mM CaCl 2 , incubated at 65°C for 8 min and dialyzed again. After this procedure the enzyme solutions were homogeneous as judged by SDS-PAGE. The absence of low molecular weight fragments was checked by size-exclusion chromatography using a Superdex 75 (10/30) column (Pharmacia,

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT 7
Uppsala, Sweden) with 50 mM Tris/HCl buffer, pH 7.5, containing 5 mM CaCl 2 and 20% (v/v) isopropanol as eluent.

For the preparation of solutions containing denatured enzymes at protein concentrations of 5 to 20 µg mL -1 (for reactivation experiments) the native enzymes were adjusted to a protein concentration of 0.1 to 0.4 mg mL -1 , diluted 1:20 into 50 mM Tris/HCl buffer, pH 7.5, containing 5 mM CaCl 2 and 8.4 M GdnHCl and incubated for at least 3 hours. For the preparation of solutions containing denatured enzymes at protein concentrations higher than 20 µg mL -1 (for fluorescence spectroscopy, SDS-PAGE analysis and aggregation experiments) the native enzymes were first precipitated by TCA in the presence of sodium deoxycholate, washed with acetone and ethanol, and dried under vacuum as described in [START_REF] Dürrschmidt | Differentiation between conformational and autoproteolytic stability of the neutral protease from Bacillus stearothermophilus containing an engineered disulfide bond[END_REF]. Then they were dissolved in 50 mM Tris/HCl buffer, pH 7.5, containing 5 mM CaCl 2 and 8.4 M GdnHCl to the protein concentration needed for the subsequent experiments (typically 5 -10 mg mL -1 ).

Protein concentration

The protein concentration was determined using the BCA protein assay (Pierce, Germany) with bovine serum albumin as standard.

Activity toward casein

The activity was determined at 37°C by 30 min incubation in 50 mM Tris/HCl buffer, pH 7.5, containing 5 mM CaCl 2 and 0.5% (w/v) casein as described previously [START_REF] Fujii | Molecular cloning of a thermostable neutral protease gene from Bacillus stearothermophilus in a vector plasmid and its expression in Bacillus stearothermophilus and Bacillus subtilis[END_REF][START_REF] Kunitz | Crystalline soybean trypsin inhibitor. II. General properties[END_REF]. To determine the reactivation of pWT and G8C/N60C, the denatured enzymes (in GdnHCl-containing solution)

were diluted into 50 mM Tris/HCl buffer, pH 7.5, containing 5 mM CaCl 2 . Samples were taken from this renaturation mixture at different times and assayed for activity.

A C C E P T E D

M A N U S C R I P T ACCEPTED MANUSCRIPT 8 2.

Activity toward Abz-AGLA-Nba

The activity was measured at a substrate concentration of 20 µM in 50 mM Tris/HCl buffer, pH 7.5, containing 5 mM CaCl 2 at 25°C. The increase in fluorescence emission at 415 nm after excitation at 340 nm was recorded over at least 10 min [START_REF] Nishino | Pseudomonas aeruginosa elastase. Development of a new substrate, inhibitors, and an affinity ligand[END_REF]. To determine the reactivation of pWT and G8C/N60C with half life > 5 min, the denatured enzymes (in GdnHCl containing solution) were diluted into 50 mM Tris/HCl buffer, pH 7.5, containing 5 mM CaCl 2 . After several periods of time samples were taken from this renaturation mixture and assayed for activity. The progress curves (change of activity vs. renaturation time) were fitted to a single exponential function yielding the rate constants k obs .

Figure 2

For monitoring fast reactivation kinetics (half life < 5 min) the substrate was directly added to the renaturation mixture and the change in the fluorescence signal was followed as function of time as exemplified in Fig. 2. The linear section of the curve representing the final recovered activity was fitted using the function y = m • x + n, where y is the fluorescence signal and x is time.

Subtraction of the corresponding y-values from the fluorescence data in the non-linear initial phase yielded exponential curves (Fig. 2, insert) from which the rate constants of reactivation (k obs ) were calculated.

Autoproteolysis

SDS-PAGE was used to follow autoproteolysis of the enzymes at a protein concentration ≥ 20 µg mL -1 . To concentrate the protein solution and to remove GdnHCl, the samples taken from the renaturation mixture were treated with sodium deoxycholate [START_REF] Dürrschmidt | Differentiation between conformational and autoproteolytic stability of the neutral protease from Bacillus stearothermophilus containing an engineered disulfide bond[END_REF]. The precipitates were dried

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

under vacuum, dissolved in SDS-sample buffer and separated on 15 % polyacrylamide gels according to Laemmli [START_REF] Laemmli | Cleavage of structural proteins during the assembly of the head of bacteriophage T4[END_REF]. The proteins were stained with Coomassie Brilliant Blue G-250 for quantitative evaluation of the bands of intact enzyme by scanning with a CD 60 densitometer (Desaga, Darmstadt, Germany), or with silver [START_REF] Nesterenko | A simple modification of Blum's silver stain method allows for 30 minute detection of proteins in polyacrylamide gels[END_REF] for evaluation of the fragmentation pattern.

Fluorescence measurements

Fluorescence spectroscopy was carried out in 1 cm quartz cuvettes at 25°C with a FluoroMax-2 spectrofluorometer (Yvon-Spex, Grasbrunn, Germany) with excitation at 278 nm. The slit width was 5 nm for both excitation and emission. Protein spectra were corrected by blank spectra.

Results

Despite the strong tendency of being degraded by autoproteolysis, the variants of the neutral protease from B. stearothermophilus pWT and G8C/N60C can be obtained in the intact, unfolded state at high concentrations (≥ 7.5 M) of GdnHCl [START_REF] Dürrschmidt | Differentiation between conformational and autoproteolytic stability of the neutral protease from Bacillus stearothermophilus containing an engineered disulfide bond[END_REF][START_REF] Dürrschmidt | An engineered disulfide bridge mimics the effect of calcium to protect neutral protease against local unfolding[END_REF]. This finding enabled us to investigate refolding of these proteases. Other denaturants such as urea or high temperature were not appropriate for observing refolding because they were not able to yield completely unfolded protein without autoproteolytic fragmentation. The transfer of the denatured enzymes ([GdnHCl] ≥ 7.5 M) to native-like conditions ([GdnHCl] ≤ 1.0 M) was investigated by following enzyme activity, intrinsic fluorescence, autoproteolysis, and aggregation. Native pWT and G8C/N60C are stable toward autoproteolysis and keep their structural integrity for several weeks [START_REF] Dürrschmidt | An engineered disulfide bridge mimics the effect of calcium to protect neutral protease against local unfolding[END_REF].

In all the studies, Ca 2+ ions (5 mM) were added to stabilize the protein. Zn 2+ ions, although essential cofactors and exhibiting extremely high affinity, were without effect on the stability and were not added exogenously. For technical reasons, different protein concentrations had to be chosen in the various experiments, which must be considered in the interpretation of the results. GdnHCl. As concluded from light-scattering studies and SDS-PAGE analyses at higher enzyme concentrations (see below), aggregation and autoproteolysis strongly compete with reactivation and are probably the reason for incomplete reactivation even at low protein concentrations.

Reactivation

Spectroscopic characterization of refolding

For spectroscopic measurements of the refolding process higher enzyme concentrations (≥ 5 µg mL -1 ) were necessary than for the determination of reactivation (0.1 µg mL -1 ). To study refolding, fluorescence spectra of pWT and G8C/N60C were recorded in the protein concentration range of 5 -30 µg mL -1 . Immediately after transferring the unfolded enzymes to native-like conditions (≤ 1 M GdnHCl) native-like fluorescence spectra were observed (Fig. 5).

The characteristic separation of the tyrosine and tryptophan fluorescence, which is characteristic for the denatured state of the enzymes (Fig. 5) and has been studied in detail previously [START_REF] Dürrschmidt | Differentiation between conformational and autoproteolytic stability of the neutral protease from Bacillus stearothermophilus containing an engineered disulfide bond[END_REF], was abolished.

Figure 5

These spectral changes occur within 10 -15 seconds, which is a time-scale where the reactivation

is not yet completed. Therefore, an inactive intermediate with an at least partly folded structure, allowing resonance energy transfer from the tyrosine to the tryptophan residues, must be formed within the first seconds after the initiation of refolding. The fluorescence spectra of the refolded enzymes showed lower intensities than those of the native enzymes under the same conditions which may be attributed to concomitant aggregation (section 3.4). The maximum fluorescence of refolded G8C/N60C was somewhat larger than that of pWT, reflecting a higher refolding yield as

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

13 was also found in the reactivation study (section 3.1). Examination of the shape of the spectra reveals, however, that there were only minor deviations between the two enzyme variants.

Figure 6

With increasing time, the native-like spectra changed and the fluorescence signals of tyrosine and tryptophan split again (Fig. 6A andB). This slow process becomes measurable at a higher residual GdnHCl concentration in the renaturation assay (e.g. 0.64 M) as demonstrated in Fig. 6.

A kinetic evaluation of the fluorescence signals revealed first order kinetics (Fig. 6C andD) with rate constants in the same range for pWT (2.6 • 10 -4 s -1 ) and G8C/N60C (3.0 • 10 -4 s -1 ).

Autoproteolysis

The decrease in fluorescence intensity as described above was accompanied by autoproteolysis as revealed by SDS-PAGE analysis in the protein concentration range of 5 -50 µg mL -1 . Obviously, the intermediate state occurring immediately after the transfer of denatured enzyme from high to low GdnHCl concentrations (section 3.2.) was susceptible to autoproteolysis. In contrast, both native pWT and G8C/N60C were very stable and did not undergo autoproteolysis. Interestingly, distinct fragments accumulated during this autoproteolysis. The optimal conditions to visualize fragments were found to be at a residual GdnHCl concentration of 0.42 M and a final protein concentration of 50 µg mL -1 (Fig. 7). There were some slight differences in the fragmentation pattern between pWT and G8C/N60C (Fig. 7A andB), which was probably caused by differences in the accessibility of the region around the amino acids 8 and/or 60. Under the conditions shown in Fig. 6, the kinetics of autoproteolysis determined from the decreasing band intensities of the intact proteases in SDS-PAGE analysis after Coomassie staining could be described by a first-order reaction (Fig. 6C andD). The resulting rate constants (1.6 • 10 -4 s -1 for pWT and 1.9 • 10 -4 s -1 for G8C/N60C) were very similar to the rate constants as determined from the decrease of the fluorescence signals (section 3.2). In both cases, the degradation of the autoproteolytically sensitive folding intermediate was observed.

Aggregation

Aggregations occurring during renaturation of 5 -60 µg mL -1 pWT and G8C/N60C were detected by light scattering at 360 nm (Fig. 8). In relation to the light-scattering of solutions containing the native enzymes precipitated by trifluoroacetic acid, the aggregation of pWT and G8C/N60C amounted to 12 and 22 %, respectively, at a final GdnHCl concentration of < 100 mM.

Figure 8

Solutions of native pWT and G8C/N60C did not show any light-scattering. Analytical ultracentrifugation confirmed that the enzymes remained monomeric up to protein concentrations of at least 0.5 mg mL -1 (data not shown).

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

Discussion

Reactivation competes with autoproteolysis and aggregation

In contrast to the general view on the folding of proteases, pWT and G8C/N60C can be renatured (Fig. 3), after denaturation in highly concentrated GdnHCl solutions (where no significant autoproteolysis occurs [START_REF] Dürrschmidt | Differentiation between conformational and autoproteolytic stability of the neutral protease from Bacillus stearothermophilus containing an engineered disulfide bond[END_REF]) followed by subsequent dilution with buffer to native-like conditions.

The reactivation competes with autoproteolysis and aggregation processes, which can be suppressed at low protein concentrations. At protease concentrations of 0.1 µg mL -1 maximal reactivation yields of 50 % (pWT) and 70 % (G8C/N60C) could be obtained, indicating that the proteases can attain the correctly folded tertiary structure without the help of folding helpers such as the natural propeptide or chaperones.

Based on the first-order kinetics of the reactivation reaction (Fig. 4), it can be concluded that the observed rate constants reflect the folding step in which the enzymes obtain their native state. In autoproteolysis and aggregation reactions, where two or more molecules are involved, rate constants are of second or higher order.

The reactivation rate constants give linear branches in the chevron plot (Fig. 4). From fluorescence experiments we learned that refolding proceeds via a rapidly formed intermediate (sections 3.2. and 4.2.). Therefore, the first-order reaction observed in activity measurements is suggested to reflect the second, slower step in which the active state is formed from the intermediate. The derived Gibbs free activation energies of this folding step (ΔG # ) and the slopes of the straight lines (m) in the plots of Fig. 4 are summarized in Table 1. The smaller slope of the refolding branch of G8C/N60C compared to that of pWT reflects the more rigid structure of the disulfide variant. Moreover, the differences in the slopes of the straight lines indicate that the disulfide bond in G8C/N60C influences the observed refolding.

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

Refolding proceeds via intermediates sensitive to autoproteolysis

As fluorescence measurements in the protein concentration range of 5 -30 µg mL -1 showed, refolding of pWT and G8C/N60C proceeds via intermediates with native-like structure. Upon dilution of the enzymes previously unfolded using 8 M GdnHCl, these intermediates are formed within 10 to 15 seconds and are spectroscopically very similar to the native states (Fig. 5). The finding that the total fluorescence intensities of the corresponding amounts of native enzymes are not reached can be attributed to accompanying aggregation of the proteases, as indicated by results of light-scattering experiments (Fig. 8). The folding intermediates are slowly degraded as concluded from the changes in the fluorescence spectra with time (Fig. 6) and the generation of peptide fragments as monitored by SDS-PAGE (Fig. 7). Obviously, traces of the proteases, that have already attained their native, autoproteolytically-resistant state, are sufficient to attack and degrade the remaining molecules in the intermediate state. Only at very low protein concentrations as used in the reactivation experiments (0.1 -5 µg mL -1 ), does a major percentage of the molecules have the chance to fold to the native structure.

The occurrence of peptides originating from the degradation of the intermediates during refolding as visualized by SDS-PAGE (Fig. 7), is in contrast to the unfolding experiments where no conditions could be found in which fragment accumulation due to autoproteolysis [START_REF] Dürrschmidt | An engineered disulfide bridge mimics the effect of calcium to protect neutral protease against local unfolding[END_REF] occurred.

It is likely that fragments can accumulate during refolding because the recovered activity is very low compared to the conditions during unfolding where the produced susceptible molecules are rapidly degraded by the large amount of still active species. Differences in the fragment patterns

A C C E P T E D M A N U S C R I P T
of pWT and G8C/N60C (Fig. 7) indicate that the folding intermediates of the two enzymes are differentially accessible to autoproteolysis.

For GdnHCl-induced denaturation it has been shown recently [START_REF] Dürrschmidt | Differentiation between conformational and autoproteolytic stability of the neutral protease from Bacillus stearothermophilus containing an engineered disulfide bond[END_REF][START_REF] Dürrschmidt | An engineered disulfide bridge mimics the effect of calcium to protect neutral protease against local unfolding[END_REF] that unfolding of pWT also proceeds via a native-like intermediate, which is sensitive to autoproteolysis and that autoproteolysis apparently promotes unfolding. In contrast, during the denaturation of G8C/N60C this intermediate does not seem to play any role in autoproteolysis and unfolding. The present results suggest that the same intermediates are formed during folding and unfolding. For pWT and G8C/N60C, these intermediates differ in their sensitivities to degradation. In the intermediate of pWT the region 56-69 is highly flexible and prone to autoproteolysis, whereas in G8C/N60C this region is fixed and insensitive to autodegradation. While we assumed in our earlier studies of unfolding of G8C/N60C that autoproteolysis starts from the completely unfolded enzyme [START_REF] Dürrschmidt | An engineered disulfide bridge mimics the effect of calcium to protect neutral protease against local unfolding[END_REF],

the present results of folding show that autoproteolysis also starts from an intermediate with native-like tertiary structure (Fig. 6B). Therefore, we assume now that this intermediate is also formed but not populated during unfolding of G8C/N60C, and instead only accumulates during refolding.

Combination of the refolding and unfolding data allows the thermodynamic characterization of the intermediate -native state equilibrium

Unfortunately, no conditions could be found to characterize the transition of the intermediate to the native state in equilibrium. With some caution, however, the present and previous kinetic data can be used to estimate the Gibbs free energy, ΔG 0 , of this equilibrium. Assuming that the firstorder rate constants of reactivation reflect the folding from the intermediate to the native state, these data can be combined with the corresponding data of unfolding [START_REF] Dürrschmidt | An engineered disulfide bridge mimics the effect of calcium to protect neutral protease against local unfolding[END_REF], which have been inserted in Fig. 4 and Table 1. Disregarding the incomplete reactivation and making the

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

assumption that the intermediate is the same during folding and unfolding, this allows an energetic characterization of the intermediate-native state equilibrium to be undertaken (Table 1).

The midpoint of the transition [GdnHCl] 0.5 obtained from the intersection of the corresponding refolding and unfolding branches is higher by 1.5 M for G8C/N60C than for pWT. Using

Eyring´s equation [START_REF] Tanford | Protein denaturation. C. Theoretical models for the mechanism of denaturation[END_REF][START_REF] Tanford | Kinetics of unfolding and refolding of proteins. 3. Results for lysozyme[END_REF] the Gibbs free energy of activation under native conditions ΔG # can be obtained from the rate constants extrapolated to 0 M GdnHCl. Comparison of ΔG # with the Gibbs free energy ΔG 0 for both enzyme variants shows that the high stability of these proteins mainly is due to a high kinetic barrier. From the difference of the ΔG 0 values it can be concluded that the contribution of the disulfide bond to the stabilization is 8.1 kJ mol -1 .

Concluding remarks

For the first time, the interplay of folding and autodegradation processes of an unspecific protease

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

The denatured enzymes (20 µg mL -1 in 50 mM Tris/HCl buffer, pH 7.5, containing 5 mM CaCl 2 and 8 M GdnHCl) were diluted 1:20 into the same buffer but without GdnHCl. After several periods of time samples were removed and assayed for activity using the casein assay as described in section 2.4. Each value is the mean of three independent measurements. The error bars represent the standard deviation. The regained activity was related to the activity of the native enzyme at the same GdnHCl concentration, which was used as control. The insert shows the data in combination with the rate constants of unfolding for pWT (▲) and G8C/N60C (Δ) as a function of the GdnHCl concentration (data taken from [START_REF] Dürrschmidt | An engineered disulfide bridge mimics the effect of calcium to protect neutral protease against local unfolding[END_REF]). The parameters were calculated from the data of Fig. 4 according to Tanford [START_REF] Tanford | Protein denaturation. C. Theoretical models for the mechanism of denaturation[END_REF][START_REF] Tanford | Kinetics of unfolding and refolding of proteins. 3. Results for lysozyme[END_REF] and Eyring The data marked by a star were taken from [START_REF] Dürrschmidt | An engineered disulfide bridge mimics the effect of calcium to protect neutral protease against local unfolding[END_REF].

10 %

 10 Kinetic measurements were started by dilution of the denatured enzymes (in GdnHCl-containing solutions) into 50 mM Tris/HCl buffer, pH 7.5, containing 5 mM CaCl 2 . If necessary, GdnHCl was added to adjust the final GdnHCl concentration. The progress curves (change of signal vs. renaturation time) were fitted to a single exponential function yielding the rate constants k obs .2.8 Light scatteringLight scattering to monitor aggregation was measured in 1 cm quartz cuvettes at 25°C with a FluoroMax-2 spectrofluorometer (Yvon-Spex, Grasbrunn, Germany) and with excitation and emission at 360 nm. The slit width was 1 nm for both excitation and emission. First, the baseline was recorded by measuring the renaturation buffer for 5 minutes with gentle stirring. Then the denatured enzyme (in GdnHCl-containing solutions) was added and the change in the signal was recorded for 15 minutes while gently stirring the solution. The amplitude of the increase in the signal was used for quantifying the extent of aggregation. The calibration of the signal was performed by measuring the light-scattering of the protein completely precipitated by TCA (2.5A C C E P T E D M A N U S C R I P TACCEPTED MANUSCRIPT , w/w) and setting the amplitude of this signal as 100 % aggregation. The completeness of precipitation was checked by SDS-PAGE analysis.

Figure 3 A 11 Fig. 3

 3113 Figure 3

Figure 4

 4 Figure 4

Figure 7 A

 7 Figure 7

21 [

 21 have been analyzed quantitatively. The results of refolding show that in a first very fast folding step a native-like structure arises. At low protein concentrations (0.1 -5 µg mL -1 ) this intermediate structure can slowly fold to the native active protease, whereas it is degraded at higher protein concentrations by autoproteolysis. Moreover, this study is the first attempt to obtain thermodynamic parameters for the stability of a neutral protease and the influence of an engineered disulfide bond.A C C E P T ED 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 ACCEPTED MANUSCRIPT 26] R. Lumry and H. Eyring, Conformation changes of proteins, J. Phys. Chem. 58 (1954) 110-120.

Figure 1 .

 1 Figure 1. Model of the three-dimensional structure of the neutral protease of B.

Figure 2 .

 2 Figure 2. Progress curve of pWT reactivation.

Figure 3 .

 3 Figure 3. Reactivation of pWT ( ) and G8C/N60C ( ) followed by the casein assay.

Figure 4 .

 4 Figure 4. Rate constants of reactivation of pWT ( ) and G8C/N60C ( ) as a function of the final

Figure 5 .Figure 6 .

 56 Figure 5. Fluorescence spectra of native, unfolded and refolded pWT (A) and G8C/N60C (B).The enzymes were incubated at 5 µg mL -1 in 50 mM Tris/HCl buffer, pH 7.5, containing 5 mM

Figure 7 .and 8 . 4 M

 784 Figure 7. Fragment patterns monitored by SDS-PAGE during refolding of pWT (A) and

Figure 8 . 26 Table 1 :

 8261 Figure 8. Aggregation during refolding of pWT ( ) and G8C/N60C ( ).

[ 26 ]

 26 using the linear extrapolation model ln k obs = m • [GdnHCl] + n, where m represents the slope (∂ln k obs /∂[GdnHCl]) of the folding or unfolding branch and n is the logarithm of the observed rate constant extrapolated to 0 M GdnHCl (ln k 0 ). k obs is given in s -1 . [GdnHCl] 0.5 is the GdnHCl concentration at the point of intersection of the refolding and unfolding branches. ∆G # stands for the Gibbs free energy of activation of unfolding or refolding in the absence of GdnHCl. ∆G 0 was calculated from the difference of the ∆G # values of unfolding and refolding and represents the Gibbs free energy in the absence of GdnHCl.

1 A

 1 Figure 1

Figure 2 A

 2 Figure 2

Table 1

 1

Acknowledgment

The authors thank G. Vriend (Radboud University Nijmegen Medical Centre, Nijmegen, The Netherlands) for preparation of Figure 1, V. G. Eijsink (University of As, Norway) for the plasmid pGE501 encoding the gene of the pseudo-wild type sequence of the neutral protease from B. stearothermophilus, H. Lilie (Martin-Luther-University Halle-Wittenberg, Halle, Germany) for the assistance in ultracentrifugation, R. Golbik Halle-Wittenberg, Halle, Germany) for stimulating discussions, and G. Sawers (Martin-Luther-University Halle-Wittenberg, Halle, Germany) for reading the manuscript. The financial support by the Kultusministerium des Landes Sachsen-Anhalt (Magdeburg, Germany) and the Deutsche Forschungsgemeinschaft (Bonn, Germany) is gratefully acknowledged.


C Fagiano
email: c.fagiano@tudelft.nl

M M Abdalla

C Kassapoglou

Z Gürdal

Interlaminar Stress Recovery for Three-Dimensional Finite Elements

Keywords: A. Layered structures, B. Mechanical properties, C. FEA, Interlaminar stresses

An accurate evaluation of interlaminar stresses in multilayer composite laminates is crucial for the correct prediction of the onset of delamination. In general, three dimensional finite element models are required for acceptable accuracy, especially near free edges and stress concentrations.

Interlaminar stresses are continuous both across and along layer interfaces. Nonetheless, the continuity of interlaminar stresses is difficult to enforce in C 0 interpolated elements. Nodal values of the stresses are usually retrieved using extrapolation techniques from super-convergent points, if known, or Gauss points inside the element. Stress fields within an element can be deduced using either constitutive relations or variationally consistent procedures. In either case, spurious oscillations in stress fields may be encountered leading to a reduced accuracy of the recovered stresses at nodes. In this paper, an efficient interlaminar stress recovery procedure for three-dimensional finite element formulations is presented. The proposed procedure does not rely on extrapolation techniques from super-convergent or integration points. Interlaminar stress values are retrieved directly at nodes and stress continuity at the inter-element boundary is automatically satisfied.

Several benchmark problems were analysed. Comparisons with finite element software and available solutions in the literature confirmed the accuracy of the procedure. Accurate interlaminar stresses were obtained using coarser meshes compared to customary recovery procedures.

Introduction

An accurate evaluation of interlaminar stresses in multilayer composite laminates is crucial for the correct prediction of delamination initiation. Due to its importance, this subject has been the focus of numerous researchers. Several levels of modelling can be used for the evaluation of interlaminar stresses. The simplest approach is based on using Classical Laminated Plate Theory (CLPT). While, CLPT assumes a state of plane stress within each lamina of a multi-directional laminate, interlaminar stresses are most commonly obtained by integrating the three-dimensional equilibrium equations of elasticity through the thickness [START_REF] Reddy | A review of literature on finite element modelling of laminated plates[END_REF][START_REF] Kant | Estimation of transverse/interlaminar stresses in laminated composites -a selective review and survey of current developments[END_REF]. Nonetheless, three dimensional stress states cannot, in general, be accurately analysed, and three dimensional theories are required for acceptable accuracy. This is especially true in regions such as free edges and regions of stress concentration where the basic assumptions of CLPT are no longer valid. Transverse shear stresses are directly accounted for in the First Order Shear Deformation Plate Theory and higher order theories [START_REF] Kant | Estimation of transverse/interlaminar stresses in laminated composites -a selective review and survey of current developments[END_REF][START_REF] Reddy | An evaluation of equivalent-single-layer and layerwise theories of composite laminates[END_REF]. These Equivalent Single Layer (ESL) theories have limited accuracy in predicting interlaminar stresses because of their failure to account for both Zig-Zag effects (rapid change of slope across layer interfaces due to through the thickness discontinuity of mechanical properties) of the displacement fields in the thickness direction and interlaminar continuity of the transverse stress field [START_REF] Carrera | Theories and finite elements for multilayered anisotropic, composite plates and shells[END_REF]. Thus, for accurate evaluation of interlaminar stresses, three-dimensional elasticity is the tool of choice. This is reflected in the profusion of layerwise (LW) theories [START_REF] Reddy | An evaluation of equivalent-single-layer and layerwise theories of composite laminates[END_REF][START_REF] Carrera | Theories and finite elements for multilayered anisotropic, composite plates and shells[END_REF][START_REF] Carrera | Theories and finite elements for multilayered plates and shells: A unified compact formulation with numerical assessment and benchmarking[END_REF] and efficient solid shell elements [START_REF] Kant | Estimation of transverse/interlaminar stresses in laminated composites -a selective review and survey of current developments[END_REF][START_REF] Carrera | Theories and finite elements for multilayered anisotropic, composite plates and shells[END_REF].

Within the class of three-dimensional elements, an important issue is the accurate calculations of stresses. Stress recovery is highly dependent on the design of the particular element and on the formulation used in deriving the element. Mixed formulations provide more flexibility in the construction of efficient elements compared to classical formulations. Although primarily motivated by avoidance of the locking phenomena in finite element solutions, mixed formulations based on the Reissner or Hu-Washizu functional [START_REF] Washizu | Variational methods in elasticity and plasticity[END_REF][START_REF] Felippa | On the original publication of the general canonical functional of linear elasticity[END_REF], with displacements, stresses and/or strains as field variables, are so closely related to stress recovery that the formulations can also be viewed as stress recovery methods [START_REF] Mota | On mixed finite element formulations and stress recovery techniques[END_REF]. Mixed Layer-Wise theories are exemplified by the work of Carrera [START_REF] Carrera | Evaluation of layer-wise mixed theories for laminated plates analysis[END_REF], where a Reissner's mixed variational equation [START_REF] Reissner | On a certain mixed variational theorem and a proposed application[END_REF][START_REF] Reissner | On a mixed variational theorem and on a shear deformable plate theory[END_REF] is used to derive the governing equations in terms of displacements and transverse stress variables [START_REF] Carrera | Classical and advanced multilayered plate elements based upon pvd and rmvt. part 1: Derivation of finite element matrices[END_REF][START_REF] Carrera | Classical and advanced multilayered plate elements based upon pvd and rmvt. part 2: Numerical implementations[END_REF]. On the one hand assumed stress methods do not require any post-processing for the accurate recovery of stresses. On the other hand, inverse constitutive relations which define strains in terms of stress measures are employed greatly complicating the extension to the nonlinear case. Another shortcoming is the great added burden of solving simultaneously for stress and displacement degrees of freedom.

A particular form of the three field Hu-Washizu variational principle has been presented by Simo and Rifai [START_REF] Simo | A class of mixed assumed strain methods and the method of incompatible modes[END_REF], who developed a class of assumed mixed finite element methods to allow the systematic development of low order elements possessing good coarse mesh and distortion insensitive properties. The strain approximation is split into two parts: one is the usual displacement-gradient term, and the second is an enhanced strain part. Of particular note is the solid shell element developed by Vu-Quoc and Tan [START_REF] Vu-Quoc | Optimal solid shells for non-linear analyses of multilayer composites. i. statics[END_REF][START_REF] Vu-Quoc | Optimal solid shells for non-linear analyses of multilayer composites. ii. dynamics[END_REF]. For these three-dimensional enhanced strain elements, C 0 interpolated displacement fields are used. Physically, interlaminar stresses are continuous both across and along layer interfaces. Nonetheless, the continuity of interlaminar stresses is difficult to enforce in C 0 interpolated elements. Stress recovery procedures are necessary to find acceptably accurate stress fields.

Many researches have focussed on developing reliable stress recovery procedures. Early attempts to construct effective procedures included: the interpolation-extrapolation from super-convergent points [START_REF] Barlow | Optimal stress locations in finite element models[END_REF], the L 2 projection [START_REF] Oden | On the calculation of consistent stress distributions in finite element approximations[END_REF], the stress smoothing [START_REF] Hinton | Local and global smoothing of discontinuous finite element functions using a least squares method[END_REF][START_REF] Chen | Interfacial stress estimation using least-square extrapolation and local stress smoothing in laminated composites[END_REF], and the integral stress technique [START_REF] Argyris | Some considerations for the evaluation of finite element models[END_REF]. Zienkiewicz and Zhu [START_REF] Zienkiewicz | The superconvergent patch recovery and a posteriori error estimates. part 1: The recovery technique[END_REF][START_REF] Zienkiewicz | The superconvergent patch recovery and a posteriori error estimates. part 2: Error estimates and adaptivity[END_REF] made a very significant breakthrough towards an efficient postprocessing technique when they proposed the super-convergent patch recovery (SPR) procedure.

Many investigators have modified this procedure to include satisfaction of boundary conditions [START_REF] Blacker | Superconvergent patch recovery with equilibrium and conjoint interpolant enhancements[END_REF][START_REF] Wiberg | Enhanced superconvergent patch recovery incorporating equilibrium and boundary conditions[END_REF][START_REF] Lee | A super-convergent stress recovery technique with equilibrium constraint[END_REF]. For a more accurate prediction of the transverse stresses in laminated composites and shells, a Modified Super-convergent Patch Recovery (MSPR) technique has been derived to obtain accurate nodal in-plane stresses which, subsequently, are used in the integration along the thickness of the equilibrium equations for evaluating the transverse shear and normal stresses [START_REF] Yuan | Stress projection, layerwise-equivalent formulation for accurate predictions of transverse stresses in laminated plates and shells[END_REF].

Patch recovery procedures depend crucially on the accurate evaluation of stresses within elements.

Stress fields within an element can be deduced using either constitutive relations or variationally consistent procedures. Stress fields deduced using constitutive relations can show spurious oscillations due to the retention of higher order (inconsistent) terms which do not contribute to the determination of the displacements and strains, which then get reflected as extraneous stress oscillations [START_REF] Prathap | Consistent force resultant distributions in displacement elements with varying sectional properties[END_REF][START_REF] Prathap | Consistent thermal stress evaluation in finite elements[END_REF]. As far as variationally consistent procedures are concerned, it is worth mentioning the weighted residual approach developed by De Miranda and Ubertini [START_REF] De Miranda | Recovery of consistent stresses for compatible finite elements[END_REF]. Nevertheless, even variationally consistent procedure suffer from spurious oscillations in stress fields as shown by Dakshina-Moorthy and Reddy [START_REF] Dakshina Moorthy | Recovery of interlaminar stresses and strain energy release rates in composite laminates[END_REF].

The purpose of this paper is to present an efficient interlaminar stress recovery procedure within three dimensional FE formulations that is able to overcome the problems usually encountered using customary procedures, i.e. severe oscillations of stress distributions, without relying on extrapolation techniques from super-convergent points. Thus, the accuracy of the recovered interlaminar stresses is not dependent on the knowledge of superconvergent point nor is it sensitive to the stress recover method employed to obtain element stress distributions. Moreover, the procedure attains as much accuracy as assumed stress methods using coarser meshes and without the need to include stress degrees of freedoms in the solution process. The proposed procedure was developed based on equilibrium considerations, and interlaminar stresses are recovered directly at nodes. The numerical formulation is presented and explained in detail in the next section. The reliability of the approach is tested analysing several benchmark problems, and the results are reported in section 3. Conclusions and future developments are addressed in section 4.

Interlaminar stresses recovery procedure

The present procedure is an extension of the interlaminar stress recovery procedure developed by Dakshina Moorthy and Reddy [START_REF] Dakshina Moorthy | Recovery of interlaminar stresses and strain energy release rates in composite laminates[END_REF] in the context of studying delamination in multilayered composites. In their approach, each ply is modelled as a separate body and the interlaminar boundary is treated as a contact surface. The interlaminar forces are obtained using an interface model based on the penalty method. Using these contact loads the interlaminar stresses are recovered.

The recovery procedure partitions the contact surface into a set of non-overlapping patches corresponding to groups of elements. The traction distribution is interpolated over each patch in terms of nodal values. Static equivalence between the tractions and the contact forces is used to calcu-late the nodal value of the tractions, hence the interlaminar stresses. The procedure was found in practice to lead to oscillatory interlaminar stresses. Thus, a final step is the application of a smoothing technique to obtain more physically meaningful interlaminar stresses.

The present procedure extends the procedure of Dakshina Moorthy and Reddy in two ways. First, each ply is still modelled separately but the compatibility between the plies is enforced using Lagrange multipliers following the procedure commonly employed in domain decomposition methods.

Thus, connecting forces are obtained without resorting to a penalty formulation. Second, the traction distribution is interpolated over the complete interlaminar surface. Thus, the full distribution of interlaminar stresses is obtained simultaneously. In the following, we give more details about the proposed formulation.

Consider a typical interface between two plies as shown in figure 1. The system of linear equilibrium equations and compatibility conditions can be written as,

      K (1) 0 B (1) T 0 K (2) B (2) T B (1) B (2) 0             u (1) u (2) λ       =       f (1) f (2) 0       (1)
where K (1) and K (2) are the stiffness matrices of the two plies, u (1) and u (2) the displacement vectors of the two plies, f (1) and f (2) the load vectors of the two plies, B (1) and B (2) are the Boolean matrices providing the equality conditions between displacements of the two plies, and λ are the Lagrange multipliers (interlaminar forces) introduced to enforce the compatibility constraints between the plies. Solving the system of equations ( 1), the interlaminar forces between the plies are obtained as well as displacements.

The traction distribution t over a surface element is interpolated using C 0 iso-parametric shape functions,

t = N t f (2)
where t f is the vector of the nodal traction values and N is the matrix of C 0 shape functions. The surface element connectivity is directly inherited from the solid element faces.

Displacements over a surface element are interpolating using the same C 0 iso-parametric shape functions. Thus, the internal work W in done by the tractions over the element surface can be written as,

W in = Ω f ( t f N ) T traction N ∂q f disp. ∂Ω f = t T f M f δq f . (3)
where q f is the vector of nodal displacements, and M f is the matrix of areas,

M f = Ω f N T N ∂Ω f . (4)
The surface element contributions are assembled over the complete interlaminar surface in a matrix M , and considering the static equivalence between the tractions and the connecting forces, the following relation is obtained,

M t = λ ( 5
)
where t is the global vector of nodal traction. The process is repeated for every interlaminar surface of interest. The accuracy of the present approach is demonstrated in the next section.

Numerical Results

Simply supported plate subjected to bisinusoidal pressure load

The present procedure was validated considering a benchmark problem analysed by Pagano [START_REF] Pagano | Exact solutions for rectangular bidirectional composites and sandwich plates[END_REF], where three dimensional exact elasticity solutions of idealised simply supported cross-ply

[0 • /90 • /0 • ],
square plates under bisinusoidally distributed pressure load of intensity p z , are provided. The length and thickness of the plate is denoted by "a"and "H"respectively. The laminate is made of material plies that are idealized to be homogeneous, elastic and orthotropic. The following material properties are used:

E 11 = 25 GP a, E 22 = E 33 = 1 GP a, G 12 = G 13 = 0.5 GP a, G 23 = 0.
2 GP a, and ν 12 = ν 13 = ν 23 = 0.25. Subscripts 1,2 and 3 denote the fibre, transverse and thickness directions, respectively. The rectangular Cartesian coordinate system used is such that the origin is located at the center of the middle layer of the plate. Stresses are normalized according to the following formulae,

(τ ′ xz , τ ′ yz ) = 1 p z S (τ xz , τ yz ) , σ ′ zz = 1 p z σ zz ,
where the S is the laminate length to thickness ratio S = a/H.

The discrete model is generated using a simple linear solid-shell element developed by Vu-Quoc and Tan [START_REF] Vu-Quoc | Optimal solid shells for non-linear analyses of multilayer composites. i. statics[END_REF][START_REF] Vu-Quoc | Optimal solid shells for non-linear analyses of multilayer composites. ii. dynamics[END_REF]. This element possesses good coarse mesh and distortion insensitivity properties for a large range of aspect ratios. Unless differently stated, the in-plane mesh consists of 14 × 14 elements. The minimum required number of elements in the thickness direction is three. While this is sufficient for accurate interlaminar stresses, it doesn't allow the detailed evaluation of their through the thickness variation. In the case under consideration the values of transverse stresses in the middle of every layer are also reported and, as a consequence, three more mathematical interfaces are required in these locations to give a total number of six elements through the thickness.

Laminates of length to thickness ratios S = 20, 50, 100 are analysed and the results calculated at several points of interest are reported in table 1. Compared to Pagano's exact solutions, excellent agreement is obtained for all the considered cases. Through-thickness distributions obtained using ABAQUS 6.8 TM are also reported for comparison. The adopted elements are respectively, the 8-node linear brick incompatible modes C3D8I, the 20-node quadratic brick C3D20R, and the linear 8-node C3D8R. Depending on the element formulation, stresses are recovered at integration points using constructive models derived from either variational principles or other energy laws [START_REF]Abaqus analysis user's manual-abaqus version 6[END_REF]. Stresses are extrapolated at nodes successively using the shape functions. Both the linear C3D8R and C3D8I elements do not provide single valued results at the interfaces. A through-the-thickness refinement of the mesh is required.

The quadratic C3D20R element does not produce an accurate estimation of the transverse shear stress τ ′ yz . The percentage error at the point of maximum value, i.e. z = 0, compared to the exact solution is 9.46%.

The aforementioned results are obtained adopting one additional mathematical interface in the middle of every ply. However, these mathematical interfaces are not required to have appropriate interlaminar stress distributions using the present procedure.

The averaged applied stress is calculated as:

σ 0 = 1 aH a 0 H 0 σ xx (a, y, z) dy dz ( 7
)
where H is the complete laminate thickness.

Owing to the symmetry of loading and lay-up, only half of the total thickness of the laminate is modelled. The in-plane mesh is shown in figure 8, and different number of elements along the thickness direction are considered. The designation [ner/rsrnet] indicates: "ner"is the number of elements extending in a radial direction from the hole edge to the end of the circular region, r = 80 mm, "rsr"is the radial spacing ratio of the elements in the circular region, and "net"is the number of elements through the half thickness. The number of elements around a quarter of the hole edge is 18.

Interlaminar stress distributions along the radial r coordinate at θ = 90 • , evaluated in a cylindrical coordinate system r, θ, z with the origin in the center of the plate, are shown in figure 9. These distributions are plotted starting from the hole edge up to one laminate thickness r = H + D/2, at the 45 • /-45 • interfacial surface, and compared with the stress fields provided by Iarve and Pagano (not reported) in the singular neighborhood of the ply interface and the hole edge. The present approach shows excellent agreement with Iarve and Pagano's distributions for a distance from the hole edge greater than one ply thickness. Convergence is reached using coarse meshes (mesh labelled as 24/30 -2 solves a system of 69120 linear equations), however, a mesh-dependent influence of the stress singularities is encountered within one ply-thickness from the hole edge.

This influence is present in terms of oscillations that seem to be reduced introducing mathematical interface in the middle of every ply as far as the interlaminar shear stress τ θz is concerned, and using a refined mesh close to the hole edge as far as the interlaminar normal stress σ zz is concerned.

These oscillations may be related to the high value of aspect ratio of the solid-shell element used for meshing this region, and probably generating some inaccuracy in presence of high gradients.

To understand how this mesh-dependent behaviour of the stress singularities effects the reliability of the present procedure close to the hole edge, a symmetric cross ply [90 • /0 • ] s laminate similar to the previous one and subjected to an uniaxial tensile load σ 0 was also considered. This laminate has been analysed by Hu et al. [START_REF] Hu | Interlaminar stresses in composite laminates with a circular hole[END_REF] using a three-dimensional finite-element (FE) analysis based on displacement formulation employing a curved isoparametric 20-node element. The total length of the panel is 60 mm, the total width 30 mm, the hole radius R is 2.5 mm, and the ply thickness h is equal to 0.125 mm. Material properties are the same used for the [45 • / -45 • ] s previously analysed.

Owing to the symmetry of loading, hole location and lay-up, only one-eight of the laminate is modelled. The in-plane mesh structure is the same used previously, with the addition of a coarse parte extended in the x direction since the laminate is not anymore square. Henceforward, if not differently stated, the mesh adopted solves a system of 43200 linear equations. In their work, Hu et al. adopt a FE model consisting of 4000 elements and solving a system of 56000 linear equations.

They provide stress distributions along radial lines away from the hole and around the hole at the 90 • /0 • interface, obtained by averaging the 90 • and 0 • ply values at the interface, a customary finite element practice. These distributions are obtained near but not at the hole edge because of the mathematical interlaminar stress singularity.

As far as the present procedure is concerned, the radial σ zz distributions away from the hole at three angular positions θ = 0 • , 45 • , 90 • are presented in figure 10(a). The σ zz distributions have step gradients near the hole edge and approach to zero within a ply thickness from the hole. Due to the more refined mesh adopted compared to the one used in the [45 • / -45 • ] s laminate, oscillations close to the hole edge are not anymore present, but smooth behaviours are generated.

These distributions are in excellent agreement with the results (not reported) of Hu et al. [START_REF] Hu | Interlaminar stresses in composite laminates with a circular hole[END_REF]. The maximum difference between the two distributions is within 4% at θ = 75 • . The shear stress component τ rz is very small compared to τ θz and can be neglected.

The minimum distance r/R = 0.000082 from the hole to report the circumferential interlaminar stress distributions was selected based on two considerations. Firstly, there is a limit in the refinement of the mesh close to the hole edge since the aspect ratio of the solid element used can be increased until a certain threshold. Secondly, stress oscillations begin after the considered point and a certain error in the distributions can be introduced. To verify if these oscillations effect the prediction of delamination initiation, a stress failure criteria needs to be selected.

From the literature [START_REF] Brewer | Quadratic stress criterion for initiation of delamination[END_REF] it can be seen that the approach for predicting failure in such laminates has been that of averaging the interlaminar stresses over a distance from the hole edge, suggesting that the exact values of the stresses at the free edge are not too important. As failure stress criterion was assumed the Tsai-Wu criterion [START_REF] Tsai | Introduction to composite materials[END_REF] which takes into account the interaction of all six stress components in a quadratic equation. Since it is assumed that the delamination initiation is mainly attributed to interlaminar stress effects, only the interlaminar stresses are considered. Moreover, since the interlaminar shear stress τ rz can be neglected compared to τ θz , the Tsai-Wu criterion in a cylindrical coordinate system can be simplified as,

σ zz Z 2 + τ θz S 2 = e 2      e < 1 no failure e > 1 failure (8)
where Z is the interlaminar normal strength and S is the interlaminar shear strength. For positive interlaminar normal stress σ z the uniaxial tensile strength Z t should be used while for negative σ z the compressive strength Z c should be employed.

The average stress failure criteria assumes that delamination initiates when the stresses at a char-acteristic distance a 0 from the discontinuity meet the failure criteria [START_REF] Mota | On mixed finite element formulations and stress recovery techniques[END_REF]. The average of a stress component is defined as,

σ ij = 1 a 0 R+a 0 R σ ij dr (9)
The following strength properties, Z t = 50.6M P a, Z c = 200M P a and S = 103M P a, were used for the carbon fiber-epoxy laminates. Figure 12 It is worth mentioning that the characteristic distance a 0 is experimentally determined and can vary with lay-up and hole size. Average stress failure criteria are appropriate for predicting onset of delamination but they are not suitable for failure analysis of composites, properly analysed using fracture mechanics based methods instead [START_REF] Kashtalyan | Strain energy release rate for off-axis ply cracking in laminated composites[END_REF][START_REF] Kashtalyan | Analysis of local delaminations induced by angle ply matrix cracks[END_REF].

Conclusion

An efficient interlaminar stress recovery procedure for three-dimensional finite elements formulations is presented. This procedure, based on equilibrium considerations, retrieves interlaminar stress values directly at the interfacial nodes between the plies and stress continuity at the interelement boundary is automatically satisfied. As a consequence, smooth distributions are easily generated. Plates of various geometries were considered and the results were compared with available exact and Finite Element solutions. These comparisons indicated that excellent agreement was obtained with exact solutions and convergence was reached using considerably less degrees of freedom compared to other finite elements procedures. Consequently, the procedure is more suitable for design purposes. Special emphasis is placed on the problem of loaded plate with an open circular hole. The present procedure is able to produce converged mesh-independent average interlaminar stresses even though some oscillations are encountered close to the hole edge. As such it can be effectively combined with average failure stress criteria for the prediction of delamination initiation in presence of curved free edges and stress concentrations. The present procedure is currently being implemented in the commercial Finite Element software ABAQUS 6.8 TM .

A

  laminate of length to thickness ratio S = 100 is considered in detail to demonstrate the efficacy of the method for thin laminates. The convergence of the interlaminar stresses τ ′ xz and σ ′ zz at the [0 • /90 • ] interface is plotted in figure 2. Excellent agreement with Pagano's exact solution and fast convergence are obtained. The percentage error at points of maximum values compared to the exact solution is within 0.2%. Smooth distributions of interlaminar stress are observed as shown in figure 3. The same laminate has been analysed in the work of Dakshina-Moorthy and Reddy [31]. Regarding the structural symmetry, only a quarter of the plate is taken into account. An in-plane mesh of [8 × 8] nine-node quadratic elements and a mesh of linear Lagrange elements through the thickness of the laminate is considered. The interlaminar stresses are calculated and compared using three different approaches: (1) their equilibrium based stress recovery procedure, at nodes; (2) from constitutive relations, at integration points; and (3) variationally optimal stress recovery procedure, at integration points. As shown in figure 4, all the aforementioned procedures generate severe oscillations. Smoothing techniques are required. Comparing these distributions with the ones obtained using the present procedure it can be stated that the present procedure is able to circumvent the shortcomings of the recovery procedures developed by Dakshina-Moorthy and Reddy. A laminate of length to thickness ratio S = 50 is also considered to demonstrate the accuracy of the proposed procedure in analyzing moderately thick plates. Through-thickness distributions of interlaminar stresses produced using the proposed method, labelled as Equilibrium, are compared with Pagano's exact solution in figure 5 at points of major interest. Excellent agreement and fast convergence are again achieved.

Figure 6

 6 shows surface plots of the interlaminar shear stress τ ′ xz at the [0 • /90 • ] interface obtained respectively without and with additional mathematical interfaces. An almost identical response is generated, accuracy and smoothness are satisfied regardless of the introduction of additional mathematical interfaces, implying that convergence is only related to the in-plane mesh. The present procedure uses a model with 5400 degrees of freedom and attains converged interlaminar stress results, whereas the ABAQUS's model using quadratic elements has 17595 degrees of freedom and still fails to converge. Stress concentration phenomena play an important role in the design of layered structures. A well known stress concentration problem is the problem of a plate loaded in-plane and containing a circular open hole. This problem combines strong in-plane stress gradients with free edge effects and is characterized by the occurrence of strongly three-dimensional singular stress fields at the free edges in the interface between two layers of composite laminates. As a consequence, stress calculation at the interlaminar surfaces in the vicinity of the hole edge is a difficult problem. The accuracy of the present recovery procedure is tested examining a square [45 • / -45 • ] s plate with a circular hole as shown in figure 7. This problem has been analysed by Iarve and Pagano [34] using a method developed for the superposition of a hybrid and displacement approximation. They use asymptotic analyses to derive the hybrid stress functions and the displacement approximation is based on the polynomial B-spline functions. The laminate consist of the planar dimensions a = 508 mm, diameter of the hole D = 50.8 mm, coordinate of the hole center x c = y c = a/2, ply thickness h = 2.54 mm. Each ply is treated as an homogeneous, elastic and orthotropic material with the following properties: E 11 = 138 GP a, E 22 = E 33 = 14.5 GP a, G 12 = G 13 = G 23 = 5.86 GP a, and ν 12 = ν 13 = ν 23 = 0.21. The uniaxial loading u 0 /a = 0.001 is applied via displacement boundary conditions at the lateral sides (x = 0, a),-u x (0, y, z) = u x (a, y, z) = u 0 , u y (0, y, z) = u y (a, y, z) = 0.

Figure 10 (

 10 Figure 10(b) shows the interlaminar normal stress distributions at the 90 • /0 • ply interface around the hole. As a mathematical interlaminar stress singularity exists at the free edge between the 90 • and 0 • plies, the computed stresses are presented near but not at the hole boundary. As the distance from the edge (r -R) increases, the interlaminar stress σ zz is rapidly decreased. When (r -R) = 0.1R, i.e. two-ply thickness away from the hole boundary, σ zz becomes almost zero. σ zz is compressive for most of the region around the hole with a small tensile region near θ = 90 • . The largest compressive σ zz occurs at about 60 • from the loading axis. The distribution obtained by Hu et al. at (r -R)/R = 0.0001 is also plotted for comparison. The maximum difference between the two distributions is within 3% at θ = 60 • . The radial variation of τ θz at θ = 10 • , 45 • , 75 • is shown in figure 11(a). The decay ratio to zero varies with θ; the maximum decay ratio is at θ = 75 • and, as pointed out in the plot, oscillations are still present really close to the hole edge for all the distributions. Nevertheless, these distributions are in very good agreement with the ones obtained by Hu et al. The circumferential interlaminar shear stress distributions τ θz at various distances, (r -R)/R, from the hole boundary are shown in figure 11(b). Similar to the normal stress, σ zz , the interlaminar shear stress τ θz decreases as the distance (r-R) from the hole boundary increases and becomes small within two-ply thickness

  shows the distributions of the e-index around the hole obtained assuming different values of a 0 . The maximum value of e varies with a 0 and occurs in the region θ = 67.5 • -75 • to the loading direction, indicating that these are the critical locations for delamination growth. e-indexes distributions obtained by Hu et al. are also reported for comparison. Slightly different values of the failure index e are obtained assuming a 0 = 0.0005 between 70 • and 80 • , but the percentage difference is always within 3%. Consequently, the oscillations of the interlaminar shear stress τ θz close to the hole edge do not seems to influence the failure index distribution. The introduction of mathematical interfaces in the middle of every layer produces negligible variations on the results, indicating that convergence is already reached.

Figure 1 :

 1 Figure 1: Connecting forces λ f at the interfacial nodes of each element obtained by the FETI method

  Convergence of the interlaminar normal stress σ ′ zz with different in-plane meshes.

Figure 2 :

 2 Figure 2: Convergence of the recovered interlaminar transverse stresses at the [0 • /90 • ] interface and y = 0 of a [0 • /90 • /0 • ] laminate, S = 100.

  Surface plot of the interlaminar shear stress τ ′ xz . Surface plot of the interlaminar normal stress σ ′ zz .

Figure 3 :

 3 Figure 3: Surface plot of the recovered interlaminar stresses at the [0 • /90 • ] interface of a [0 • /90 • /0 • ] laminate, S = 100.

Figure 4 :

 4 Figure 4: Shear stress distribution τ ′ xz at the interlaminar surface [0 • /90 • ] and y = 0 of a [0 • /90 • /0 • ] laminate obtained by Reddy and Dakshina-Moorthy [31] considering different approaches, S = 100.

Figure 5 : 5 ( 5 (

 555 Figure 5: Comparisons between through-thickness distributions of the transverse stresses obtained using the present equilibrium based procedure, using different ABAQUS's elements, and Pagano's exact solution, S = 50.

Figure 6 :Figure 7 :YFigure 8 : 10 Figure 9 :

 678109 Figure 6: Surface plot of the recovered interlaminar shear stress τ ′ xz at the [0 • /90 • ] interface of a [0 • /90 • /0 • ] laminate, S = 50.

  R)/R = 0.01 (r -R)/R = 0.001 (r -R)/R = 0.0001 (r -R)/R = 0.00008 Hu et al. (r -R)/R = 0.0001 (b) Normalized interlaminar normal-stress distributions around the hole.

Figure 10 :

 10 Figure 10: Normalized interlaminar normal stress distributions σzz/σ0 at the 90 • /0 • interface of a [90 • /0 • ]s symmetric laminate.

  interlaminar shear stress distributions at different angular positions.

  R)/R = 0.1 (r -R)/R = 0.01 (r -R)/R = 0.001 (r -R)/R = 0.0001 (r -R)/R = 0.00008 Hu et al.(r -R)/R = 0.0001 (b) Normalized interlaminar shear stress distributions around the hole.

Figure 11 :

 11 Figure 11: Normalized interlaminar shear stress distributions τ θz /σ0 at the 90 • /0 • interface of a [90 • /0 • ]s symmetric laminate.

Figure 12 :

 12 Figure 12: e-index distributions at the 90 • /0 • interface of a [90 • /0 • ]s symmetric laminate near the hole edge determined by the average stress criterion.

Table 1 :

 1 Recovered stress fields considering different values of length to thickness ratio S.

	S	τ ′ xz (-a 2 , 0, 0) τ ′ yz (0, -a 2 , 0) σ ′ zz (0, 0, H/6)
	20		0.3846	0.0938	0.7398
		EQUILIBRIUM	0.3845	0.0938	0.7429
	50	EXACT	0.3934	0.0842	0.7406
		EQUILIBRIUM	0.3928	0.0841	0.7436
	100	EXACT	0.3946	0.0828	0.7407
		EQUILIBRIUM	0.3926	0.0823	0.7434
		Connecting	W f	W
		forces l f

Acknowledgments

The useful collaboration of Prof. D.J.Rixen (Delft University of Technology, Department of Mechanical,Maritime and Materials Engineering) in the implementation of the Finite Element Tearing and Interconnecting method, is gratefully acknowledged.