Scalable and Robust Dual-Primal Newton–Krylov Deluxe Solvers for Cardiac Electrophysiology with Biophysical Ionic Models

The focus of this work is to provide an extensive numerical study of the parallel efficiency and robustness of a staggered dual-primal Newton–Krylov deluxe solver for implicit time discretizations of the Bidomain model. This model describes the propagation of the electrical impulse in the cardiac tissue, by means of a system of parabolic reaction-diffusion partial differential equations. This system is coupled to a system of ordinary differential equations, modeling the ionic currents dynamics. A staggered approach is employed for the solution of a fully implicit time discretization of the problem, where the two systems are solved successively. The arising nonlinear algebraic system is solved with a Newton–Krylov approach, preconditioned by a dual-primal Domain Decomposition algorithm in order to improve convergence. The theoretical analysis and numerical validation of this strategy has been carried out in Huynh et al. (SIAM J. Sci. Comput. 44, B224–B249, 2022) considering only simple ionic models. This paper extends this study to include more complex biophysical ionic models, as well as the presence of ischemic regions, described mathematically by heterogeneous diffusion coefficients with possible discontinuities between subregions. The results of several numerical experiments show robustness and scalability of the proposed parallel solver.


Introduction
Modern computing platforms allow the study and simulation of complex biophysical phenomena with very high accuracy, thanks to the availability of large computing and memory resources on last generation machines, as well as to advanced mathematical and numerical models. As a consequence, the development of accurate and efficient large-scale numerical solvers have become a very important research topic in the parallel scientific computing community. In the cardiac research field, the opportunities given by computational studies of the heart functions, have increased the interactions among clinicians, mathematicians, physicists and engineers, working together to better understand the behavior of the different parts of the heart and to use these studies to predict and simulate several cardiac pathologies [25,30].
The mathematical modeling of the heart involves systems of partial differential equations (PDEs) and ordinary differential equations (ODEs), which are combined together to model the three main cardiac functions: the electrophysiology [4], the induced cardiac mechanical activity [27], and the blood circulation inside the heart [8].
In this paper, we focus on the first function and design efficient numerical solvers for cardiac electrophysiological models which are also scalable and robust. Among the recent works devoted to the study of parallel cardiac solvers as well as the coupling of different electrical and mechanical cardiac models, we refer to [2,6,23,32] and to the references therein. The main objective of this work is to investigate the efficiency and robustness of a decoupled dual-primal Newton-Krylov solver for implicit time discretizations of the Bidomain model, already introduced in [15] but only coupled with simplified ionic models. This paper extends this study to include both more complex and biophysically detailed ionic models, as well as the presence of ischemic regions.
The Bidomain model describes the propagation of the electric signal in the cardiac tissue by means of a system of parabolic PDEs, coupled with a system of ODEs representing the ionic currents dynamics. By using a fully implicit time scheme for the discretization of the temporal variable, we obtain at each time step a nonlinear algebraic system, which is solved by the Newton method. We choose to solve the Jacobian linear system arising at each Newton step with a Krylov iterative solver, together with Balancing Domain Decomposition with Constraints (BDDC) [9,10] and Dual-Primal Finite Element Tearing and Interconnecting (FETI-DP) [12] preconditioners, in order to accelerate convergence. These algorithms belongs to the class of dual-primal Domain Decomposition (DD) methods, where the degrees of freedom (dofs) are classified into those internal to each subdomain and into those on the interface, which are further divided into dual and primal dofs [29].
In previous work, the authors have analyzed and presented a theoretical convergence rate estimate for the preconditioned solver, together with several preliminary parallel tests using a simple phenomenological ionic model. In particular, the proposed solution strategy relies on a staggered approach, where the two PDE and ODE systems are solved successively, in contrast to a monolithic or coupled approach (e.g. [16,22]).
We extend here the numerical results by studying the robustness of the proposed solver both in case biophysical ionic models are employed, such as the Luo-Rudy phase one [19] and the Ten Tusscher-Panfilov [28] human ionic models, and in case ischemic regions are considered. The inclusion of ischemic regions in the computational domain is modeled mathematically by introducing jumps in the diffusion coefficients; in turn, the discontinuity of the diffusion coefficients on the boundaries of the ischemic region impact on the conditioning of the linear systems, thus requiring robust preconditioned iterative solvers. This paper is structured as follows: in Section 2 we review the micro and macroscopic models of cardiac electrophysiology, introducing the application we consider for our dualprimal solver, presented in Section 3, together with the adopted discretization and solution schemes. Section 4 provides extensive parallel numerical tests that show scalability and robustness of the proposed dual-primal Newton-Krylov solver, posing the basis for several possible future extensions, discussed in the conclusive Section 5.

