Natural disasters, such as volcanoes, floods, landslides, hurricanes, earthquakes, and climatic changes, represent a real danger for the conservation of the existing cultural heritage. Nevertheless, insufficient programs have been applied, aimed at reducing the vulnerability, and the related risk of cultural heritages. An idea of the problem dimensions is given by analyzing the costs required in the last natural disasters in order to manage the emergencies, repairs, and reconstruction. According to an estimation provided by World Bank (World Bank Indipendent Evaluation Group, 2006), natural disaster damage costs are increasing, and they have achieved about 652 billion US dollars in the 1990s. These costs result 15 times higher than the ones registered in the 1950s, even for natural disasters. Another piece of important information is obtained: if one considers that one-third of the cost to the global economy (about 50 billion US dollars) is spent on predicting, preventing, and mitigating disasters and the other two thirds represent the direct costs of the damage (Alexander, 2017).

The ability to react to these disasters is directly related to the public awareness about the future destructive events that may occur. In this context, it is necessary to design and to apply prevention measures addressed for protecting the cultural heritage, where experts of different fields (engineering, statistics, chemistry, seismology, etc.) are involved. Thanks to this cooperation, different documents have been proposed such as, among the others, the Mexican risk identification atlas (CENAPRED, 2014), and the territorial information systems proposed in ISCR (2017).

In line with these premises, this paper presents a comparative study of the seismic performance of ancient masonry churches. To this scope, the simplified methods reported within the Italian Directive (GU. n. 47, 09/02/2011) are considered, with particular attention to Level of Valuation 1 (LV1, for an evaluation at a territorial dimension) and Level of Valuation 2 (LV2, considering the macro-elements approach). Furthermore, a new simplified method proposed in (CENAPRED, 2006) and validated in D'Amato et al. (2018), useful at a territorial scale for preliminary scoring and classifying the seismic behavior of churches, is applied. This methodology has also been recently validated for ancient adobe Chilean Churches (Fuentes et al., 2019).

The seismic performance evaluation in this study is conducted by considering in detail the case study of the Cathedral of Matera (Italy), an ancient masonry church dedicated to SS Maria Della Bruna. The obtained results are then compared with other churches, similar to the Cathedral with respect to materials, geometrical characteristics, and constructions details. All the churches considered fall within the historic centre of Matera, named “Sassi of Matera,” recognized since 1993 as a World Heritage Site by UNESCO.

The main feature of the considered methods is that they can be used within a multilevel seismic vulnerability approach. In fact, they are based only on simple surveys characterized by visual and geometric detections, implying limited costs. Thus, these methods can also be applied as decision making tools in order to rank priorities and to proceed to further material and structural investigations for performing advanced structural analyses. Within this multi-level approach, the new simplified method recently proposed can be intended as a “Level of Valuation 0” (“LV0 method”), by allowing a very rapid seismic assessment, even at a larger territorial scale. Moreover, the comparisons among the obtained results show that the simplified methods may also overestimate the seismic actual response. Therefore, they remain useful for comparing and ranking the case studies, but they cannot substitute more refined methods for realistically simulating the seismic behavior of churches.

The Italian directive, specific for Cultural Heritage (G. U. n. 47, 26/02/2011), defines three levels of valuation, namely the LV1, LV2, and LV3 method for assessing the seismic performance of different historical constructions. These methods have an increasing complexity, requiring in parallel an increasing amount of information regarding structural details and materials properties.

As regards to the LV1 method, it is suitable when a comparative seismic performance at a territorial level is conducted. Basically, the idea of this method is to relate the seismic performance to a global vulnerability index of a structure. No intervention may be designed with this method. Recent applications and similar approaches to this method may be found, for instance, in Lourenço and Roque (2006); Lourenço et al. (2013); Caprili et al. (2017); Formisano et al. (2017, 2018); Marotta et al. (2018); Sarhosis et al. (2018);Valente and Milani (2018a).

