Simone Pezzuto
email: simone.pezzuto@usi.ch

Ali Gharaviri

Ulrich Schotten

Mark Potse

Giulio Conte

Maria Luce Caputo

François Regoli

Rolf Krause

Angelo Auricchio

Caputo Ml

wave morphological variability in patients with paroxysmal atrial fibrillation: an in silico study

Keywords: P-wave variability, P-wave morphology, Paroxysmal atrial fibrillation, inhomogeneous conduction, fibrosis, cosine distance

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Pulmonary vein (PV) isolation is the treatment of choice in patients with paroxysmal and persistent atrial fibrillation [START_REF] Haïssaguerre | Spontaneous Initiation of Atrial Fibrillation by Ectopic Beats Originating in the Pulmonary Veins[END_REF] . However, recurrence rate is still significant especially at mid-and long-term follow-up [START_REF] Verma | Approaches to Catheter Ablation for Persistent Atrial Fibrillation[END_REF] . Apart from reconnection of PV(s), alternative mechanism of atrial fibrillation recurrence may include persisting presence and extension of heterogeneous conduction areas (the substrate for reentry) and/or non-PV dependent triggered activity [START_REF] Zahid | Patient-derived models link re-entrant driver localization in atrial fibrillation to fibrosis spatial pattern[END_REF] . Areas of abnormal, heterogeneously distributed conduction within the atria have been well described in pre-clinical [START_REF] Akoum | Atrial fibrosis helps select the appropriate patient and strategy in catheter ablation of atrial fibrillation: A DE-MRI guided approach[END_REF][START_REF] Marrouche | Association of Atrial Tissue Fibrosis Identified by Delayed Enhancement MRI and Atrial Fibrillation Catheter Ablation[END_REF] and clinical setting [START_REF] Fukumoto | Association of Left Atrial Local Conduction Velocity With Late Gadolinium Enhancement on Cardiac Magnetic Resonance in Patients With Atrial Fibrillation[END_REF] , and can possibly represent the substrate for prolonged atrial conduction time (P-wave duration). Previous studies have indicated that both abnormal Pwave axis and long P-wave duration are strongly associated with an increased risk of AF [START_REF] German | Atrial Fibrillation Predictors: Importance of the Electrocardiogram[END_REF] . More recently, dynamic morphological alteration of P wave (beat-to-beat P-wave morphology variability) has been reported by our group [START_REF] Conte | Usefulness of P-Wave Duration and Morphologic Variability to Identify Patients Prone to Paroxysmal Atrial Fibrillation[END_REF] and others [START_REF] Filos | Beat-to-beat P-wave morphology as a predictor of paroxysmal atrial fibrillation[END_REF][START_REF] Huo | Variability of P-wave morphology predicts the outcome of circumferential pulmonary vein isolation in patients with recurrent atrial fibrillation[END_REF] , and the combination of abnormal P wave duration and beat-to-beat variability has been shown to more effectively discriminate patients with history of AF from healthy age-matched subjects. This finding however has not been mechanistically well explained. P-wave morphological variability, consistent with a beat-to-beat timeframe, could be related to a shift of the earliest activation site (EAS) in sinoatrial node (SAN) region [START_REF] Loewe | Influence of the earliest right atrial activation site and its proximity to interatrial connections on P-wave morphology[END_REF] . For instance, a shift in the EAS can be induced by exercise or drugs [START_REF] Löbe | Sinus node modification with an ultra high-density electroanatomical mapping system in inappropriate sinus tachycardia[END_REF] . SAN remodeling has been reported in patients with congestive heart failure [START_REF] Sanders | Remodeling of Sinus Node Function in Patients With Congestive Heart Failure: Reduction in Sinus Node Reserve[END_REF] . Therefore, a possible mechanism explaining the difference in P-wave morphological variability between healthy subjects and patients with history of AF could be based on the interaction between EAS variability, physiologically always present, and the heterogeneity in conduction, generally increased in patient with history of AF.

By using a generic computer model of the human atria and torso, we aimed at reproducing the beat-to-beat variation of P-wave and the duration of P wave as clinically observed. We simulated P-waves as recorded on the 12-lead ECG, considering variation in the exit zone of the SAN, and modeling distribution and size of area of conduction heterogeneities generated by a novel automatic method in which spatial correlation, distribution, and size was prescribed.

Methods

Reference atrial model

We considered a reference anatomical model of the human atria previously built using MRI data from a normal subject [START_REF] Potse | P-wave complexity in normal subjects and computer models[END_REF] . The model contained several anatomical features such as the Bachmann's bundle, the crista terminalis, the interatrial bundles, and the pectinate muscles. All were designed by following anatomical studies [START_REF] Ho | Atrial structure and fibres: morphologic bases of atrial conduction[END_REF] . Orientation bundles were manually added to define the fiber direction in the atrial tissue. A computational 3d structured grid with uniform spacing of 0.2mm was obtained from the anatomical model, with roughly 5 million voxels in total. Each voxel was automatically marked with the corresponding anatomical region for the electrophysiological model parametrization (see Figure 1B).

Model of slow conducting regions

The slow conducting area in the model was defined by thresholding an automatically generated random field with given spatial correlation. More specifically, the random field was the average of two random Gaussian fields with different spatial correlation lengths, respectively set to 2cm and 4mm, thus to model organ-level and tissue-level inhomogeneity. The organ-level random field was sampled by using the truncated Karhuhen-Loève expansion for a squared-exponential covariance function, with geodesic distance computed with the eikonal model [START_REF] Quaglino | Fast uncertainty quantification of activation sequences in patient-specific cardiac electrophysiology meeting clinical time constraints[END_REF] . The tissue-level random field was sampled by solving the linear fractional stochastic partial differential equation

( 1 𝜌 2 𝐼 -Δ) 𝛼/2
𝑢 = 𝑊 , where 𝑊 was white-noise, 𝛼 = 𝜈 + 3 2 with 𝜈 related to the spatial regularity of the field, and 𝜌 > 0 related to the spatial correlation length. In practice, the partial differential equation smoothed the white-noise by introducing some spatial correlation of the order of 𝜌 scale [START_REF] Lindgren | An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach[END_REF][START_REF] Stein | Interpolation of spatial data: some theory for kriging[END_REF] . The model was iteratively solved for increasing values of 𝜈 = 4𝑘-3 2 , 𝑘 = 1,2, …, and 𝜌 = 𝜃 2√𝜈 ⁄ . In the limit case, for 𝑘 → ∞, the random field has squared-exponential spatial correlation of parameter 𝜃. The superposition of the random fields was eventually rescaled from 0 to 1. An illustrative example of such random field is provided in Figure 2A.

The area of slow conduction was finally defined as the region of the domain such that the random field was less than a threshold 𝑇. The threshold was selected by a bisection argument to ensure that the total volume of the slow conducting region was equal to a given fraction 𝑓 of the total atrial volume (see Figure 2B).

Bachmann's bundle, interatrial bundles and a portion of the RA (the intercaval region) were excluded from the threshold procedure but included in the random field generation.

Eikonal model for atrial activation

The electric activation of the atria was modeled by using the eikonal model [START_REF] Franzone | Spread of excitation in 3-D models of the anisotropic cardiac tissue. II. Effects of fiber architecture and ventricular geometry[END_REF] with no curvature-correction:

𝜃 √𝛽 √𝑮 m ∇𝜏 ⋅ ∇𝜏 = 1,
where 𝜏 = 𝜏(𝑥) was the activation time at some location 𝑥 of the atria, 𝜃 was the conduction velocity scaling factor, 𝛽 was the surface-to-volume ratio and 𝑮 m was the monodomain electric conductivity tensor. The activation was initiated by setting 𝜏 to zero at a unique EAS. The electric conductivity tensor was transversely isotropic with major axis in the local fiber direction. The longitudinal and transverse intra-and extra-cellular electric conductivity were modeled with regional heterogeneity. The model contained 4 of such regions: an isotropic intercaval area in the right atrium (0.75 mS/cm monodomain conductivity for both longitudinal and transversal); the Bachmann's bundle and atrial connections (4.5 mS/cm and 0.72 mS/cm respectively longitudinally and transversely); tissue with oriented bundles (1.5 mS/cm and 0.72 mS/cm); and a slow conducting area (1 mS/cm and 0.28 mS/cm). The surface-to-volume ratio was set to 800 cm -1 everywhere except in the slow conducting region, where it was doubled. The conduction velocity scaling factor (units cm ⋅ ms -1 ⋅ mS -1/2 ) and the transmembrane potential template were pre-computed with a single fiber monodomain simulation using the Courtemanche ionic model [START_REF] Courtemanche | Ionic mechanisms underlying human atrial action potential properties: insights from a mathematical model[END_REF] . A summary of the parameters adopted in the simulation is found in Table 1.

P-wave simulation

The 12-lead surface ECG was simulated by using the formula