The Cardiac Electrical Model
In the following Section, we briefly review our cardiac reaction -diffusion model, by introducing the assumptions needed for its formulation. Moreover, we introduce the ionic models employed both in the theoretical analysis and in the numerical experiments.

The Bidomain Model
The mathematical description of the electrical activity in the cardiac tissue, known as myocardium, is provided by the Bidomain model.
The myocardium can be represented as the composition of two ohmic conducting media, named intra-( i ) and extracellular ( e ) domains, separated by the active cellular membrane ( ); the latter acts as insulator between the two domains, as otherwise there would be no potential difference across the domain. These anisotropic continuous media are assumed to coexist at every point of the tissue and to be connected by a distributed continuous cellular membrane [5]. Additionally, the cardiac muscle fibers rotate counterclockwise and their arrangement is modeled as laminar sheets running radially from the outer (epicardium) to the inner surface (endocardium) of the heart.
This setting influences the mathematical definition of the conductivity tensors needed for the formulation of the Bidomain equations. In particular, at each point x of the cardiac domain we can define an orthonormal triplet of vectors a l (x), a t (x) and a n (x), respectively parallel to the local fiber direction, tangent and orthogonal to the laminar sheets, and transversal to the fiber axis ( [18]). By denoting with σ i,e l,t,n the conductivity coefficients in the intra-and extracellular domains along the corresponding directions, we define the conductivity tensors D i and D e of the two media as e n a n (x)a T n (x). Structural inhomogeneities in the intra-or extracellular spaces due to the presence of gap junctions, blood vessels and collagen are generally included in the conductivity tensors D i and D e as inhomogeneous functions of space.
Thanks to the hypothesis on the cardiac tissue, the electric potential is defined in each point of the two domains as a quantity averaged over a small volume: consequently, every point of the cardiac tissue is assumed to belong to both the intracellular and the extracellular spaces, thus being assigned both an intra-and an extracellular potential. From now on, we will denote by the cardiac tissue volume represented by the superposition of these two spaces.
The parabolic-parabolic formulation of the Bidomain model can be stated as follows: find the intracellular and extracellular potentials u i,e : × (0, T ) → R, the transmembrane potential v = u i − u e : × (0, T ) → R, the gating variables w : × (0, T ) → R M , the ionic concentration variables c : given I i,e app : × (0, T ) → R intra-and extracellular applied currents and initial values v 0 : → R, w 0 : → R M and c 0 : in .
Here, C m is the membrane capacitance for unit area of the membrane surface and χ is the membrane surface to volume ratio. The Neumann zero boundary conditions in the system (1, last row) represent mathematically the assumption that the heart is electrically insulated. The nonlinear reaction term I ion and the ODEs system for the gating variables w (which represent the opening and closing process of the ionic channels) and the ionic concentrations c are given by the chosen ionic membrane model. Results on existence, uniqueness and regularity of the solution of system (1) have been extensively studied, see for example [5,31].

Ionic Models
In this work, we consider three different ionic models. First, a phenomenological ionic model, derived from a modification of the renowned FitzHugh-Nagumo model [13,14], named the Rogers-McCulloch (RMC) ionic model [24]. This model overcomes the hyperpolarization of the cell during the repolarization phase by adding a nonlinear dependence between the transmembrane potential and the gating variable and by neglecting the ionic concentrations variable. For this model, I ion (v, w) and R(v, w) are given by where G, v th , v p , η 1 and η 2 are given coefficients. We then consider two more biophysically detailed ionic models, namely the Luo-Rudy phase one (LR1) and the Ten Tusscher-Panfilov (TP06) ionic models, which provide a more detailed description of the ionic currents in cardiac cells. For the several equations of these more complex ionic models we refer to the original papers [19,28].
The theoretical analysis presented in [15] considered only the class of phenomenological ionic model (neglecting the ionic concentration variables). This work extends our previous study to include the LR1 and TP06 ionic models, as well as the presence of ischemic regions.