A more refined seismic performance evaluation may be obtained with the LV2 method. It consists of evaluating the local mechanisms of architectural parts of a structure, named “macro-elements,” and the global vulnerability corresponds to the activation of the most vulnerable mechanism considered. The macro-elements approach has been proposed after the damage surveys suffered by ancient churches during some recent Italian earthquakes. The conducted surveys have highlighted a systematic repetition of the detected damages, suggesting that the church structural system should be considered as a group of architectural parts, showing independent response mechanisms under seismic lateral loading. The macro-elements are identified as architectural parts (such as the main façade, lateral walls, longitudinal and transverse colonnade of main nave, triumphal arch, bell tower, dome, etc.), evaluating possible interactions with adjacent elements. Applications of this approach may be found, among the others, in Betti et al. (2018), Brandonisio et al. (2013), D'Ayala and Paganoni (2011), Doglioni et al. (1994), Formisano and Marzo (2017), Lagomarsino and Podestà (2004a,b), while a recent investigation on the estimation of main frequency of ancient masonry churches may be found in Lopez et al. (2019).

Non-linear finite element analyses are considered in the LV3 method (G. U. n. 47, 26/02/2011), and they may regard limited parts or a whole structure. The adopted models must reproduce the real distribution of stiffness and mass, as well as the material non-linear behavior. Of course, the results obtained may realistically simulate the structural behavior only if the numerical assumptions made are true. This is strictly dependent on the knowledge of structural details (such as, for instance, wall connections, connections with the roof parts, and clear knowledge of the structure evolution during the past). The amount of these available data, and consequently, the reliability of the obtained results may increase if an in-depth and critical investigation campaign, through in situ tests and inspections, is performed. However, it should be carefully considered that in-situ tests should be focused at first on structural details more than on material mechanical properties, since the former dominate the response influencing the elements interaction more than the latter. Discussion about this issue may be found in Clementi et al. (2016), Castori et al. (2017), and Valente et al. (2017). In addition, the use of non-linear methods, such as the conventional pushover approach, implies the decisive choice of the control point on which the non-linear global response and the evaluation of the seismic performance depend. A discussion of this aspect may be found, among others, in Betti et al. (2018). While a general discussion on the application of non-linear analyses and the comparisons with the macro-element approach may be found in Castellazzi et al. (2013) and Valente and Milani (2018b).

Although it is the most simplified one, the LV1 method is useful for evaluations at a territorial scale capable of providing, through a vulnerability index, an estimation of the ground acceleration related to the collapse. It requires only visual inspections and a qualitative judgment of some structural details. On the other hand, a local analysis may be performed with the LV2 method (macro-element approach) also designing local interventions, where the potential interaction among the structural part may also be taken into account. The most refined approach is represented, without any doubt, by the LV3 method, where non-linear FE models are required. However, as discussed previously, this approach required a high amount of data and it may be extremely high time consuming from a computational point of view.

In the following, only the LV1 and LV2 methods are described in detail and, later on, are applied to some case studies. In addition, in this study a new simplified method for seismic assessment of the considered cases is also considered, proposed as LV0 in continuity with the Italian Directive approach. The obtained results will be compared, in order to establish if a general correspondence in terms of the predicted seismic performance among the examined methods exists.

In Dìaz (2016), a new simplified methodology capable of assessing the seismic risk of ancient masonry churches has been proposed. Subsequently, this methodology has been preliminarily validated in D'Amato et al. (2018), and recently also extended to the case of Chilean Churches in Fuentes et al. (2019).

The simplified methodology provides a seismic risk score R, applying the known symbol equation (UNDRO, 1979; FEMA, 2004):

where E, H, and V represent the exposure of elements or assets at risk, the seismic hazard, and the vulnerability, respectively.

The seismic risk score R is obtained by multiplying the scores E, H, and V, separately evaluated by applying three different tools, that are:

As proposed in Dìaz (2016), the simplified method may also be applied as the scoring method by neglecting the exposition value, that is the E score within the Equation (17). Hence, the seismic score R is given by:

obtained by considering only the scores separately obtained with the Tool 2 and the Tool 3. In the previous equation, H is increased as a unity for obtaining a seismic score R more major than 1. One may note that as R increases, the seismic risk also increases. The proposed method is particularly suitable for conducting, at a territorial scale, a seismic comparative analysis, with the simplicity of requiring few information to be applied. Instructions on how to calculate the vulnerability (V) and seismic hazard (H) scores may be found in Dìaz (2016) and in D'Amato et al. (2018).