𝑉(𝑡) = ∫ 𝑮 𝒊 ∇𝑉 𝑚 ⋅ ∇𝑍 𝑑𝑥, Ω
where Ω was the atrial tissue, 𝑉 𝑚 (𝑥, 𝑡) was the transmembrane potential, 𝑮 𝑖 was the intracellular electric conductivity, and 𝑍(𝑥) was a precomputed lead field. For each lead, the bidomain problem in the full torso was solved for the lead field 𝑍 = 𝑍(𝑥) by injecting a unit current at the electrode locations. The torso, heterogeneous and anisotropic, included atrial and ventricular tissue, blood cavities, lungs, skeletal muscles.

The transmembrane potential was set to 𝑉 𝑚 (𝑥, 𝑡) = 𝑈(𝑡 -𝜏(𝑥)), with 𝑈 a pre-computed action potential

template. An efficient implementation of the ECG computation in combination with the eikonal solution, on GPGPU architecture, was provided by the software Propeiko 21 .

Randomized early activation site location

The electric activation was triggered at a single EAS location for each beat, located in the sinoatrial node (SAN) area, on the right atrium (RA). From an anatomical viewpoint, the SAN is functionally isolated for the atrial myocardium, except for several sinoatrial connection pathways (SACP) uniformly distributed around the SAN. A total of 12 SAN locations were considered in the atrial model, as depicted in Figure 1A. These were manually placed on the RA starting from the superior vena cava (SVC) and moving towards the inferior vena cava (IVC). The selected locations were consistent with histological studies [START_REF] Csepe | Novel application of 3D contrast-enhanced CMR to define fibrotic structure of the human sinoatrial node in vivo[END_REF] and previous in-silico investigations [START_REF] Loewe | Influence of the earliest right atrial activation site and its proximity to interatrial connections on P-wave morphology[END_REF] .

To enable a beat-to-beat variability of the activation, the EAS location was modeled as a uniform random variable of spherical shape and intersected with the right-atrial tissue. The sphere was centered at a specific SAN locations with a predefined radius. The sampling radius was set from 1mm up to 5mm.

P-wave analysis

For each sampled EAS location, the corresponding P-wave was computed. The P-wave morphological variability was evaluated with two approaches: one based on the spectral radius of the cosine distance matrix, and the second based on the standard deviation of the P-waves. Denoting by 𝑉 𝑖 (𝑡) the 𝑖-th P-wave, 𝑖 = 1, … , 𝑁, the cosine distance between the 𝑖-th and 𝑗-th P-wave was as follows

𝑑 cos (𝑉 𝑖 , 𝑉 𝑗 ) = 1 𝜋 cos -1 ( ∫ 𝑉 𝑖 (𝑡)𝑉 𝑗 (𝑡) 𝑑𝑡 𝑇 0 (∫ 𝑉 𝑖 2 (𝑡) 𝑑𝑡 𝑇 0 ) 1/2 (∫ 𝑉 𝑗 2 (𝑡) 𝑑𝑡 𝑇 0 ) 1/2 ) .
Then, the variability 𝑣 cos was evaluated as the largest absolute value of the eigenvalues of the matrix 𝐷 𝑖𝑗 = 𝑑 cos (𝑉 𝑖 , 𝑉 𝑗 ), i.e. its spectral radius: max 𝑖=1,…,𝑁 |𝜆 𝑖 (𝐷)|, divided by the number of eigenvalues or, equivalently, the number of beats.

With the second methodology, the variability was assessed from the 𝑝-norm (with 𝑝 = 1 or 𝑝 = ∞) of the standard deviation of the P-wave sample

𝑣 𝑝 = ‖ √ 1 𝑁 -1 ∑ (𝑉 𝑖 (𝑡) -𝑉 ̅ (𝑡)) 2 𝑁 𝑖=1 ‖ 𝐿 𝑝 (𝑂,𝑇) ‖𝑉 ̅ ‖ 𝐿 𝑝 (𝑂,𝑇)
,

where 𝑉 ̅ (𝑡) = 1 𝑁 ∑ 𝑉 𝑖 (𝑡) 𝑁 𝑖=1
was the average P-wave. Both measures of variability were non-dimensional.

Monodomain simulations

An activation map and corresponding surface ECG was computed with the monodomain model and compared to the eikonal model. The monodomain solver was Propag [START_REF] Krause | Hybrid Parallelization of a Large-Scale Heart Model[END_REF] , using the same computational grid and the same parameters, when applicable. The comparison was performed in the control case and with heterogeneous conduction, pacing from a fixed EAS location (SAN01). In the monodomain model, a current of 150A/cm 3 was applied for 2ms in a 1mm 3 volume centered at the EAS location. In the eikonal model, we simply set the activation time to the onset value at the EAS location. The P-wave was computed with the lead field approach for the eikonal solution, and with the bidomain model in the torso for the monodomain solution. The torso geometry for the bidomain model was as for the lead field computation. Sampling time was 0.2ms in both cases (time-step for explicit time integration in the monodomain model was 0.01ms).

Results

Comparison between the eikonal and the monodomain model

An activation map and the corresponding P-wave was computed by using either the eikonal or the monodomain model, both in control case and with slow conducting areas (30% area, pattern 3). In the eikonal model, the parameter 𝜃 was globally adjusted to reflect the average conduction velocity observed in the monodomain case (𝜃 = 1.519 in control case and 𝜃 = 1.4273 for heterogeneous conduction). A fixed delay in activation of 2ms was introduced in the eikonal model to accommodate the delay in activation in the monodomain case, due to the current injection. Results are presented in Figure 3.

In control case, the correlation of between the activation maps of the two models was 0.999, the RMSE was 1.2ms and the 99 th percentile of the absolute error was 2.9ms. In patchy conduction case, the correlation, the RMSE and the 99 th percentile of the error were respectively 0.998, 2.3ms and 6.4ms. By using the reference parameter 𝜃 = 1.6176, obtained with a high resolution 1d monodomain simulation, the rootmean-square error (RMSE) in control case was increased to 5.1ms. In all cases, the monodomain activation resulted slightly longer than the eikonal one.

The cosine distance between the P-waves in control case was between 0.02 (lead V6) and 0.08 (lead V3). In lead II (maximum amplitude lead), the cosine distance was 0.02. For the heterogeneous conduction case, the cosine distance was 0.016 in lead II, and ranged from 0.015 (aVL) to 0.093 (V3). In terms of amplitude and duration, P-waves were indistinguishable. The RMSE was 5.4 ⋅ 10 -3 mV on average in both cases. The monodomain solution took 100s on CSCS Piz Daint cluster using 4 nodes. The eikonal solution took 20sec in total on a single GPGPU (Nvidia GTX1080).

P-waves simulation

A total of 60 beats for each of the 12 SAN exit locations were simulated for the control case (no slow conduction region) and for 15 different heterogeneous conduction patterns. For each random pattern case (see Figure 2C-E, on the right), the fraction of slow conducting tissue was set from 20% up to 60%, with increments of 10% (Figure 2B). In total, 16 configurations were simulated (including the control case) for a total of 11'520 P-waves. The EAS locations, sampled around each SAN, were different from case to case.

Simulations were performed with the same setup as above (20sec each on a local workstation).

Variability in control and with heterogeneous conduction

In control case, the average P-wave amplitude and duration in lead II were 0.210.07 mV and 140.25.2 ms, respectively. The P-wave variability using the cosine distance was 0.0750.043. Variability was comparable among SAN exit locations, except for EASs around SAN11 and SAN12, for which the variability was 2-fold larger. On average, lead V3 reported the largest variability (0.1240.036), while V6 the smallest one (0.0390.013). The 𝑣 1 and the 𝑣 ∞ distance were respectively 0.200.137 and 0.250.138. A similar trend to cosine distance was observed.

In Figure 4 we report a comparison between the control and case 3 with 40% area of heterogeneous conduction. A similar analysis was performed for all the other cases. Overall, the variability was larger in the presence of slow conducting regions, of a factor of 2 on average. In Figure 4A, we report the distribution of the entries of the cosine distance matrix, that is the cosine distance between all possible pairs of simulated P-waves. The histogram showed a significantly increased spread of cosine distance in the heterogeneous case, in all the leads. The distribution of distances was bell-shaped (as in control) but shifted to the right and with significantly fatter tails. The increase in variability was also confirmed by the box-plot analysis (Figure 4B).

In Figure 4C and D we superimpose the P-waves to highlight the variability. The average P-wave in the control case was very close (in amplitude and morphology) to the P-wave produced by pacing exactly at the SAN exit location. In the presence of slow conducting regions, the average P-wave was smoother and with lower amplitude than the P-wave obtained as above.

Finally, in Figure 5A, we report the cosine distance between all P-wave pairs, as a function of the (Euclidean) distance between the corresponding EASs. In the heterogeneous case, the variability versusEAS distance relationship was steeper and less sharp than the control case.

Effect of distribution and area of inhomogeneous conduction