Dual-Primal Newton-Krylov Solvers
In this section, we briefly present our discretization choices for the model, by providing a finite element discretization in space, a fully implicit time discretization in time and a decoupling (or staggered) solution strategy, in the same fashion as in [20,26]. Then, we introduce our dual-primal solver for the arising nonlinear algebraic system and we mention the theoretical convergence result, which can be found in its extended form in [15].

Space and Time Discretizations
We discretize in space the cardiac domain with Q 1 finite elements, leading to the semidiscrete system with A i , A e , M the stiffness and mass matrices arising from the finite element discretization.
Regarding the discretization in time, instead of using implicit-explicit (IMEX) schemes [6,32], where the diffusion term is treated implicitly, while the remaining terms explicitly, or more generally operator splitting strategies [3], we consider here a fully implicit time discretization in a decoupling (or staggered) approach.
In this procedure, at each time step, the ODEs system representing the ionic model is solved first; then, the nonlinear algebraic Bidomain system is solved and updated. In a very schematic way, this strategy can be summarized as: for each time step n, 1. given the intra-and extracellular potentials at the previous time step, define v := u n i −u n e and compute the gating and ionic concentrations variables 2. solve and update the Bidomain nonlinear system. Given u n i,e at the previous time step and given w n+1 and c n+1 , compute u n+1 = (u n+1 i , u n+1 e ) by solving the system We observe that the Jacobian linear system associated to the nonlinear problem in step 2 is symmetric. This strategy is usually adopted in contrast to a monolithic approach, where the PDEs and ODEs systems are solved altogether, and the computational workload is higher due to the presence of the gating and ionic concentrations variables in the nonlinear algebraic system. Nevertheless, this strategy has been extensively studied and several scalable parallel preconditioners have been designed and analyzed (e.g. [16,21,22]).

Dual-primal Preconditioners
In this approach, a linear system has to be solved at each Newton step: to this end, we employ an iterative method, preconditioned by a dual-primal substructuring algorithm.
Let us consider a partition of the computational domain into N nonoverlapping subdomains of diameter H j , such that = ∪ N j =1 j , and we define the subdomain interface as = ∪ N j =1 ∂ j \∂ . Let W (j ) be the associated local finite element spaces. We introduce the product spaces Using this notation, we can decompose W into a primal subspace W which has continuous elements only and a dual subspace W which contains finite element functions which are not continuous. In this work, we will denote with subscripts I , and the interior, dual and primal variables, respectively.
In the substructuring framework, the starting problem Kw = f is reduced to the interface by eliminating the degrees of freedom (dofs) interior to each subdomain, obtaining the Schur complement system where S = K − K I K −1 I I K I and g = f − K I K −1 I I f I are obtained by reordering the dofs in interior and interface.
The reordering of the degrees of freedom allows to define algorithms where each subproblem is solved independently from the others except for the primal constraints, where the variables are assumed to be continuous across the subdomains. For a detailed presentation, we refer to [29].
On these premises, it is possible to built two of the most used dual-primal iterative substructuring algorithms, namely the Balancing Domain Decomposition with Constraints (BDDC) and dual-primal Finite Element Tearing and Interconnecting (FETI-DP) preconditioners.
In order to ensure a correct continuity of the solution across the subdomains, an appropriate interface averaging is needed. In our work, we focus on both standard ρ and deluxe scalings of the dual variables.