The city of Matera is located at the south of Italy, in the region of Basilicata. Specifically, the SS Maria della Bruna church is located in the city center, called “Sassi of Matera,” consisting of two areas protected, since 1993, by UNESCO and reported within the World Heritage List. As shown in Figure 1, the church is characterized by a main body composed of three naves, having a latin cross configuration in plan. As in the Cristian tradition, the church main façade is oriented toward the west while the altar is oriented toward the east. The Cathedral construction started between 1226 and 1231 (Morelli, 1970). An inscription on the bell tower door reports that in 1270 the construction finished. Between the 15 and 16th centuries some chapels were annexed to the main body on the north side. As for the bell tower, two distinct parts may be clearly identified, since the upper zone was built later than the rest of the tower but not after 1709, given that this part appears in a fresco. The choir was originally completed in 1729 and, later, completely reconstructed in 1738. Additionally, it is supposed that due to a partial collapse, the dome was reconstructed. Regarding the bell tower, it is possible that the upper section was built later than the rest of the tower but not after 1709, given that this part appears in a fresco after this year (Ragone et al., 2017). The façade of the Cathedral presents several ornamental elements with religious meanings. Above the main entrance, a statue of SS Maria della Bruna, the Saint Patron of the City, is situated. At the top of the main façade, a finely decorated oculus is situated. A plan of the church is reported in Figure 1.