The effect of slow conduction distribution on the P-wave variability was explored by increasing the fraction of slow conducting region for 3 different patterns. Results are reported in Figure 5B and Figure 6. In general, an increasing area of slow conduction led to an increase of variability. For some pattern (e.g. pattern 1 and pacing around SAN02, or pattern 2 with pacing around SAN03), the variability 𝑣 cos was maximum for a fraction of 20% and reducing for larger values. The cosine variability was also significantly affected by the distribution of the slow conducting region: pacing around SAN03, variability in lead II was 2-fold the control case in case 3 with 40% area, while it was 1.5-fold in case 2 and comparable to control in case 1. Only very few cases reported a reduced variability with respect to control case. This occurred when activation around SAN12, especially for case 3.

Effect of filtering on variability

The analysis was repeated with filtering applied to all the P-wave samples prior the evaluation of the variability. The filter was a 5 th order Butterworth filter with cut-off frequency of 40Hz. With filtering, the cosine variability was 33% lower in control case, and from 16% to 46% lower in the inhomogeneous cases. In general, the cosine variability with filter applied was still larger than the control case, corroborating the above analysis. More specifically, the ratio in variability (heterogeneous over control) was higher with filter applied (from 1.3x to 2.8x versus 1.1x to 2.6x in absence of filtering).

Discussion

This manuscript presents several novelties. First, we developed a novel approach to generate areas of conduction heterogeneities in the atria with a given spatial correlation, distribution, and size. Secondly, we could simulate a beat-to-beat P-wave variability as reported by clinical observation in healthy subjects and in patients with paroxysmal AF [START_REF] Conte | Usefulness of P-Wave Duration and Morphologic Variability to Identify Patients Prone to Paroxysmal Atrial Fibrillation[END_REF] , by adopting a physiologically-motivated variation of SAN exit location which has not been previously accounted for in the literature, to the best of our knowledge. Finally, we compared the eikonal model with the lead field approach for P-wave computation against the monodomain model.

Our model suggested that the beat-to-beat variability is related to presence and quantity of heterogeneous conduction, in combination with the intrinsic variability in the EAS location in the SAN region [START_REF] Li | Redundant and diverse intranodal pacemakers and conduction pathways protect the human sinoatrial node from failure[END_REF] . The numerical investigation showed, as expected, that the variability increases with respect to the EAS variability, in both control case and with heterogeneous conduction. In the latter case, however, the variability grows more rapidly than the former, as certified by the results in Figure 5. We also observed that for the sampling radius of the EASs had a modest effect on variability for values larger than 4 mm.

The EAS physiological variability alone can only partially account for the difference in P-wave morphological variability between healthy subjects and patients with history of AF. We found that the P-wave variability is also a complex function of distribution and area of slow conducting region. Moreover, its magnitude also depends from which lead is measured, and whether a filtering procedure was performed prior the analysis.

In addition, there is no clear definition of variability, mathematically speaking. Despite these difficulties, we have found a consistent increase of variability in the presence of slow conducting regions, compared to the control case, except in a very few cases. For instance, the variability in the presence of slow conducting regions was comparable or slightly lower than the control case when activation occurred around SAN11 and SAN12, which are close to IVC. Interestingly, the variability in control was already very large, probably because these locations are close to the interatrial bundles (IAB). Moreover, none of the analyzed patterns had an extensive slow conducting area in the region of these SAN exit locations.

The presence of slow conducting region around the SAN exit location only partially explains the increase in P-wave variability. In the experiments, another important factor was the presence of such region close to the IABs and the Bachmann's bundle (BB), and its distribution on the left atrium. In such cases, the P-wave variability could be explained by a shift of some portions of the P-wave, and an increased fragmentation. In fact, the slow conducting region was present in all the patterns around SAN01 to SAN03 and the BB.

Another interesting observation is that the P-wave variability increased with the fraction of slow conducting region but, in some cases, for very large values of fraction the variability declined (see Figure 5B, case 2). It is reasonable to argue that under such circumstances the variability is comparable to the control case, because heterogeneity at organ-scale is reduced, albeit heterogeneity at tissue-scale may still be large. This may also indicate that intermediate fractions of slow conducting region (or severe heterogeneity at organ-scale) might characterize the electric substrate of patients with paroxysmal AF. In this respect, the P-wave morphological variability in sinus rhythm may provide additional insight on the stage of AF, from paroxysmal to persistent [START_REF] Schotten | Pathophysiological Mechanisms of Atrial Fibrillation: A Translational Appraisal[END_REF] .

Importance of heterogenous conduction

Heterogeneous conduction in the atria has been already reported and explored by other authors in the literature, especially as a possible mechanism to sustain and enhance AF 3 . These regions have been identified

as fibrotic, and they are generally obtained from imaging techniques [START_REF] Marrouche | Association of Atrial Tissue Fibrosis Identified by Delayed Enhancement MRI and Atrial Fibrillation Catheter Ablation[END_REF] , in patient-specific context. Our novel efficiency: for the atrial model adopted in this work, a full P-wave simulation took less than half a minute on a single GPGPU. Incidentally, the total cost could have been reduced to 5sec by using a coarser computational grid, but in this case the comparison to the reference monodomain solution would have been less favorable (not shown). In a previous study from our group we compared the eikonal-based method against the full monodomain model, in terms of activation map and ECG, for the ventricular activation, obtaining very good agreement [START_REF] Pezzuto | Evaluation of a Rapid Anisotropic Model for ECG Simulation[END_REF] . In this work we extended the analysis to the atrial activation and the P-wave in sinus rhythm (SR), confirming the appropriateness of the methodology for both the control case and with patchy slow conduction. The suitability of the eikonal model was also explored previously by other group with convincing results [START_REF] Jacquemet | An Eikonal Approach for the Initiation of Reentrant Cardiac Propagation in Reaction-Diffusion Models[END_REF][START_REF] Krueger | Towards personalized clinical in-silico modeling of atrial anatomy and electrophysiology[END_REF] , albeit they focused only on the activation map.

The parameter 𝜃 in the eikonal simulations required an adjustment for matching the corresponding monodomain simulation. The reason is the relatively large numerical error in the monodomain case. Using a standard finite difference discretization at 0.2mm voxel size, the conduction velocity is generally 10% to 20% slower (depending on the wave-front thickness) than the converged value [START_REF] Pezzuto | Space-discretization error analysis and stabilization schemes for conduction velocity in cardiac electrophysiology[END_REF] . This is consistent with 𝜃 values we found. In this respect, we believe that the eikonal model was more accurate (in terms of conduction velocity) than the monodomain model for the presented simulations. However, if the slow conduction region were modeled by removing computational cells from the active tissue, or by substituting a fraction of cells with fibroblasts [START_REF] Sridhar | Effect of myocyte-fibroblast coupling on the onset of pathological dynamics in a model of ventricular tissue[END_REF] , the classic eikonal approach would have not been suitable.

The rapidity of eikonal model for atrial simulations potentially enables multi-fidelity parameter identification and uncertainty quantification in SR [START_REF] Quaglino | Fast uncertainty quantification of activation sequences in patient-specific cardiac electrophysiology meeting clinical time constraints[END_REF][START_REF] Konukoglu | Efficient probabilistic model personalization integrating uncertainty on data and parameters: Application to Eikonal-Diffusion models in cardiac electrophysiology[END_REF] . The personalization process of the atrial model is of paramount importance for clinical application, and this can be performed in SR.

Nonetheless, the eikonal model comes with significant approximations with respect to the monodomain or bidomain model. In its original form, it is not suitable for atrial fibrillation (AF) simulations, lacking reexcitability of the tissue, thus limiting its clinical applicability. However, modifications of the original model have been proposed to include action potential refractoriness [START_REF] Jacquemet | An Eikonal Approach for the Initiation of Reentrant Cardiac Propagation in Reaction-Diffusion Models[END_REF] and, more recently, conduction velocity dependency on the cycle length [START_REF] Corrado | A conduction velocity adapted eikonal model for electrophysiology problems with re-excitability evaluation[END_REF] . Finally, the lead field approach is fully equivalent to simulating the P-wave using the bidomain model in the torso with a precomputed transmembrane potential. It includes both heterogeneity and anisotropy of the torso and heart, and it is a computationally efficient method when the extracellular potential is required only at a few given locations, such as the standard 12-lead ECG.

Limitations

This study certainly comes with some limitations. The effect of the distribution of heterogenous conduction was analyzed for only 3 patterns. No distinction was made between LA and RA in the generation of the random field, although a difference in area and distribution between the two chambers might be present. In the model, we assumed a sharp reduction in the intracellular conductivity within the slow conducting regions, thus to produce a reduced conduction velocity. The assumption was loosely based on experimental data. Some groups [START_REF] Zahid | Patient-derived models link re-entrant driver localization in atrial fibrillation to fibrosis spatial pattern[END_REF] have proposed a similar approach to model fibrotic regions, but there is no consensus on how to select the parameters.

Conclusions