Restriction and scaling operators
Before introducing the preconditioners, it is helpful to understand how the scaling procedure works. We define the restriction operators The ρ-scaling can be defined for the Bidomain model at each node x ∈ (j ) as where N x 1 is the set of indices of all subdomains with x in the closure of the subdomain. Conversely, the deluxe scaling (introduced in [7, 10]) computes the averagew = E D w for each face F or edge E equivalence class.
Suppose that F is shared by subdomains j and k . Let S (j ) F and S (k) F be the principal minors obtained from S (j ) and S (k) by extracting all rows and columns related to the degrees of freedom of the face F. Denote with u j,F = R F u j the restriction of u j to the face F through the restriction operator R F . Then, the deluxe average across F can be defined as It is possible to extend this definition when considering the deluxe average across an edge E. Suppose for simplicity that E is shared by only three subdomains with indices j 1 , j 2 and j 3 ; the extension to more than three subdomains is equivalent. Let u j,E = R E u j be the restriction of u j to the edge E through the restriction operator R E and define S (j 123 ) E ; the deluxe average across an edge E is defined as The averageū is constructed with the contributions from the relevant equivalence classes involving the substructure j . These contributions will belong to W , after being extended by zero to \F or \E. The sum of all contributions R T * ū * are then added from the different equivalence classes to obtainū where E D is a projection and is its complementary projection. We define the scaling matrix for each subdomain j as being k 1 , . . . , k j ∈ Ξ * j , a set containing the indices of the subdomains that share the face F or the edge E and where the diagonal blocks are given by D (j ) in case the deluxe scaling is considered. Conversely, if the ρ-scaling is taken into account, the j th diagonal scaling matrix D (j ) contains the pseudo-inverses (3) along the diagonal.
Lastly, we define the scaled local restriction operators BDDC preconditioners BDDC are two-level preconditioners introduced in [9] for the Schur complement system where S = R T S R and f = R T f are obtained with the operator R which is the sum of local operators R (j ) that return the local interface component.
They can be considered as evolution of balancing Neumann-Neumann algorithms, where local and coarse problems are treated additively. In these algorithms, the choice of primal constraints across the interface is important, since it influences the structure and size of the coarse problem and hence the overall convergence of the method.
It is possible to define BDDC preconditioners using the scaled restriction operators R T D, as where the action of the inverse of S can be evaluated as The first term is the sum of local solvers on each subdomain j , while the second is a coarse solver for the primal variables, where are the primal problem and the matrix which maps the primal degrees of freedom to the interface variables respectively. Then, BDDC algorithm for the solution of the Schur complement problem (4) can be defined as: find u ∈ W such that Once the interface u is computed, we can retrieve the internal solution u I by solving local Dirichlet problems. FETI-DP preconditioners FETI-DP preconditioners were first proposed in [12] as an alternative to one-level and two-level FETI and are based on the transposition from the Schur complement problem (2) to a minimization problem on W with continuity constraints on the dual degrees of freedom, by means of additional variables (named Lagrange multipliers). By introducing a jump matrix B, needed to ensure a correct transmission of the solution between subdomains, the resulting saddle point system is further reduced to a problem only in the Lagrange multipliers unknowns Once the Lagrange multipliers are computed and the interface variables w retrieved, it is possible to compute the internal variables w I by solving local problems on each subdomain.

Convergence analysis
We have proved in [15] a theoretical convergence rate estimate for these two dual -primal preconditioned operators, since BDDC and FETI-DP have been proven to share the same essential spectrum in [17], if the same coarse space is chosen.