The results obtained from the considered case study are also compared in this work with four similar ancient masonry churches, all located within the “Sassi of Matera” area. The four additional churches, examined in a previous work (D'Amato et al., 2018) are: San Pietro Caveoso, San Rocco, San Francesco d'Assisi, and San Giovanni Battista. The considered churches are illustrated in Table 1, where for each of them, the main facade, an interior view and the floor plan is reported. For each considered church, original drawings reporting the geometrical dimensions were retrieved. In addition, specific detailed surveys were conducted in order to control any possible existing cracking pattern, material deterioration or structural instability. No in-situ tests were carried out since the calculation methods applied in this study do not depend on the mechanical properties of the involved materials. In the following, a summary of the description of the considered churches is reported.

With the aim of requalifying the highest number of manufacts of the Italian cultural heritage, the Italian directive (G.U. n. 47, 26/02/2011) allows for improving the seismic performance without imposing a complete retrofitting. This consequently allows us to design small and cheap interventions, instead of heavy, widespread and invasive ones. In this way each manufact can be considered as seismically protected for a lower action level (i.e., for a lower Return Period, T_R). In other words, it is permitted to reduce the observation time in which the seismic action is evaluated. This interval, defined as reduced nominal life V_{N, Red}, has to be assumed equal to or higher than 20 years (G.U. n. 47, 26/02/2011).

In line with this criterion, the seismic analyses presented herein are performed considering a reduced nominal life V_{N, Red} assumed corresponding to 20 years, and a standard nominal life V_N = 50 years, as requested for newly designed buildings. The coefficient of use adopted in both the cases is C_U = 1.5, and, therefore, the observation time is adopted to define the seismic action results V_R = 30 years for V_{N, Red} = 20 years and V_R = 75 years for V_N = 50 years.

The city of Matera falls in a moderate seismicity zone, and many moderate seismic events hit it in the past years. Specifically, in the Italian historic catalog (CPTI, 2015), some events with macro-seismic intensity equal to 7 and many events with intensity equal to 6 are found. The seismic action is defined in accordance with the Italian seismic code (Ministerial Decree, 14/01/2008) by considering the Damage Limit State (DLS) and the Life Safety Limit state (LSLS). By assuming the observation periods discussed above, we obtain a return period T_R equal to 30 and 75 years for the DLS, and T_R equal to 285 and 712 years for the LSLS.

In Figure 2 the considered elastic and design spectra for the horizontal component of seismic action are shown. It should be noted that, according to what observed in Paulay and Priestley (1992) for masonry structures, the elastic spectra is calculated by considering a damping ratio equal to 10%. Furthermore, a behavior factor q = 2 is considered for the design spectra, according to the Italian code suggestions.

In order to apply the LV1 method, Table 2 numerically reports, for each macro-element considered, the score assigned to the actual vulnerability (v_ki) and seismic-resistant device (v_Kp). With these scores, it has been possible to calculate for each church the vulnerability index i_v (Equation 1), and then the seismic capacity measured through a_LSLS and a_DLS (Equations 2, 3). The obtained results are summarized in Table 3, where they are reported as the acceleration factor f_a (Equation 7) and the seismic safety index I_S (Equation 4). These parameters are calculated by considering two observation periods, that are V_R = 30 years (obtaining I_{s, 30} and f_{a, 30}) and V_R = 75 years (obtaining I_{S, 75} and f_{a, 75}). The seismic verifications are conducted by referring to both the LSLS (having a P_VR = 10 %) and to the DLS (having instead a P_VR = 63 %). Moreover, in comparing the results, the following additional assumptions have been made: FC = 1.35 (factor of confidence), S = 1 (soil stratigraphic factor).

Firstly, the results are presented in terms of return periods (T_LSLS and T_DLS, for LSLS and DLS, respectively) obtained on the predicted accelerations (a_LSLS/FC and a_DLS/FC) of each case study. It must be noted that the expected return period T_DLS results are always, except for San Pietro Caveoso church, lower than 30 years, corresponding to the first point available of the site seismic hazard (Figure 2B). Therefore, in these cases, by following the Italian Design Code, a_{g, DLS30} has been considered equal to the one provided for a return period of 30 years (in this case it results a_{g, DLS30} = 0.038g). For sake of completeness, Figure 3 illustrates the numerical results summarized in Table 3 for the two considered limit states.

Figure 3A shows the considered indexes obtained for LSLS (i.e., I_s and f_a), evaluated for all the analyzed churches and for the two considered observation periods. For V_R = 30 years, the SS Maria della Bruna church shows a low vulnerability, having both indexes exceeding the unity: I_s30 = 1.65 and f_a30 = 1.26. On the contrary, these indexes are significantly lower than unity if the observation period is equal to 75 years: in this case it results as I_s75 = 0.66, and f_a75 = 0.89. It is worth noting that SS Maria della Bruna church shows a low vulnerability if compared with others, except for San Pietro Caveoso church, where it is found I_s30 = 2.11 (I_s75 = 0.85) and f_a30 = 1.33 (f_a75 = 0.94). In Figure 3B, the resulting indexes for DLS are shown. In this case, the SS Maria della Bruna church provides, in all cases, numerical values that are always lower than unity: I_s30 = 0.92 (I_s75 = 0.37) and f_a30 = 0.95 (f_a75 = 0.59), even though higher than the other ones, except again for San Pietro Caveoso church where it is found that I_s30 = 1.01 (I_s75 = 0.40) and f_a30 = 1.00 (f_a75 = 0.62).

One may note that, as for the LSLS, the seismic safety index I_s constantly reduces by about 2.5 times in all the churches (about 1.4 times in the case of f_a), passing from V_R = 30 years to V_R = 75 years. This constant value may be justified by the fact that the ratio I_s30/I_s75 (or f_a30/f_a75) is dependent only on the return periods of the seismic actions expected on the site (accelerations expected on the site). The same conclusion may be done for the DLS, where I_s30/I_s75 results at about 2.5 (f_a30/f_a75 ≅ 1.6).

The results of Figure 3 clearly show that all the plotted curves for LSLS and DLS have qualitatively the same shape. In other words, the ranking of the seismic performances is always the same independently from the Limit State and from the index considered.

Seven failure mechanisms have been considered in this study: four directly connected to the ground, and the other three are not directly ground-connected. The considered mechanisms are the following:

For all the macro-elements, since the local masonry disintegration may be excluded, a monolithic behavior has been assumed. The macro-elements considered, as previously described, have been schematized starting from only visual investigations and geometric reliefs, since no clear cracking patterns have been encountered. In addition, no stabilizing contribution for each mechanism has been considered (such as interaction of transversal walls, or additional restrains). Therefore, the calculated activation multipliers underestimate the real ones. Furthermore, for simplicity, only the most vulnerable mechanisms have been considered in this study, not considering the in-plane response ones since they typically show higher activation multipliers.

Figure 4 illustrates the failure mechanisms in the case of SS Maria della Bruna church. It must be clarified that, for the comparison with the other churches, in the case of San Rocco and San Giovanni Battista churches, the triumphal arch mechanism has not been considered due to the absence of this architectonic element.

The spectral accelerations a0* (derived from α₀ according to the Equation 9), related to the activation of the considered macro-elements mechanism connected to the ground are shown in Table 4 and graphed in Figure 5A. Moreover, Table 4 reported the acceleration factors f_a for a V_R = 30 years and a V_R = 75 years, calculated as follows:

In Equation (19) for each church, it is considered the minimum values of a0* for a given limit state (i.e., the most vulnerable failure mode). Whereas, a_demand corresponds to the related expected seismic demand equal to: a_{g, LSLS30} S/q (or a_{g, LSLS75} S/q), or to a_{g, DLS30} S (a_{g, DLS75} S). Again, S has been assumed equal to 1 and q equal to 2.

The façade overturning results in the most vulnerable mechanism, with an activation acceleration of a0* = 0.078 g, while the others have acceleration significantly higher (i.e., transversal response of triumphal arch a0* = 0.121 g, about 55% higher than the value obtained for the façade overturning). It results that the most vulnerable mechanism for all churches is the simple overturning of the main façade, except for San Francesco d'Assisi church. In this case, the lowest spectral acceleration corresponds to the nave transversal response (a0* = 0.046 g). Furthermore, it should be remarked that in the case of San Giovanni Battista church, the accelerations obtained for all the mechanisms are significantly higher than the others, with a minimum value a0* of 0.162 g that is more than twice of the one obtained for the SS Maria della Bruna Church. Moreover, Table 5 and Figure 6A satisfy the seismic performance requirements in the case of V_R = 30 years for both LSLS and DLS, and in the case of V_R = 75 years for DLS only while for LSLS f_a75it is slightly < 1 (I.e. f_a75 = 0.975). Again, as observed for the LV1 method, the Cathedral is close to satisfying all the required seismic safety obtaining indexes, and this satisfaction level is similar or higher than the other considered churches, except for San Giovanni Battista.

In Table 5 and Figure 5B the accelerations relative to the macro-elements that are not directly ground-connected are shown. For SS Maria della Bruna church, the most vulnerable mechanism is the gable breakout with an activation acceleration of 0.131 g. By comparing the results with the others, it should be noted that only San Francesco d'Assisi and San Giovanni Battista churches have low a0*, and therefore, the Cathedral is the third most vulnerable regarding these mechanisms. The seismic demand a_demand is calculated in these cases as S_e(_T1)30ψ(Z)γ/q for V_R = 30 years [or S_e(_T1)75ψ(Z)γ/q for V_R = 75 years] for the LSLS, and as S_e(_T1)30ψ(Z)γV_R = 30 years [or S_e(_T1)75ψ(Z)γ for V_R = 75 years] for the DLS. Given this seismic demand, the acceleration factors can be calculated as the ratio of a0* and such demands (i.e., f_a30 = a0*/S_e(_T1)30ψ(Z)γ/q and f_a75 = a0*/S_e(_T1)75ψ(Z)γ/q). By analyzing these acceleration factors (Figure 6B), it is worthy of noting that the SS Maria della Bruna church has values that are close or beyond the unity (i.e., for LSLS f_a30 = 1.074, f_a75 = 0.736; while for DLS f_a30 = 2.382, f_a75 = 1.284), confirming the need for only light interventions in order to obtain the required seismic performances. Furthermore, it should be remarked that the most vulnerable church regarding these mechanisms is San Giovanni Battista (i.e., a0* = 0.064 g top façade overturning), which is less vulnerable regarding the mechanisms directly connected to the ground.

The comparisons of the results calculated with the LV1 and LV2 methods of the Italian directive (G.U. n. 47, 26/02/2011) are illustrated in Figure 7. It should be specified that, in the case of the LV2 method, only the ground connected mechanisms are considered in these comparisons. Precisely, Figures 7A,B show the results for LSLS by referring to V_R = 30 years and V_R = 75 years, respectively. It is worth to note that for the Matera Cathedral the LV1 method, which is simpler and less accurate, overestimates the results if compared with the LV2 one. This trend is also found for San Pietro Caveoso church. However, for SS Maria della Bruna the scatter is low and, at least for V_R = 30 years, the ratios are always beyond the unity. Whereas, in Figures 7C,D, the comparisons between the two adopted methods for the considered DLS are shown. Contrarily to the LSLS, in this case the LV1 method is, for the considered churches, always conservative with respect to the LV2 one. With the LV1 method only, the San Pietro Caveoso church exhibits, for V_R = 30 years, a seismic capacity closer to the demand. As for the LV2 method, the capacity of all churches is always higher than demand except for San Rocco and San Francesco d'Assisi for V_R = 75 years. Anyway, it should be remarked that the vulnerability classification provided by considering the two methods separately would be very similar since, as illustrated in the Figure 7, the trends found are quite the same. Nevertheless, by increasing the detail level of the analysis (from LV1 to LV2), it is found that in the analyzed cases, the LV1 overestimates the seismic performance, as for SS Maria della Bruna and San Pietro Caveoso church. However, it should be kept in mind that in this study, for simplicity, the LV2 results have been obtained by neglecting, as previously described, all the stabilizing contributions. Therefore, the actual results will be higher than those discussed here, and higher than the ones provided by the LV1 method. This confirms that after a first and fast numerical evaluation, useful for classifying the case studies, there is always the need to implement more realistic numerical models that cannot be generalized since they are related to the boundary conditions of the analyzed problem.

Finally, in this study a comparison with the method proposed in (Dìaz, 2016) is illustrated and then validated in D'Amato et al. (2018). Its application is numerically reported in Table 6. With this method, the SS Maria della Bruna church obtains a score of R = 22.15, resulting as the lowest value. All the obtained results are compared in relative terms in Figure 8: at first, each seismic score (R_i) is divided by the minimum value found (R_min); then, all resulting ratios (R_i/R_min) are represented in ascending order. For comparison, in Figure 8 the values calculated with the LV1 method are also reported. In detail, the dimensionless return period T_{i, LSLS}/T_LSLSmin, where T_{i, LSLS} is the returning period corresponding to the LSLS achievement and T_LSLSmin is the minimum value found among the churches (corresponding to 181 years for the San Francesco d'Assisi church), are plotted. In addition, in Figure 8 the dimensionless ground acceleration capacity a_{i, LSLS}/a_LSLSmin is also reported, where a_{i, LSLS} is the ground acceleration related to the LSLS achievement and a_LSLSmin is the lowest value found (that is 0.127 g corresponding, again, to the San Francesco d'Assisi church). One may consider these ratios as relative seismic vulnerability indexes: the ratio T_{i, LSLS}/T_LSLSmin (or a_{i, LSLS}/a_LSLSmin) increases as the relative seismic vulnerability (R_i/R_min) decreases, where both indexes are calculated with respect to the most vulnerable church (in this case San Francesco d'Assisi church). The results found confirm the applicability of the new simplified method proposed in Dìaz (2016) and validated in D'Amato et al. (2018). This method may be intended, in a multi-level frame-work approach, as a “Level of Evaluation 0” (LV0) since it permits rapidly ranking the seismic performance at a territorial level. In this way, useful information for individuating the priorities may be found, to be investigated in more detail with more refined approaches (such as LV1, LV2, or LV3).