In this work, we showed that the morphological P-wave variability could be explained by an intrinsic variability in the SAN exit points in combination with a diffuse and heterogenous conduction in the atrial tissue. The larger the slow conducting area is, the larger the variability. The heterogenous conduction may be associated to fibrotic regions.

By using the same anatomy and slow conduction area in the monodomain AF simulations, we observed an increase in the variability of AF conduction pattern. An increase in the slow conducting area led to an increase in AF conduction pattern complexity. Hence, this might be an important step forward to help the outcome of AF computer models become closer to clinical and experimental recordings.

Paroxysmal AF has been traditionally assumed to originate from a combination of triggered activity and reentry circuits mediated by the atrial refractory period and the conduction velocity. Computer modeling is a recognized tool of utmost importance to provide new mechanistic insights of arrhythmias behavior and to identify the coexistence of different mechanisms responsible for the arrhythmia development. Indeed, for paroxysmal AF, the combination of fibrosis and variability in the SAN exit points might provide information on the arrhythmia stage, which is important while assessing pulmonary vein isolation outcomes. Moreover, P-wave beat-to-beat variability can be used to identify patients with asymptomatic episodes of paroxysmal AF and can drive the subsequent diagnostic and therapeutic management.

Further computational studies and patient-specific approaches are needed to address prediction of optimal pharmacological drug agents and ablation sites able to decrease P-wave beat-to-beat morphological variability and avoid episodes of paroxysmal AF.