Parallel Numerical Results
We extend here the parallel numerical experiments presented in [15] by testing the efficiency, scalability and robustness of the proposed dual-primal Newton-Krylov solver in the following settings: -when biophysical ionic models, such as the LR1 and TP06 models, are considered in a physiological situation; -in the presence of an ischemic region, modeled mathematically by introducing jumps in the diffusion coefficients and other ionic parameters.
We consider two simplified geometries, a thin slab and a portion of half truncated ellipsoid, modeling an idealized left ventricle tissue geometry. The parametric equations of the prolate ellipsoid are given by the system with a 1,2 , b 1,2 and c 1,2 given coefficients defining the main axes of the ellipsoid. Regarding the cardiac fibers, we consider a intramural rotation linearly with the depth of 120 • proceeding counterclockwise from the surface corresponding to the outer layer of the tissue (epicardium, r = 1) to the surface corresponding to the inner layer (endocardium, r = 0), see e.g. Fig. 1.
As already stated, we refer to the original papers [19,24,28] for the equations and parameters of Luo-Rudy phase one (LR1), Rogers-McCulloch (RMC) and Ten Tusscher-Panfilov (TP06) ionic models, respectively. Details about the implementation of the ischemic region  will be given in the related paragraph. Additionally, we refer to [5] for the parameters related to the Bidomain model. We apply an external stimulus I app = 100 mA/cm 3 for 1 ms on a small region of the epicardium. Instead, when the slab geometry is considered, the stimulus is applied in one corner of the domain, over a spheric volume of radius 0.1 cm.
We consider a time interval of 2 ms, with fixed time step τ = 0.05 ms, for a total of 40 time steps. This assumption is not restrictive, since with larger time steps it is possible to lose accuracy in the representation of the wavefront propagating in the tissue, while smaller time steps would only increment the computational workload to marginally increase the accuracy of the solution.
The parallel C codes have been developed using the Portable Extensible Toolkit for Scientific Computation (PETSc) library [1] and the numerical tests have been carried out on Indaco cluster from University of Milan (https://www.indaco.unimi.it/), a Linux Infiniband cluster with 16 nodes, each carrying 2 processors Intel Xeon E5-2683 V4 2.1 GHz with 16 cores each, for a total amount of 512 cores. We assign to each processor one local problem, thus resulting in a correspondence between subdomains and processors. We employ the SNES package from PETSc library, which provide a ready-to-use environment for the solution of nonlinear algebraic systems. We compare the performance of BDDC and FETI-DP preconditioners with the one of algebraic multigrid, from the Hypre implementation (boomerAMG, or bAMG, [11]) and from the built-in PETSc implementation (GAMG), both with default parameters. For the optimality tests, we compare both the standard ρ-scaling and the deluxe scaling. In order to test the efficiency of our solver on parallel architectures, we also compute the parallel speedup S p = T 1 T p , the ratio between the runtime T 1 on 1 processor and the average runtime T p on p processors.

Normal Physiological Tests
With start with considering a normal physiological situation for the LR1 and TP06 ionic models and study the weak scalability and optimality of our dual-primal solvers.  The number of nonlinear iterations do not increase with the number of subdomains, and present lower values for TP06 in both geometries. Additionally, this parameter seems to be affected by the complexity of the geometry, since it is higher for the truncated ellipsoid.
The performance of BDDC and FETI-DP in terms of average CPU time show robustness of the preconditioned solver, since this quantity remains bounded while increasing the number of subdomains, except for the case of BDDC for LR1 with 16 processors for which we do not have any clear explanation; this trend cannot be found for bAMG, which presents higher and increasing values.
Regarding the average number of linear iterations, for the slab geometry and for both LR1 and TP06 ionic models, BDDC and FETI-DP present lower and bounded values, while for bAMG these values increase with the number of processors. On the ellipsoidal geometry instead, we have some fluctuations for BDDC and FETI-DP, although the average number of linear iterations remains bounded; the multigrid presents higher values with respect to its corresponding cases on the slab. Tables 3 and 4. As in the preliminary tests with the simple RMC ionic model reported

Optimality tests We report the results of optimality tests for increasing ratio H/h in
in [15], these results are independent of the scaling employed (the ρ-scaling and deluxe scaling, see [10,29] for more details). The average number of nonlinear iterations settles around 2-4 for the LR1 ionic model and around 2-3 for the TP06 model, with higher values for the ellipsoidal domain. The number of nonlinear iterations seems to be independent of the coarse space employed (V, V+E, V+E+F). As confirmed by Fig. 2, the average number of linear iterations per time step deteriorates while increasing the local problem size when the coarse space consists of only subdomain vertex constraints (V). Instead, by enriching the primal space by adding edge (V+E) and face (V+E+F) constraints, this quantity remains bounded and with lower values.
The average CPU times are almost identical if we consider the richest primal spaces (V+E and V+E+F).

Transmural Ischemic Tests
In order to test the robustness of our dual-primal Bidomain decoupled solver, we add jumps in the diffusion coefficients, modelling pathological conditions such as myocardial ischemia. In particular, we consider a small and regular transmural ischemic region located at the center of both geometries (see Fig. 1). We consider only the RMC and TP06 ionic models, investigating the behavior of the dual-primal solver both in case of a simple phenomenological ionic model and in case of a human ventricular ionic model.
In the ischemic region, we decrease the diffusion coefficients σ i l and σ i t along and across fiber as shown in Table 5. Furthermore, we reduce the ionic current by 30% for the RMC ionic model; in case of the TP06 ionic model, we increase the potassium extracellular concentration K o from 5.4 mV to 8 mV, and we decrease the sodium conductance G Na by 30%, simulating a region with moderate ischemia. Others settings can be found in [5,Chapter 9] and the references therein.
We represent in Figs. 3 and 4 the time evolution of the transmembrane potential v and the extracellular potential u e from the epicardial surfaces respectively, with a transmural ischemic region in the middle of the slab geometry.
In these experiments, we employ the PETSc implementation of algebraic multigrid (GAMG), with default parameters.

Weak scalability. Transmural ischemic region
The results of the weak scalability tests can be found in Tables 6 and 7, where we compare the performance of algebraic multigrid and dual-primal preconditioners, with both ionic models and geometries. The local mesh size is  As already observed in the previous comparison in normal physiological cases, the average number of nonlinear iterations for the TP06 model is higher than the RMC model. Additionally, we again observe higher nonlinear iteration counts for the more complex ellipsoidal domain.
Regarding the slab geometry, BDDC and FETI-DP scale well in terms of average number of linear iterations, since this quantity remains bounded while increasing the number of processors, while GAMG's iterations increase. For the ellipsoidal geometry, we obtain more linear iteration fluctuations for all preconditioners, but the dual-primal preconditioners present lower values than GAMG.
BDDC and FETI-DP are slower in terms of average CPU time than GAMG, probably due to the higher need of interprocessors communication. In contrast to these higher values, GAMG presents an increasing trend for increasing processor counts, while the dual-primal preconditioners show a slower CPU time growth.

Strong scalability. Transmural ischemic region
We then consider strong scalability tests, where we fix the global mesh to 192 · 192 · 32 elements (for a total of 2,458,434 dofs) for the slab geometry and we increase the number of subdomains from 32 to 256. In contrast, we fix the global mesh to 128 · 128 · 64 elements for the portion of truncated ellipsoid (thus resulting in 2,163,3360 dofs). The average number of nonlinear iterations increases with the complexity of the ionic current: as reported in Tables 8 and 9, this parameter increases from 1-2 per time step using the RMC model to 2-3 per time step with TP06 model.
The robustness of the solver is confirmed by the average number of linear iterations, which is comparable for all the preconditioners and for both ionic models; moreover, this indicates that our dual-primal solver retains its good convergence properties even for more complex ionic models. Again, as a consequence, the CPU times of TP06 model increase (with respect to RMC CPU times) due to the increase of nonlinear iterations, but the associated parallel speedups of the models are comparable. Since we are working with a low number of processors (thus meaning larger local problems), BDDC and FETI-DP outperform the ideal speedup (both with respect to 32 and 64 processors), since the factorization of the matrices takes most of the computational time.
Optimality tests. Transmural ischemic region Lastly, we report in Tables 10 and 11 the results of optimality tests for increasing ratio H/h. We fix the number of processors (and subdomains) to 64 and we increase the local size H/h from 4 to 24, thus reducing the finite element size h.
Since the dual-primal preconditioners have been proven to be spectrally equivalent, we focus only on the performance of BDDC. For both ionic models, we again consider both scalings (ρ-scaling on top, deluxe scaling at bottom of each table) and we test the solver for   increasing primal spaces including only vertex constraints (V), vertex and edge constraints (V+E) or vertex, edge and face constraints (V+E+F). Similar results hold for both geometries, independently of the ionic model employed or the scaling chosen. Despite higher average CPU timings for the deluxe scaling, all the other parameters are similar between the two scalings (the average number of nonlinear iterations are the same, for each ionic model). If we consider only vertices V in the coarse space, we observe that the number of linear iterations increases with the local size, as shown in Fig. 5. On the other hand, by adding edge and (V+E) and face (V+E+F) constraints in the primal space, we obtain a quasi-optimality condition, where the linear iterations remains bounded, except for the slab geometry with TP06 ionic model, where the ρ-scaling with the full primal space (V+E+F) behaves as the coarsest primal space (V).

Conclusions
We have reviewed BDDC and FETI-DP preconditioners for fully implicit time discretizations of the Bidomain system, solved through a staggered strategy. We have presented extensive parallel numerical experiments testing the efficiency, scalability and robustness of the solver, both in case of complex ionic models and in the presence of regions with moderate ischemia. The results confirm the validity of the proposed solvers, enlarging the class of available methods for the numerical solution of complex reaction -diffusion models. Future studies should investigate alternatives to the Newton method, by exploring the potentiality of quasi-Newton and inexact-Newton algorithms.
Funding Open access funding provided by Università degli Studi di Pavia within the CRUI-CARE Agreement.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.