In this work, it has been conducted as a comparative analysis among the seismic performances of ancient masonry churches, by the means of different simplified methods. At first, the seismic vulnerability of the Matera Cathedral, called SS Maria della Bruna of Matera, has been investigated, and then the obtained results have been compared with the ones of four churches located in the Matera city center. The great advantage of these methods is that they don't imply advanced structural analyses and investigations. Thus, they can be used as decision making tools for identifying the priorities, and consequently, for performing further analyses and designing the potential interventions.

The analyses performed with the methods, indicated as LV1 and LV2 within the Italian directive (G.U. n. 47, 26/02/2011), have highlighted that, among the churches analyzed, the Matera Cathedral is one of the less vulnerable. Specifically, according to the LV2 method, it satisfies the seismic protection level required in the case of V_R = 30 years, both for LSLS and for DLS. On the contrary, for V_R = 75 years, few and light interventions are requested in order to achieve the required seismic demand for LSLS.

Finally, the comparisons with the new simplified method for seismic risk assessment confirm that it may be considered as a preliminary appraisal method for comparing the seismic performances of ancient churches at a territorial level. It may be proposed as a LV0 approach, since it requires limited and qualitative information for ranking the seismic performance at a territorial scale. However, by analyzing the obtained results, it has emerged that, in general, simpler methods may overestimate the actual seismic performance of a church. Therefore, the simplified methods, although useful for comparing and ranking the churches seismic performances at a territorial level, cannot substitute the refined ones for realistically assessing the seismic behavior of a structure.

All data generated or analyzed during this study are included in this published article.

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.