FiguresFigure 1 (

 1 Figures

Figure 2

 2 Figure 2 Generation of slow conduction regions (SCRs) for various patterns and densities.. On the top left, the background random field for pattern 1 is represented, with low values towards blue and high values towards the red. On bottom left, the result of the thresholding procedure to extract the SCR at different level of total area (20%, 40% and 60%).. On the right, all the patterns considered in this study, represented only for 30% of SCR. . SCR is in bright green in all figures.

Figure 4

 4 Figure 4 Comparison between P-wave beat-to-beat variability in the control case (no slow conduction regions) and in the presence of regions of heterogeneous conduction(pattern 3, 40% of slow conduction area). (A) Histogram of the entries of the cosine distance matrix in control (blue) and heterogeneous case (orange), for each ECG lead. (B) Box-plot of the cosine distance in control andheterogeneous case, for each ECG lead. (C) Lead II average P-wave with shaded  2 standard deviation and (D) superimposed P-waves, with EAS around SAN01 an variability of 4 mm.

Figure 5 P

 5 Figure 5 P-wave in lead II variability as a function of the EAS radius sampling. (A) Scatter plot between the EAS pairwise Euclidean distance, and the P-wave cosine distance, for both control and heterogeneous case. (B) P-wave variability as a function of the fraction of the slow conduction region, for 3 different patterns.

Acknowledgments

The authors acknowledge financial support by the Theo Rossi di Montelera Foundation, the Metis Foundation Sergio Mantegazza, the Fidinam Foundation, and the Horten Foundation to the Center for Computational Medicine in Cardiology. This work was also supported by grants from the Swiss National Supercomputing Centre (CSCS) under project ID s778.

Condensed abstract

By using a generic computer model of the human atria, the P-wave variability on a beat-to-beat basis was explored. The variability increased in the presence of heterogeneous conduction, corroborating clinical observations.

What's new

• Simulation of beat-to-beat P-wave variability with physiological earliest activation site variability in the sinoatrial node region.

• Automatic generation of multi-scale, spatially-correlated heterogeneous conduction regions for realistic 3d atrial anatomy.

• In the presence of heterogeneous conduction in the atria, P-wave variability increases.

• P-wave and activation map compares excellently between the monodomain solution and the eikonal model with lead field approach.

• Robust methodology to measure the P-wave morphological variability. methodology significantly expands the previous knowledge as we proposed a general methodology to automatically generate such fibrotic regions, by requiring a specific spatial correlation and total fibrotic area.

The generated patterns closely resemble those available in the literature [START_REF] Marrouche | Association of Atrial Tissue Fibrosis Identified by Delayed Enhancement MRI and Atrial Fibrillation Catheter Ablation[END_REF][START_REF] Fukumoto | Association of Left Atrial Local Conduction Velocity With Late Gadolinium Enhancement on Cardiac Magnetic Resonance in Patients With Atrial Fibrillation[END_REF] . The key aspect of our approach was the enforcement of a spatial correlation at multiple scales. It is reasonable to argue that the formation of fibrosis is spatially-correlated at the tissue-level, due to cellular coupling. An organ-scale spatial correlation may instead be favored by mechanical coupling, for instance. Unfortunately, to the best of our knowledge, spatial correlation of such fibrotic regions has never been quantitatively investigated in the literature.

Variability measure

Conte et al. [START_REF] Conte | Usefulness of P-Wave Duration and Morphologic Variability to Identify Patients Prone to Paroxysmal Atrial Fibrillation[END_REF] suggested the use of the Euclidean distance between one beat and the successive one, normalized with respect to the current beat. An overall assessment of variability was eventually based on statistics of the Euclidean distance time series. While such measure may be practical in a clinical context, it comes with some limitations: it is not a distance in strictly mathematical terms (it misses symmetry), and it is not independent on the order of the beats. In this work, we proposed two different methods to compute the P-wave variability which they both overcame the mentioned limitations. Interestingly, we found that the 𝑣 cos variability measure is often very close (relative difference less than 5%) to the mean value of the entries of the distance matrix, which is straightforward to compute. This is a consequence of the positivity of the matrix (all entries are non-negative), and the Frobenius-Perron Theorem: the dominant eigenvalue, which is non-negative and corresponds to the spectral radius, is bounded by the minimum (from below) and the maximum (from above) of the sum by row of the matrix entries; dividing by the sample size, we bound 𝑣 cos with minimum and maximum average (by row).

We observed minor differences in the choice of the distance, albeit they have different properties. The cosine distance is suitable to measure morphological (or shape) differences between two P-waves; it is invariant under scaling, i.e. 𝑑 cos (𝑎𝑉 𝑖 , 𝑏𝑉 𝑗 ) = 𝑑 cos (𝑉 𝑖 , 𝑉 𝑗 ) for any 𝑎, 𝑏. However, the cosine distance is not invariant under shift in amplitude and time. This may represent an issue in clinical setting.

The 𝑣 𝑝 -variability measure may highlight localized differences in P-wave morphology. An example is provided by Figure 4. The shaded area is region spanned by ±2 standard deviations around the average P-wave. An interesting interpretation of the 𝑣 𝑝 -variability measure is the following: for 𝑝 = 1, the variability is half of the shaded area; for 𝑝 = ∞, the variability is the maximum deviation in absolute value from the average P-wave.

Clearly, the 𝑣 𝑝 -variability is unaffected by a constant change of the baseline, but it still may be affected by shifts in the time domain.

Suitability of the eikonal model

We assessed the P-wave variability analysis by using the eikonal model for the activation map and the lead field approach for the P-wave simulation. The main advantage of the eikonal approach is its computational


Ques%on: predictability of evolu%on at the molecular level

Jump then Climb: can rearrangements predict the occurrence of muta%onal bursts?

Guillaume Beslon 1 , Vincent Liard 1 , San6ago F. Elena 2,3 1: INRIA-Beagle team (INSA-Lyon), Lyon, France 2: IBMCP (CSIC-UPV), Valencia, Spain; 3: Santa Fe Ins6tute, Santa Fe NM, USA

• Due to the stochas6c nature of muta6ons, evolu6on is generally supposed to be unpredictable at the molecular level.

• But muta6ons are filtered-out by selec6on which may introduce correla6ons in the muta6onal paSerns.

• There are many different kinds of muta6onal events (switches, InDels, rearrangements, HGT…).

• Some of these events may poten6ate the occurrence of others, resul6ng in a non-random fixa6on.

! How to study this process?

• Modeling and simula6on can be used to study how a random spontaneous muta6onal process can turn into a non-random process when looking at fixed muta6ons.

• We need a model in which muta6onal paSerns can account for the variety of molecular events that can alter real genomes.

• The model should include a complex genotype-to-phenotype map.

• Both proper6es are at the core of the Aevol model (www.aevol.fr).

!

Experimental framework:

• We evolved 30 viral wild-types by simula6ng 200,000 genera6ons of evolu6on under a high muta6on pressure (10 -4 mut.bp -1 .gen -1 ).

• Each WT has been cloned 30x and the 900 clones were further evolved for 30,000 genera6ons. • We analyzed the sequence of fixed muta6ons in terms of (1) effect on fitness, genome size, robustness and evolvability (2) wai6ng 6me between two muta6onal events.

Genome size

Favorable Deleterious Neutral

Switch/InDel Rearrangement

Evolu%onary dynamics:

• 215 of the 900 clones significantly improved their fitness. • Fitness gain ocen occurs during short muta6onal bursts with rapid fixa6on of muta6onal events. • These bursts are characterized by a strong increase of evolvability. • More than 50% of the bursts start with a segmental duplica6on. • Compared to spontaneous rates, muta6ons are rare, except during the bursts.

Wai%ng %me between muta%ons:

Delays from the previous fixa6on event (top) and to the next one (boSom) are es6mated per kind of muta6on for the 215 clones that significantly gain fitness (Hodges-Lehmann es6mator). Segmental duplica6ons show a strong skew: they are fixed acer a "muta6onal desert" and are likely to be immediately followed by another muta6on fixa6on event. InDels are also skewed although the skew is less pronounced.

In our experiment evolu%on proceeds by "jump-and-climb" steps:

1. The viruses climb their local fitness peak. This process mainly relies on subs6tu6ons. 2. At the top of the fitness peak, no more favorable subs6tu6ons are available. S6ll many rearrangements remain to be tested. 3. A rearrangement is fixed; viruses jump to a new peak where new favorable subs6tu6ons are available. The climbing process starts again. This sequen6al process enables par6al predic6on at the molecular level: fixa6on of a rearrangement opens the path to new adapta6ons…

The jump-and-climb process is rooted in the combinatorics of muta%onal events In compacted genomes, like viral ones, the combinatorics of point muta6ons is quickly exhausted. Yet, the combinatorics of rearrangements is much larger and cannot be explored in a reasonable 6me. When fixed, they open new paths in the fitness landscape that enable fixa6on of previously impossible point muta6ons… Similar processes have been observed in viruses (e.g. Chikungunya) and bacteria (Blount et al., 2012). Our results open three important ques6ons: (1) is this process restricted to short, compact, genomes or can it be generalized, e.g. to cancer evolu6on? (2) Are there other "jumping" muta6onal events (3) can these events be used to predict disease emergence or evolu6on of drug resistance?

  Here we used Aevol to test the interac%ons between the different kind of muta%ons… Method: In Silico experimental evolu%on with the Aevol model Results: Random spontaneous events don't fix independently Discussion: Can rearrangements be used as predictor of molecular evolu%on? The Aevol model: Aevol is an In Silico Experimental E v o l u 6 o n p l a [ o r m t h a t m o d e l s microorganisms evolu6on with explicit selec6on and replica6on processes (A). Aevol uses a realis6c genome structure (B.1) and a sound genotype-tophenotype map (B). All func6onal levels are modeled as mathema6cal func6ons (B.2-3). Fitness is computed by comparing the phenotype with a predefined target (in red on B.3). Muta6on operators include chromosomal rearrangements (C. 1), switches and Indels (C.2).


Experimental design and complexity measures Experimental design:

• To unravel the origin of molecular complexity, we evolved popula8ons in the simplest possible environment: the Aevol target is a triangle. • We evolved 300 popula8ons of 1024 individuals for 250,000 genera8ons under 3 muta8on rates and monitored the evolu8on of genomic and func8onal complexity.

The complexity ratchet: Stronger than selec;on! Vincent Liard, Jonathan Rouzaud-Cornabas, David P. Parsons, Guillaume Beslon INRIA-Beagle team (INSA-Lyon), Lyon, France Using the Aevol digital gene8cs plaWorm we designed an in silico experiment to study the rela8onship between molecular complexity and phenotypic complexity: We evolved popula8ons of digital organisms in an environment designed to allow survival of the simplest possible organism: one which genome encodes a single gene. By repeatedly evolving popula8ons in this experimental framework, we observed that ≈1/3 of the lineages quickly found this simple genotype and were then stable for the rest of the experiment. At the same 8me, most lineages were not able to find this simple solu8on and showed slow gene acquisi8on along the 250,000 genera8ons of the experiment. Importantly, simple organisms ended up with a very high fitness while complex genotypes ended up with a ≈10x lower fitness. This shows that, even in a simple environment, evolu8on leads to a complexity ratchet: each gene acquisi8on creates the poten8al for the acquisi8on of further genes, ul8mately pushing evolu8on towards complex solu8ons even in a simple environment. Moreover, organisms engaged in this complexifica8on process were never able to outcompete the simple ones, showing that selec8on is not able to invert the complexity ratchet.

Aevol (www.aevol.fr) is an In Silico Experimental Evolu8on (ISEEaka digital gene8cs) plaWorm developed by the Beagle team to study the evolu8on of genome structure. Aevol is based on three principles that makes it perfectly suited to study the evolu8on of complexity:

  Complexity measures:• Qualita8ve measure: "simple" organisms are those encoding only proteins with the same m and w values. • Genomic complexity: quan8ty of informa8on encoded on the genome (total amount of coding sequences). • Func8onal complexity: quan8ty of informa8on encoded on the proteome (number of different parameters).

Discussion

Results

(1) Organisms evolved complex func;onal structures in 66% of the simula;ons Whatever the muta8on rate, ≈1/3 of the simula8ons led to "simple" organisms with few genes and a low func8onal complexity (A). ≈2/3 of the simula8ons led to "complex" organisms despite the simplicity of the target func8on (B).

(2) Complex organisms accumulate more informa;on at the genomic and func;onal levels Genomic complexity is strongly bounded by muta8on rates (A) due to robustness constraints on the genome (Knibbe et al., 2007;Fischer et al., 2014). Muta8on rates also constrain the func8onal complexity (B) but this effect is less stringent at the func8onal level.

(3) Simple organisms are fiKer than complex ones Whatever the complexity measure, we observe a clear trend for simple organisms to be fifer than complex ones ager 250,000 genera8ons. This demonstrates that in our simula8ons complexity is not driven by selec8on.

On the opposite, complex func8onal structures have evolved in spite of selec8on. (4) Despite the advantage of being simple, complex organisms evolve greater complexity on the long term The simple/complex iden88es are determined early on in the simula8on and generally conserved thereager (A).

Complex organisms evolve greater complexity (B); their fitness grows but remains far below simple organisms.

A. Genome and proteome of a simple organism B. Genome and proteome of a complex organism A. Distribu8on of genomic complexity for complex (top) and simple (bofom) organisms. The higher the muta8on rate, the lower the genomic complexity. Genomic complexity is strongly limited by muta8onal robustness.

B. Distribu8on of func8onal complexity for complex (top) and simple (bofom) organisms. The higher the muta8on rate, the lower the func8onal complexity.

A. Fitness at genera8on 250,000 vs genomic complexity. The higher the genomic complexity, the lower the fitness. Simple organisms approach the op8mum fitness ( f opt = 1 ). Mean fitness of complex organisms: f = 0.38.

A. Star8ng from simple organisms at genera8on 0, organisms' iden8ty (simple vs complex) is determined before genera8on 10,000 and generally maintained for the rest of the simula8on.

B. Long-term evolu8on of func8onal complexity in a complex organism. Func8onal complexity and fitness con8nuously grow during the 250,000 genera8ons but the fitness remains far below that of simple organisms. 250,000 vs func8onal complexity. The higher the func8onal complexity, the lower the fitness.

The emergence of complex organisms in a simple environment is a strong argument in favor of a complexity ratchet, i.e. an irreversible mechanism that adds components to a system but that cannot get rid of exis8ng ones, even though this could be more favorable. Indeed, in our experiments this ratchet clicks and goes on clicking despite the selec8ve advantage of being simple. Evolu8on of fitness in complex organisms shows that the ratchet is empowered by nega8ve epistasis. Our results show that complex biological structures can flourish in condi8ons where complexity is not needed and that, reciprocally, the global func8on of complex structures could very well be simple.

Experiment-specific target func8on: this triangular target func8on can be fifed by a single gene/protein.


Alain Pecker
email: alain.pecker@geodynamique.com

An estimate of maximum ground surface motion for non zero surface velocity

Keywords: Soils, Earthquake, Wave propagation, Soil strength Mots-clés : Sols, Séisme, Propagation d'ondes, Résistance du sol

An estimate of maximum ground surface motion for non zero surface velocity

Introduction

The increasing use of Probabilistic Seismic Hazard Assessment (PSHA) for critical facilities leads to the consideration of earthquake scenarios with a probability of occurrence as low as 10 -7 per year. For such low probabilities the computed ground accelerations, based on extrapolations of statistical attenuation relationships, may reach values as high as a few g, which poses tremendous difficulties for the earthquake resistance design [START_REF] Bommer | The challenge of defining upper limits on earthquake ground motions[END_REF]. Obviously such large motions cannot exist because the soil profile, through which the seismic waves travel to reach the ground surface, has a limited resistance capacity and cannot transmit any motion. When failure is reached at any depth within the soil profile, the incident motion is filtered and no motion larger than the motion reached at that stage can be transmitted to the upper strata. The soil acts as a fuse to safely limit the maximum ground surface acceleration. It is therefore necessary to develop an approach that takes into account the soil resistance to derive an upper bound estimate of the ground surface motion.

The approach presented herein is an extension of the previous approach developed in [START_REF] Pecker | An estimate of maximum ground motion[END_REF], for which is has been shown that, not only reasonable estimates of the maximum peak ground acceleration at the ground surface can be obtained, but predictions are not contradicted by observations. It is based on an analytical solution to the wave equation in an inhomogeneous soil profile: in [START_REF] Pecker | An estimate of maximum ground motion[END_REF] the stiffness of the soil column was assumed to be zero at the ground surface and to increase with depth; the condition of zero velocity at the ground surface is relaxed in the present paper. Therefore, almost any, reasonably smooth, experimental variation of the shear wave velocity with depth can be accommodated.

Equations of motion

Let us consider a soil layer of finite thickness overlying a stiff bedrock, which for the purpose of this study will be considered as a rigid boundary. The soil is assumed to be isotropic elastic, with a shear wave velocity increasing with depth according to some power law:

V (z) = V S z + d d + h p/2 (1)
where h is the layer thickness, p a real positive parameter smaller than 2 and V S the shear wave velocity at depth h; d is a strictly positive parameter that can be chosen to fit Eq. ( 1) to the experimental data. In [START_REF] Pecker | An estimate of maximum ground motion[END_REF] this parameter was equal to zero giving a zero shear wave velocity at the ground surface. It is convenient to make a change of variable and to define ζ = (z + d)/(d + h) and H = d + h. The wave equation for a plane vertically incident shear wave can be written:

V 2 S H 2 ∂ ∂ζ ζ p ∂u ∂ζ = ∂ 2 u ∂t 2 + vg (t) ( 2
)
where u is the horizontal displacement, relative to the bedrock, and v g the bedrock displacement motion. The boundary conditions express that the relative displacement at the bedrock interface and the shear stress at the ground surface are equal to 0:

ζ = d H = ζ 0 , τ(ζ 0 ) = 0 ⇒ ∂u ∂ζ ζ =ζ 0 = 0 (3) ζ = 1, u(1, t) = 0
Eq. ( 2) is solved with the mode superposition technique. The modes equation can be written:

d dζ ζ p dX dζ + H 2 V 2 S ω 2 X(ζ ) = 0 ( 4 )
and its solution is given by (p < 2):

X(ζ ) = ζ (1-p)/2 AJ ν λζ (2-p)/2 + BY ν λζ (2-p)/2 (5)
where J ν (•) and Y ν (•) are Bessel's functions of the first and second kind and of order ν

= (p -1)/(2 -p); λ = 2ωH /(2 -p)V S .
Taking into account the recurrence formula, [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF], valid for J ν (•) and Y ν (•):

C ν (z) = -C ν+1 (z) + ν z C ν (z) (6)
The derivative of Eq. ( 5) can be written:

X (ζ ) = - 2 -p 2 λζ (1-2p)/2 AJ ν+1 λζ (2-p)/2 + BY ν+1 λζ (2-p)/2 (7)
Enforcing the boundary conditions, one obtains a linear system of two equations for the unknowns A and B which has a non zero solution provided that:

J ν (λ)Y ν+1 λζ (2-p)/2 0 -Y ν (λ)J ν+1 λζ (2-p)/2 0 = 0 ( 8 )
Eq. ( 8) is the frequency equation with roots ρ i ; the soil column frequencies are given by:

f i = ρ i V S (2 -p) 4πH = ρ i V S (2 -p) 4π(h + d) (9)
The roots ρ i can be computed once and for all as functions of ζ 0 ; an example, which will serve for illustration below, is given in Fig. 1(a) for p = 1.0. Finally for convenience the mode shapes are normalized to unity at the ground surface ζ = ζ 0 ; the final solution for the mode shapes is given, in dimensionless form, by:

X i (ζ ) = π 2 ρ i ζ 0 ζ (1-p)/2 Y (p-1)/(2-p) ρ i ζ (2-p)/2 J 1/(2-p) ρ i ζ (2-p)/2 0 -J (p-1)/(2-p) ρ i ζ (2-p)/2 Y 1/(2-p) ρ i ζ (2-p)/2 0 (10)
Following the methodology of [START_REF] Pecker | An estimate of maximum ground motion[END_REF], the maximum ground surface acceleration due to the contribution of the first N modes is: where S a (ω i , ξ i ) is the spectral pseudo acceleration and α i the modal participation factor:

ümax (z = 0) = N i=1 α i S a (ω i , ξ i ) 2 1/2 (11)
α i = 1 ζ 0 X i (ζ ) dζ 1 ζ 0 X 2 i (ζ ) dζ (12)
As for the roots, the modal participation factors can be computed once and for all as a function of ζ 0 (Fig. 1(b)).

The shear strain is given, for each mode, by:

γ i (z) = ∂u i (z) ∂z = α i S a (ω i , ξ i ) H ω 2 i dX i (ζ ) dζ ( 13
)
which can be expressed in dimensionless form by:

χ i (ζ ) = γ i (z)ω 2 i S a (ω i , ξ i ) H = α i (2 -p)π 4 ρ 2 i ζ 0 ζ (1-2p)/2 J 1/(2-p) ρ i ζ (2-p)/2 Y 1/(2-p) ρ i ζ (2-p)/2 0 -Y 1/(2-p) ρ i ζ (2-p)/2 J 1/(2-p) ρ i ζ (2-p)/2 0 (14)

Soil constitutive model

Although the soil behavior is highly non linear from very small strains, as illustrated by the shear stress-shear strain curve in Fig. 2, a simplified elastic-perfectly plastic constitutive relationship is assumed for the soil layer; the shear stress-shear strain curve is defined by two parameters which may depend on the depth z:

• the shear strength τ max (z); • the yield engineering shear strain γ f (z) where the engineering shear strain is defined as twice the shear strain.

The shear modulus is then given by G = τ max /γ f . As soon as, at any depth within the soil profile, the shear strain reaches γ f , the maximum shear stress that can be transmitted is limited by τ max ; the ground surface acceleration cannot therefore exceed the value reached when γ (z) = γ f . The procedure involves the following steps:

(i) define the input motion at the rock interface by its pseudo acceleration response spectrum S * a . This is typically the result of the PSHA; (ii) compute the eigenfrequencies and mode participation factors from Eqs. ( 8) and ( 12); (iii) plot the shear strain (Eq. ( 14)) versus depth together with the yield strain γ f ; (iv) determine the depth z 1 and the scaling factor µ for which γ f = µγ (z 1 ); (v) define S a = µS * a the maximum possible pseudo acceleration from which the maximum ground surface acceleration ümax (z = 0) is determined (Eq. ( 11)).

Example of application

The procedure outlined above is illustrated with reference to an actual project for which a PSHA and truly non linear site response analyses have led to a ground surface acceleration equal to 4.8 m/s 2 for a 2000 year return period earthquake. The soil profile is composed of 100 m of alluvial deposits underlain by stiffer strata. The undrained shear strength of the top layers has been measured and can be approximated by:

S u (kPa) = 2.85z(m) + 25 (15)
A typical stress strain curve for the alluvia, measured in laboratory tests on a sample retrieved at 15 m depth, is presented in Fig. 2; it is approximated by the simplified elastic perfectly plastic model with a yield strain, determined as shown, equal to γ f = 0.03. The other relevant parameters for the analysis take the following values: p = 1.0, ζ 0 = 0.08 and V S = 71.9 m/s. The 'bedrock' motion, as determined from the PSHA, is represented by the pseudo acceleration response spectrum of Fig. 3, scaled for convenience to 10 m/s 2 . In reality, the result of the PSHA indicates that for the 2000 year return period earthquake the 'rock' acceleration is equal to 5.0 m/s 2 . The response has been computed with 10 modes to achieve a modal mass greater than 97%. The first three circular frequencies are equal to: ω 1 = 0.91 rd/s, ω 2 = 2.26 rd/s, ω 3 = 3.67 rd/s and the corresponding spectral accelerations for 20% damping, which is relevant for a near failure condition, are equal to S * a1 = 0.51 m/s 2 , S * a2 = 3.74 m/s 2 , S * a3 = 8.46 m/s 2 . The yield strain becomes equal to the induced strain for µ 0.59 at a critical depth z 1 = 6 m below the ground surface. The associated spectral accelerations of the first three modes are: S a1 = µS * a1 = 0.30 m/s 2 , S a2 = µS * a2 = 2.21 m/s 2 , S a3 = µS * a3 = 4.99 m/s 2 from which the maximum peak ground acceleration (Eq. ( 11)) is equal to: ümax = 5.6 m/s 2 . The approximation of the stress strain curve by the bilinear elastic plastic model involves some arbitrary choice. If instead of 0.03 the yield strain is taken within the range 0.025 to 0.035, which represents reasonable alternatives, the maximum peak ground acceleration varies between 5.4 and 5.8 m/s 2 . These values compare favorably with the results of the analyses for the 2000 year return period: for those analyses the PSHA indicated a "bedrock" acceleration of 5 m/s 2 and the non linear site response analyses a surface acceleration of 4.8 m/s 2 . The proposed methodology indicates that above a rock ground acceleration of 5.9 m/s 2 the surface motion is limited to 5.4 m/s 2 to 5.8 m/s 2 . These maximum values are in good agreement with the empirical relationships derived from field observations, which indicate maximum surface ground accelerations in the range 4.5 m/s 2 to 5.5 m/s 2 for cohesive soils [START_REF] Seed | Non linear dynamic response of soft and deep cohesive deposits[END_REF].

Conclusions

A simple method has been proposed to estimate the maximum ground surface motion that can be observed at the surface of a soil profile whatever the amplitude of the input rock acceleration. The proposed method is robust because the shear strength is a parameter that is routinely measured, rather reliable, and also because the results are not too sensitive to changes in the yield strain. It has been compared to more rigorous numerical site response analyses based on a sophisticated soil constitutive model from which the ground surface acceleration has been found equal to 4.8 m/s 2 for a potential event with a return period of 2000 years. The method indicates that the maximum ground surface acceleration should not exceed 5.4 to 5.8 m/s 2 for a rock motion slightly larger than that associated with the 2000 year return period event. These values are also in good agreement with the empirical relationships derived from field observations.

Fig. 1 .

 1 Fig. 1. Roots of the frequency equation; (b) mode participation factors.

Fig. 1 .

 1 Fig. 1. (a) Racines de l'équation aux fréquences propres ; (b) facteurs de participation.

Fig. 2 .

 2 Fig. 2. Soil shear stress-shear strain curve.

Fig. 2 .

 2 Fig. 2. Courbe effort-déformation en cisaillement du sol.

Fig. 3 .

 3 Fig. 3. 20% damped bedrock response spectrum.

Fig. 3 .

 3 Fig. 3. Spectre de réponse au rocher à 20 % d'amortissement.


Stéphane Andrieux
email: stephane.andrieux@edf.fr

Huy Duong Bui

Écart à la réciprocité et identification de fissures en thermoélasticité isotrope transitoire

Keywords: Mécanique des solides numérique, Fissures, Écart à la réciprocité Computational solid mechanics, Cracks, Reciprocity gap

On s'intéresse à un solide élastique isotrope renfermant une ou plusieurs fissures planes, et soumis à un chargement d'origine thermique variable dans le temps. On relève la valeur du champ de déplacement et du champ des tractions surfaciques sur le bord du solide. On propose dans cette Note la définition et l'exploitation d'un écart à la réciprocité, formulé uniquement à partir des grandeurs mécaniques connues sur le bord du solide, qui permet de déterminer explicitement le plan renfermant la ou les fissures, les conditions régnant sur celles-ci étant quelconques dès lors qu'elles assurent la continuité du vecteur contrainte et du flux de chaleur normal.

Abridged English version

This Note addresses the problem of location of planar cracks buried inside a homogeneous isotropic elastic solid. The body is submitted to a thermal loading but only mechanical quantities (displacement field and surface traction fields) are measured on the whole boundary. Equipped with these data the problem is to determine the plane where cracks are lying. Eqs. [START_REF] Andrieux | Identification of planar cracks by complete overdetermined data: inversion formula[END_REF] that are satisfied by the displacement fields and temperature fields are an evolution equation for the temperature and a quasi-static equation for the displacement field. Following a previous approach, we define a linear form [START_REF] Ben Abda | Reciprocity principle and crack identification in transient thermal problems[END_REF] acting on auxiliary displacement fields that fulfill the elastic equilibrium equations in the uncracked body [START_REF] Andrieux | Reciprocity principle and crack identification[END_REF]. This linear form is the Reciprocity Gap functional, previously introduced for similar identification problems for the stationary or transient heat (or more general conduction) equation or for the elastostatic equation [START_REF] Andrieux | Identification of planar cracks by complete overdetermined data: inversion formula[END_REF][START_REF] Ben Abda | Reciprocity principle and crack identification in transient thermal problems[END_REF][START_REF] Andrieux | Reciprocity principle and crack identification[END_REF]. It makes use of a surface integral only over the external boundary of the body. The expression [START_REF] Friedman | Determining cracks from boundary measurements[END_REF] gives the interpretation of the reciprocity gap in terms of the jump of the displacement fields across the cracks and a volume integral involving the unknown temperature field inside the solid. Then, following an idea presented in [START_REF] Bui | A path-independent integral for mixed modes of fracture in linear thermoelasticity[END_REF], a constraint (5) on the auxiliary fields is designed in order to cancel the volume integral in [START_REF] Friedman | Determining cracks from boundary measurements[END_REF]. Equipped with this result, the exploitation of the reciprocity gap for some particular polynomial auxiliary fields leads to explicit formulas for the determination of the normal to the plane containing the cracks (10)-( 13) and to the scalar constant appearing in its affine equation ( 14). Then the plane can be entirely determined without need for the resolution of any elastic or conduction equations. These results extend previous formulas given by the authors [START_REF] Andrieux | Identification of planar cracks by complete overdetermined data: inversion formula[END_REF][START_REF] Andrieux | Reciprocity principle and crack identification[END_REF], but the situation is here more interesting because a thermal loading can being applied to a stress free body, the determination of the crack plane can be achieved with only a measurement at different times of the surface displacement field, without any information about the thermal loading.

Introduction

On s'intéresse à un solide tridimensionnel contenant éventuellement une ou plusieurs fissures portées par un même plan. Sous sollicitations thermomécaniques, excluant cependant les forces de masses, un champ de température et un champ de déplacement se développent au sein du solide. On suppose connus, sur tout le bord extérieur, le champ de déplacement, le champ de densité surfacique de force. On propose dans cette Note la définition et l'exploitation d'un écart à la réciprocité, formulé uniquement à partir des grandeurs mécaniques connues sur le bord du solide, qui permet de déterminer explicitement le plan renfermant la ou les fissures. Ces résultats étendent les résultats obtenus en thermique stationnaire [START_REF] Andrieux | Identification of planar cracks by complete overdetermined data: inversion formula[END_REF] ou transitoire [START_REF] Ben Abda | Reciprocity principle and crack identification in transient thermal problems[END_REF] et en élasticité [START_REF] Andrieux | Reciprocity principle and crack identification[END_REF]. Ils peuvent cependant conduire à des méthodes de contrôle plus réalisables puisque l'on pourrait se borner à imposer à un solide, libre de tout chargement mécanique, un chargement thermique, dont il est inutile de connaître la nature (contrairement au cas thermique transitoire traité en [START_REF] Ben Abda | Reciprocity principle and crack identification in transient thermal problems[END_REF]), et de relever le champ de déplacement en surface. L'exploitation de ces données à différents instants permet par ailleurs de lever l'ambiguïté sur la normale au plan de fissure qui subsiste en élasticité lorsqu'une seule mesure est disponible [START_REF] Andrieux | Reciprocity principle and crack identification[END_REF].

En notant u et τ les champs de déplacement et de température, A le tenseur de rigidité (isotrope et homogène), α le coefficient de dilatation linéaire, k la conductivité (également isotrope et homogène) et ρc la chaleur massique, les équations d'équilibre thermoélastique et d'évolution thermique sont les suivantes dans le cadre des petites perturbations et négligeant les effets d'inertie mécaniques (pour alléger l'écriture nous omettons les mesures dans les intégrales) :

             Ω A : ε(u) -ατ I : ε(w) = S ext F .w - Γ A : ε(u) -ατ I .N .JwK ∀w ∈ H 1 (Ω\Γ ) 3 , ∀t Ω k∇τ.∇ϕ + ρc Ṫ ϕ = Ω sϕ + S ext Φϕ - Γ k∇τ.NJϕK ∀ϕ ∈ H 1 (Ω\Γ ), ∀t (1)
Cette écriture n'est valable que si sur les lèvres des fissures règnent des conditions aux limites imposant la continuité du flux normal de chaleur et du vecteur contrainte, conditions qu'il est cependant inutile ici de préciser plus. En particulier des conditions de contact avec frottement peuvent être prises en compte sans difficulté, la seule limitation est que ces conditions ne conduisent pas pour les chargements étudiés à une fermeture totale de la fissure avec blocage du glissement tangentiel puisque dans ce cas les fissures deviennent indiscernables du milieu non fissuré sur le plan mécanique.

F et Φ sont respectivement la densité de force et de flux de chaleur imposé sur le bord du solide, s la source de chaleur. La discontinuité d'un champ à la traversée de fissures est notée à l'aide du symbole J.K :

Jf K = f (x + ) - f (x -), x + ∈ Γ + , x -∈ Γ -, Γ + (respectivement Γ -) étant la face de normale extérieure -N (respectivement N ).
Enfin la température de référence du matériau constituant le solide est prise égale à zéro.

Le champ de déplacement U m et la densité de force F étant mesurés sur le bord extérieur S ext du solide à divers instants (en nombre fini), le problème d'identification des fissures consiste à déterminer tout d'abord le plan renfermant celles-ci puis d'en déterminer la forme dans ce plan, à partir des données surabondantes que constitue le couple (U m , F ). Dans cette version du problème d'identification, on ne suppose connus ni les sources s, ni les flux de chaleurs Φ imposés au solide, ni même la température en surface. De ce fait, la connaissance des conditions initiales en température pour le problème d'évolution posé sur celle-ci est également superflue. On se borne dans cette Note à donner les résultats d'identification obtenus pour le plan de fissure. Les résultats sur la géométrie de la fissure, plus longs à présenter, feront l'objet d'une publication ultérieure.

Si la question d'identification de fissures à l'aide de mesures thermiques (problème thermique seul) ou mécanique (problème élastostatique seul) a fait l'objet de nombreux travaux [START_REF] Andrieux | Identification of planar cracks by complete overdetermined data: inversion formula[END_REF][START_REF] Ben Abda | Reciprocity principle and crack identification in transient thermal problems[END_REF][START_REF] Andrieux | Reciprocity principle and crack identification[END_REF][START_REF] Friedman | Determining cracks from boundary measurements[END_REF][START_REF] Alessandrini | Determining two-dimensional cracks in three-dimensional bodies: uniqueness and stability[END_REF][START_REF] Ben Abda | Identification of 2D cracks by boundary elasticity measurements[END_REF][START_REF] Nakamura | Unique continuation for elliptic systems and crack determination in anisotropic elasticity[END_REF], fort peu en revanche se sont intéressés au problème thermoélastique. De plus, la formulation du problème inverse géométrique d'identification de fissures basé sur les équations de la thermoélasticité transitoire et utilisant uniquement des mesures mécaniques n'a jamais été abordée jusqu'ici à la connaissance des auteurs.

Écart à la réciprocité et condition sur les champs auxiliaires

Le principe de l'écart à la réciprocité consiste à écrire l'expression de réciprocité de Maxwell-Betti, non pas entre deux solutions du problème élastique posé sur un même domaine, mais entre une solution élastique sur le domaine renfermant potentiellement des fissures (le champ est alors le champ réel) et une solution du problème élastique posé sur le domaine non fissuré. Cette expression n'est alors pas identiquement nulle, sa valeur pour chaque champ en équilibre élastique dans le domaine sain fournit des informations sur les discontinuités de déplacement à la traversée des fissures. Il est possible d'exploiter ces informations dans le cas de fissures planes pour déterminer le plan contenant les fissures.

Si v est un champ de déplacement solution du problème d'élasticité dans le domaine sain, c'est-à-dire vérifiant :

Ω A : ε(v) : ε(w) = S ext A : ε(v).n.w ∀w ∈ H 1 (Ω) 3 (2)
on définit la forme linéaire ER t suivante, dite écart à la réciprocité à l'instant t :

ER t (v) = S ext A : ε(v).n.U m (t) -F (t).v (3)
L'écart à la réciprocité ER t se calcule pour tout champ v par une intégrale sur le bord du solide ne faisant intervenir que des quantités connues (F , U m ) ; il possède une propriété intéressante qui permet de faire apparaître explicitement le domaine Γ occupé par les fissures :

ER t (v) = Γ A : ε(v).N . q u(t) y + Ω 3Kατ (t) div(v) pour tout v satisfaisant (2) ( 4
)
où K est le module de compressibilité élastique du matériau qui s'exprime à partir des deux coefficients de Lamé : 3K = 3λ + 2µ. On constate aisément que l'écart est identiquement nul si le solide est à sa température de référence et ne contient aucune fissure. Dans cette expression figure le champ de température réel τ à l'intérieur du solide, champ bien entendu inconnu. Pour lever cette difficulté et ne plus faire apparaître que le premier terme portant sur le domaine occupé par les fissures, on impose une contrainte supplémentaire sur les champs auxiliaires v. Suivant une idée due à Bui [START_REF] Bui | A path-independent integral for mixed modes of fracture in linear thermoelasticity[END_REF] et utilisée dans la construction d'intégrales invariantes en mécanique linéaire de la rupture, on considère les champs v * j suffisamment réguliers qui satisfont pour l'indice j l'équation :

v * j j = 0 dans Ω (sans sommation) (5)
c'est-à-dire les champs pour lesquels la j ème composante est harmonique dans l'ouvert non fissuré. On a alors la propriété suivante (sans sommation) :

ER t v * j ,j = µ Γ ε v * j ,j : N ⊗ q u(t) y s pour tout v * j régulier satisfaisant (2) et (5) (6)
La régularité nécessaire est ici : v * j

,j ∈ H 1 (Ω). Le résultat (6) s'établit aisément grâce aux équations de Navier (forme locale de (2)) que satisfait le champ v * j .

Formules d'identification du plan de fissure

Muni de la définition (3) de l'écart à la réciprocité et de la propriété (6) obtenue pour la famille de champs auxiliaires satisfaisant :

                       Pour j = 1, 2 ou 3 v * j ,j ∈ H 1 (Ω) 3 Ω A : ε(v * j ) : ε(w) = S ext A : ε(v * j ).n.w ∀w ∈ H 1 (Ω) 3 Ω ∇v * j j .∇ϕ = S ext ∇v * j j .nϕ ∀ϕ ∈ H 1 (Ω) (7)
on peut formuler les résultats principaux de cette Note, qui permettent sans aucune résolution de problèmes thermoélastiques, d'identifier complètement le plan renfermant (par hypothèse) l'ensemble des fissures.

Formules d'identification de la normale au plan de fissure(s)

Introduisons le tenseur symétrique P d'ordre deux défini par :

P = 1 2 Γ JuK ⊗ N + N ⊗ Γ JuK def = (U ⊗ N) s (8)
et les champs de déplacements v * i et w * dont les composantes dans un repère cartésien sont données par les expressions suivantes :

v * 1 = 2x 2 1 -x 2 2 -x 2 3 -2x 1 x 2 -2x 1 x 3 , v * 2 = -2x 1 x 2 2x 2 2 -x 2 1 -x 2 3 -2x 2 x 3 , v * 3 = -2x 1 x 3 -2x 3 x 2 2x 2 3 -x 2 2 -x 2 1 , w * = 2x 2 x 3 2x 3 x 1 2x 1 x 2 (9)
alors seules les composantes de la partie déviatorique P de P peuvent être déterminées à partir du calcul de la valeur de la forme linéaire d'écart à la réciprocité sur les champs ainsi définis. Plus précisément :

       Pii = 1 12µ ER t v * i ,i , sans sommation Pij = |ε ij k | 8µ ER t w * ,k , i = j (10)
où ε ij k sont les composantes du tenseur de permutation.

Puisque N est un vecteur unitaire, on déduit de la connaissance du tenseur P les résultats suivants (seuls cas possibles pour les valeurs principales de P) :

• Si P possède une valeur principale double Λ :

Alors, la normale au plan de fissure est donnée par la direction principale associée à la valeur principale simple, l'intégrale du saut de déplacement sur les fissures est portée par N et a pour intensité 3Λ/2. • Si P possède trois valeurs principales simples (Λ 1 , Λ 2 , Λ 3 ) ordonnées par valeurs décroissantes, les directions principales normées étant (V 1 , V 2 , V 3 ), et en notant Λ ij = Λ i -Λ j , la normale N et l'intégrale du saut de déplacement U sur les fissures sont données par :

N = 1 √ Λ 13 Λ 12 V 1 -Λ 23 V 3 , U = -3Λ 2 N + 2 Λ 12 Λ 23 Λ 13 Λ 23 V 1 + Λ 12 V 3 (11)
Du fait de la symétrie du tenseur P, on ne peut discriminer les vecteurs unitaires N et Γ JuK Γ JuK = U U . Comme dans [START_REF] Andrieux | Reciprocity principle and crack identification[END_REF], il faut donc a priori deux instants distincts lors d'un chargement thermique ou deux chargements mécaniques distincts pour identifier la normale au plan des fissures qui sera le vecteur commun aux deux couples de vecteurs unitaires identifiés grâce aux formules précédentes.

Il est également indispensable que le chargement mécanique ou thermique conduise à une solution thermoélastique pour laquelle l'une au moins des trois composantes de la discontinuité moyenne de déplacement ne soit pas nulle. Cette condition peut cependant être vérifiée en examinant la nullité éventuelle du tenseur P obtenu par les Éqs. (10).

Formule d'identification complète du plan de fissure(s)

Connaissant la normale N , on effectue un changement de coordonnées tel que dans le nouveau repère (O, x, y, z) l'axe Oz soit dirigé par le vecteur N et que la composante selon Ox du saut de déplacement moyen sur les fissures ne soit pas nulle (cette dernière condition pouvant être satisfaite via un nouveau changement de repère à partir de l'identification des composantes de la moyenne du saut de déplacement fournie par l'Éq. ( 11)). L'équation du plan affine renfermant les fissures est maintenant : z -C = 0. Le scalaire C est déterminé par la formule suivante :

C = 1 6µU x ER t h * ,x
où h * = 3x(z 2y 2 ) y 3 -3z 2 y 0 (12)

Remerciements

Les auteurs remercient les rapporteurs anonymes pour leurs remarques constructives qui ont permis d'améliorer très sensiblement la présentation des résultats décrits dans cette Note.