Loading [MathJax]/jax/output/SVG/jax.js
Research article Special Issues

Towards virtual surgery planning: the modified Blalock-Taussig Shunt

  • Received: 18 March 2020 Accepted: 03 June 2020 Published: 09 June 2020
  • A modified Blalock-Taussig shunt (MBTS) is an aortopulmonary shunt to establish or augment pulmonary perfusion in congenital cardiac defects with limited pulmonary blood flow. Proper function of this shunt is of utmost importance. In clinical practice, prediction of flow in an MBTS relies on previous experience. In the research field, computational modeling techniques have been developed to simulate flow in an MBTS and predict its performance. These techniques are promising but also time consuming and prone to uncertainties; therefore not yet suitable for clinical practice. Here we present a simplified, patient-based computational model to predict mean circulatory flow characteristics after MBTS insertion. Simulations performed over a range of pulmonary vascular resistances, were compared to data from: i) previous modeling studies; ii) data from the specific patient modeled, and iii) a cohort of patients with MBTS. Model predictions were within one standard deviation from cohort data; and within 1% from results of previous (more complex) computational models. In comparison to previous studies, our model is computationally stable with significantly shorter computational time to perform simulations. We envision that our approach could be used in the future to perform virtual surgeries, quickly testing different surgical scenarios using the patient own geometrical and physiological characteristics, to aid surgeons in decision making.

    Citation: Stephen Haller, Rabin Gerrah, Sandra Rugonyi. Towards virtual surgery planning: the modified Blalock-Taussig Shunt[J]. AIMS Biophysics, 2020, 7(3): 169-188. doi: 10.3934/biophy.2020014

    Related Papers:

    [1] Valentin Duruisseaux, Melvin Leok . Time-adaptive Lagrangian variational integrators for accelerated optimization. Journal of Geometric Mechanics, 2023, 15(1): 224-255. doi: 10.3934/jgm.2023010
    [2] Elias Maciel, Inocencio Ortiz, Christian E. Schaerer . A Herglotz-based integrator for nonholonomic mechanical systems. Journal of Geometric Mechanics, 2023, 15(1): 287-318. doi: 10.3934/jgm.2023012
    [3] Álvaro Rodríguez Abella, Melvin Leok . Discrete Dirac reduction of implicit Lagrangian systems with abelian symmetry groups. Journal of Geometric Mechanics, 2023, 15(1): 319-356. doi: 10.3934/jgm.2023013
    [4] Maulik Bhatt, Amit K. Sanyal, Srikant Sukumar . Asymptotically stable optimal multi-rate rigid body attitude estimation based on lagrange-d'alembert principle. Journal of Geometric Mechanics, 2023, 15(1): 73-97. doi: 10.3934/jgm.2023004
    [5] Jordi Gaset, Arnau Mas . A variational derivation of the field equations of an action-dependent Einstein-Hilbert Lagrangian. Journal of Geometric Mechanics, 2023, 15(1): 357-374. doi: 10.3934/jgm.2023014
    [6] Marcin Zając . The dressing field method in gauge theories - geometric approach. Journal of Geometric Mechanics, 2023, 15(1): 128-146. doi: 10.3934/jgm.2023007
    [7] Jacob R. Goodman . Local minimizers for variational obstacle avoidance on Riemannian manifolds. Journal of Geometric Mechanics, 2023, 15(1): 59-72. doi: 10.3934/jgm.2023003
    [8] Francesco Bonechi, Jian Qiu, Marco Tarlini . Generalised Kähler structure on CP2 and elliptic functions. Journal of Geometric Mechanics, 2023, 15(1): 188-223. doi: 10.3934/jgm.2023009
    [9] William Clark, Anthony Bloch . Existence of invariant volumes in nonholonomic systems subject to nonlinear constraints. Journal of Geometric Mechanics, 2023, 15(1): 256-286. doi: 10.3934/jgm.2023011
    [10] Xavier Rivas, Daniel Torres . Lagrangian–Hamiltonian formalism for cocontact systems. Journal of Geometric Mechanics, 2023, 15(1): 1-26. doi: 10.3934/jgm.2023001
  • A modified Blalock-Taussig shunt (MBTS) is an aortopulmonary shunt to establish or augment pulmonary perfusion in congenital cardiac defects with limited pulmonary blood flow. Proper function of this shunt is of utmost importance. In clinical practice, prediction of flow in an MBTS relies on previous experience. In the research field, computational modeling techniques have been developed to simulate flow in an MBTS and predict its performance. These techniques are promising but also time consuming and prone to uncertainties; therefore not yet suitable for clinical practice. Here we present a simplified, patient-based computational model to predict mean circulatory flow characteristics after MBTS insertion. Simulations performed over a range of pulmonary vascular resistances, were compared to data from: i) previous modeling studies; ii) data from the specific patient modeled, and iii) a cohort of patients with MBTS. Model predictions were within one standard deviation from cohort data; and within 1% from results of previous (more complex) computational models. In comparison to previous studies, our model is computationally stable with significantly shorter computational time to perform simulations. We envision that our approach could be used in the future to perform virtual surgeries, quickly testing different surgical scenarios using the patient own geometrical and physiological characteristics, to aid surgeons in decision making.



    While the forward error of a numerical method compares the exact solution of an ODE with the numerical solution after one time-step h, to obtain qualitative statements about the long-term behaviour of numerical solutions to ODEs, it is helpful to consider a modified ODE whose exact solution closely approximates the numerical flow map at grid points. A modified equation can be obtained as an expansion of the numerical solution as a power series in the step-size h. Though the series does not converge in general, optimal truncation techniques have been established such that the flow of the modified system and the numerical method coincide at grid points up to exponentially small error terms. Computing and analysing the structural properties of modified equations is known as backward error analysis (BEA) (see, for instance, [3,§Ⅸ] or [6,§5]). Next to the analysis of long-term behaviour of numerical schemes, backward error analysis has been used to improve the initialisation of multi-step methods [2] as well as to improve physics informed machine learning techniques [13,12,10].

    If a Hamiltonian ODE is discretised by a symplectic integrator, then any truncation of the modified equation is itself a Hamiltonian system with respect to the original symplectic structure and a modified Hamiltonian. These are also called Shadow Hamiltonians. The existence of a modified Hamiltonian or a modified Lagrangian is a key ingredient to obtain statements about long-term behaviour of symplectic method, such as oscillatory energy errors over exponentially long time intervals. Moreover, just as in the exact system, symplectic symmetries of the modified system yield conserved quantities for the modified dynamics by Noether's theorem. This explains why symplectic integrators behave well on completely integrable systems. A detailed discussion can be found in [3]. Vermeeren observed that backward error analysis for variational integrators can be done entirely on the Lagrangian side [15].

    In contrast to symplectic methods, conjugate symplectic methods preserve a modified symplectic structure rather than the original symplectic structure. Conjugate symplectic methods share the excellent long-term behaviour of symplectic methods. Moreover, Noether's theorem applies such that symmetries of the modified system yield modified conserved quantities of the modified dynamics. When modified structures are explicitly known, explicit expressions of modified conserved quantities can be derived. This motivates the development of techniques to compute modified symplectic structures and Hamiltonians.

    While traditional methods for the computation of modified Hamiltonians use an Ansatz (i.e. an educated guess) of the Hamiltonian as a power series and match terms, working with an Ansatz is challenging when a modified Hamiltonian and a modified symplectic structure need to be computed simultaneously: the components of matrices representing symplectic structures are in this context not constant but depend on the state space variables, fulfil a symmetry condition, and satisfy the Jacobi identity, which is a partial differential equation [3,§Ⅶ.2]. This makes finding a suitable Ansatz difficult.

    A typical strategy [3,8] to obtain a structure preserving numerical method is to approximate the variational principle

    δS=0,S(y)=tNt0L(y(t),˙y(t))dt,y(t0)=y0,y(tN)=yN (1.1)

    which governs the exact Euler–Lagrange equations

    ddtL˙yLy=0

    by a discrete variational principle

    SΔ({yi}i)=0,SΔ({yi}i)=N1i=0hLΔ(yi,yi+1).

    Since y0 and yN are fixed in the variations considered in (1.1), the gradient above is taken with respect to all inner grid points y1,,yN1. We obtain the discrete Euler–Lagrange equations

    D1LΔ(yi,yi+1)+D2LΔ(yi1,yi)=0,i=0,,N1

    which yield approximations yiy(t0+ih) to an exact solution y. The term recursion is called a variational method. Indeed, the class of variational methods is equivalent to the class of symplectic integrators [3,8].

    While backward error analysis for discrete Lagrangians LΔ(yi,yi+1) are established, discrete Lagrangians depending on several grid-points LΔ(yi,yi+1,,yi+s) corresponding to multistep methods require different approaches because they are not symplectic but only preserve a modified symplectic structure. In other words, these methods are conjugate to symplectic methods. However, the modified symplectic structures or conjugacies, respectively, are given by a formal power series that might not be convergent. Although rigorous optimal truncation results are not available, we will demonstrate in numerical examples that truncations can be useful objects in the analysis of the numerical methods.

    In the following, we will introduce blended backward error analysis to systematically compute modified Hamiltonian and modified symplectic structures. We will prove the following theorem which applies, for instance, to series expansions LΔ of consistent discrete Lagrangians LΔ(y(t),y(t+h),,y(t+sh)).

    Theorem 1.1. Consider a power series LΔ(y[]) in a formal variable h. The series depends on the jet y[]=(y,˙y,¨y,) of a variable y such that any truncation only depends on a finite jet of y. Assume further that the truncation to zeroth order constitutes a regular Lagrangian L(y,˙y), i.e. (2L˙yi˙yj)i,j is invertible.

    There exists a 2nd order modified equation given as a formal power series in h such that for any NN a solution of the modified equation truncated to order O(hN) solves the Euler–Lagrange equations to L[N]Δ up to an error of order O(hN+1), where L[N]Δ is the truncation of LΔ to order O(hN).

    If we denote z=(y,˙y), there exists a symplectic structure matrix Jmod and a Hamiltonian Hmod given as formal power series in h such that solutions to Hamilton's equations

    ˙z=J[N]mod(z)1H[N]mod(z)

    fulfil the modified equation. Here, J[N]mod and H[N]mod are truncations to order O(hN) of Jmod and Hmod, respectively such that L[N]Δ is regular, i.e. (2L[N]Δy(M)  iy(M)  j)i,j is invertible, where y(M) is the highest derivative of L[N]Δ.

    The technique will be illustrated for linear multistep methods with matrix valued coefficients (2). These occur, for instance, when in a system of coupled ODEs components of the differential equation are discretised separately with traditional linear multistep method, when multistep methods are stabilised [4,5,14], or when discretisation schemes of PDEs are analysed [9].

    Moreover, we will analyse under which conditions modified Lagrangians exist: if the original equation is the Euler–Lagrange equation to a variational principle of the form (1.1) for a Lagrangian L(y,˙y), it is natural to ask, whether there exists a modified Lagrangian Lmod(y,˙y) such that the modified equation is governed by the Euler–Lagrange equations to Lmod. In contrast to classical variational integrators, for which Lmod(y,˙y) is known to exist and can be computed [15], for conjugate symplectic schemes Lmod may only exists in modified variables Lmod(˜y,˙˜y). We will prove that Lmod exists in the original variables (y,˙y) if Jmod has the form

    Jmod=(0).

    In particular, this shows that for classical consistent, symmetric, linear multistep methods with matrix coefficients with a central force evaluation applied to second order equations Lmod exists in the original variables (y,˙y) if all coefficients are multiples of the identity matrix, i.e. we have a multistep method with scalar coefficients. However, in the general case of matrix coefficients Lmod only exists in modified variables (˜y,˙˜y).

    The article is structured as follows: Section 2 illustrates the ideas of blended backward error analysis on linear multi-step methods. Section 3 shows how to compute the modified data Jmod and Hmod introduced in 1.1. The technique is then applied to linear multi-step methods with matrix valued coefficients in 4 and results are illustrated by numerical experiments. Additionally, for comparison of blended backward error analysis with classical backward error analysis, ?? contains an application of blended backward error analysis to a mechanical ordinary differential equation discretised with the Störmer-Verlet scheme. A formal proof of 1.1 is provided in 5. Finally, 6 discusses the existence of modified Lagrangians as formal power series and future research directions are suggested in 7.

    To illustrate the idea of blended backward error analysis, we compute modified symplectic structures and Hamiltonians of linear multistep methods.

    Consistent, symmetric linear multistep methods with a single force evaluation applied to the second order ODE

    ¨y=U(y(t)) (2.1)

    take the form

    s2j=1Aj(y(tjh)2y(t)+y(t+jh))=h2U(y(t)) (2.2)

    with

    s2j=1j2Aj=I. (2.3)

    Relation (2.3) is coming from the consistency requirement [4]. These are s-step methods, where s is even. Here we allow matrix valued coefficients Aj [4,5,9,14], h is the step-size and I denotes the identity matrix. If the coefficients Aj are scalars, then the schemes constitute classical consistent, symmetric linear multistep methods. A series expansion of (2.2) in h is equivalent to a power series expression of the form

    ¨y=U(y(t))+Ni=1hi˜gi(y(t),,y(ai))+O(hN+1) (2.4)

    for some h-independent expressions ˜gi in the ai-jet of y at t. Substituting 3rd and higher derivatives on the right hand side of (2.4) with derivatives of (2.4) itself iteratively yields an equation of the form

    ¨y=U(y(t))+Ni=1higi(y(t),˙y(t))+O(hN+1) (2.5)

    which is called the modified equation of method (2.2) applied to (2.1). We refer to [3] for optimal truncation techniques and a discussion of spurious solutions not covered by the considered modified system for the case of linear multistep methods with scalar coefficients. In the following, we will focus on the question which structural properties the modified equation (2.5) shares with the original ODE (2.1).

    Variational Structure

    The ODE (2.1) has first order variational structure as it is the Euler–Lagrange equation

    ddtL˙yLy=0

    to the variational principle

    δS=0forS(y)=L(y(t),˙y(t))dt

    with

    L(y,˙y)=12˙y2+U(y).

    Moreover, there exists a variational principle for (2.2):

    Lemma 2.1. Let the matrices AjRn×n be symmetric. For T>0 let T either be the circle T=R/TZ or the real line T=R. For y defined on T the variational principle

    δSΔ=0forSΔ(y)=TLΔ(y(t),y(t+h),,y(t+s2h))dt (2.6)

    with

    LΔ(y(t),y(t+h),,y(t+s2h))=12h2(s2j=1Aj(y(t+jh)y(t)),y(t+jh)y(t))+U(y(t)).

    implies the functional equation (2.2). Moreover, if T=[a,b] is an interval, then (2.6) implies (2.2) on the interval T=[a+s2h,bs2h]. Here we assume that the function space for y and the potential U are such that Uy and Uy constitute square integrable functions.

    Proof. Let Δτy denote the forward difference, i.e. (Δτy)(t)=y(t+τ)y(t). Then y,ΔτzL2(T,Rn)=Δτy,zL2(T,Rn) holds with Δτ=Δτ on T if T{R,R/TZ} or on the interval [a+τ,bτ] if T=[a,b]. The expression ΔτΔτy corresponds to the central difference

    (ΔτΔτy)(t)=y(t+τ)+2y(t)y(tτ).

    Let δ denote the variational derivative in the direction of a variation δy, i.e. δSΔ(y)=limϵ01ϵ(SΔ(y+ϵδy)SΔ(y)). Let δyCc(T,Rn), if T{R,R/TZ} and let δyCc(T,Rn), if T=[a,b]. We have

    δSΔ(y)=12h2s2j=1δAjΔjhy,ΔjhyL2(T,Rn)+δTU(y(t))dt=1h2s2j=1AjΔjhy,ΔjhδyL2(T,Rn)+U(y),δyL2(T,Rn)=1h2s2j=1AjΔjhΔjhy+U(y),δyL2(T,Rn).

    Now (2.2) follows from the fundamental lemma of the calculus of variations on T or T, respectively.

    To analyse structure preserving properties of the method (2.2), it might seem natural to seek a modified Lagrangian Lmod(y,˙y) given as a formal power series in the step-size h such that the modified variational principle

    δSmod=0forSmod(y)=Lmod(y(t),˙y(t))dt

    covers smooth solutions of (2.6) up to any order in the step-size h. However, we show that although a 1st order Lagrangian Lmod covering the modified equations always exists as a power series, it only exist in modified variables (˜y,˙˜y) in the most general case. Even for simple methods, the existence of an expression in closed form for the change of coordinates from (y,˙y) to (˜y,˙˜y) can not be expected. This makes it difficult to compute Lmod using an ansatz.

    Hamiltonian structure

    Another approach is to work on the Hamiltonian side. The ODE (2.1) has the form of a Hamiltonian system

    ˙z(t)=J1H(z(t)),H(z)=12˙y2U(y),z=(y˙y). (2.7)

    Here

    J=(0II0) (2.8)

    is the standard symplectic structure.

    In this paper we use a blended approach of the variational and Hamiltonian viewpoint to systematically compute a modified Hamiltonian system

    ˙z(t)=J1mod(z(t))Hmod(z(t)) (2.9)

    consisting of a modified Hamiltonian Hmod and a modified symplectic structure Jmod given as formal power series in h such that for a truncation to arbitrary order N the ODE (2.9) covers the modified equation (2.5) up to higher order terms. Here Jmod is a skew symmetric matrix which satisfies a Jacobi identity. The following theorem can be considered as an instance of 1.1 and will be proved in 5.

    Theorem 2.2. Let NN denote the considered order of the series expansion of the matrix multistep method (2.2), where the coefficients are symmetric matrices. There exists a modified symplectic structure J[N]mod, which is O(h) close to J and a modified Hamiltonian H[N]mod, which is O(h) close to H, such that

    ˙z(t)=(J[N]mod(z))1H[N]mod(z(t)) (2.10)

    with coordinate z=(y,˙y) is equivalent to the modified equation (2.5) up to terms of order O(hN+1).

    Here O(h)-closeness of J[N]mod and J means that the zeroth coefficient of the polynomial J[N]mod in the formal variable h is given by J. O(h)-closeness of H[N]mod and H has an analogous meaning.

    We observe conditions under which a modified first order Lagrangian Lmod(y,˙y) exists in the original variable y by analysing the modified symplectic structure Jmod.

    Theorem 2.3. If the matrices Aj in (2.2) are scalar multiples of the identity matrix, then there exists a modified Lagrangian L[N]mod depending on (y,˙y) that is O(h)-close to L(y,˙y) such that the modified variational principle

    δS[N]mod=0forS[N]mod(y)=L[N]mod(y(t),˙y(t))dt

    yields the modified equation (2.5) up to terms of order O(hN+1).

    Again, O(h)-closeness is to be interpreted in a formal sense, analogously to its meaning in 2.2.

    Theorem 2.3 applies to traditional multistep methods with scalar coefficients. However, we will see that for general linear multistep methods with matrix-valued coefficients the existence of a first order modified Lagrangian in the original variable y cannot be expected. A proof of 2.3 is postponed to 6.

    In this section we introduce a method to compute the modified data J[N]mod and H[N]mod of 2.2 such that (2.10) governs (2.5). We will then verify the validity of the construction method and prove 1.1 and 2.2 in 5.

    Let L[N]Δ denote the series expansion of

    LΔ(y(t),y(t+h),,y(t+s2h))

    to order N in the step-size h. The expression L[N]Δ depends on the M-jet of y at t for some MN. Notice that the order M variational principle

    δS[N]Δ(y)=0withS[N]Δ(y)=L[N]Δ(y(t),˙y(t),,y(M)(t))dt (3.1)

    recovers (2.4) up to higher order terms. We first compute a high-dimensional Hamiltonian system defined on the 2M1-jet space of y corresponding to the order M variational principle (3.1). The Hamiltonian principle is then reduced to a Hamiltonian system defined on the 1-jet space of y. It has the form (2.9) and covers (2.5) up to higher order terms.

    To construct the high-dimensional Hamiltonian system, we use Ostrogradsky's Hamiltonian description of high-order Lagrangian systems [16]. For this we define variables

    q=(y,˙y,,y(M1))

    and for i=1,,M

    pji=Mik=0(1)kdkdtkL[N]Δ(yj)(k+i).

    Here the index j enumerates the components of y and ddt denotes the total derivative operator on the jet-space of y, which acts like

    ddtρ(y,˙y,,y(M))=Mi=0y(i)ρ,y(i+1).

    on a scalar valued function ρ defined on the M-jet space of y. The high dimensional Hamiltonian system consists of the Hamiltonian

    H[N]=Mi=1pi,˙qiL[N]Δ, (3.2)

    where all expressions are expressed in y[2M1]=(y,˙y,,y(2M1)), and the symplectic structure matrix J[N]mod(y[2M1]). The skew-symmetric matrix J[N]mod(y[2M1]) is the representing matrix of the differential 2-form

    Ω[N]=Mi=1nj=1dpjidqij, (3.3)

    where pji and qij are interpreted as functions in the variable y[2M1] of the 2M1-jet space, i.e. J[N]mod is the anti-symmetrised tensor product* of the gradients y[2M1]  pij and y[2M1]  qij summed over all indices.

    *This corresponds to the command TensorWedge in Wolfram Mathematica.

    To compute the modified Hamiltonian system on the 1-jet space with variable y[1]=(y,˙y), the variables y(2),,y(2M1) in the expression (3.2) for the Hamiltonian H[N](y[2M1]) are repeatedly replaced by (2.5) until higher derivatives only occur in O(hN+1) terms. This yields H[N]mod(y,˙y). Similarly, we can consider pji and qij as functions of (y,˙y) truncating O(hN+1) terms. The matrix J[N]mod is then given as the representing matrix of the 2-form Ω[N] pulled to the 1-jet space of the variable y, i.e. interpreting y[1]=(y,˙y) as the only independent variables. Equivalently, J[N]mod is the anti-symmetrised tensor product of the gradients y[1]pij and y[1]qij summed over all indices. (As J[N]mod is constructed from a closed differential 2-form, it is automatically skew-symmetric and satisfies the Jacobi identity.) This completes the construction of the modified data.

    The system (2.9) recovers (2.5) up to higher order terms as we will prove after a computational example.

    As introduced in 2, consider the multistep method

    A2y(t2h)+A1y(th)2(A1+A2)y(t)+A1y(t+h)+A2y(t+2h)=h2U(y(t)) (4.1)

    with matrix coefficients in dimension n=2. By the consistency requirement, A2=1/4(IA1). We obtain

    J[4]mod=J+h2(J211J221J2210)+h4(J411J421J421J422)

    with

    J211=(0b1b10)

    where

    b1=14(˙y2(a12(U(2,1)U(0,3))+(a22a11)U(1,2))+14(˙y1(a12(U(2,1)U(0,3))+(a22a11)U(2,1))).

    Here and in the following U(k,l)=k+lU(y)ky1ly2.

    As the expressions of higher order terms become quite complicated, we refer the reader to the Mathematica Notebooks of our accompanying source code [11]. However, as this will be relevant in the discussion later, we are reporting J422 for the special case that A1 is a diagonal matrix: we have

    J422=(0b2b20)

    with

    b2=18(a211a222+2(a22a11))(U(2,1)˙y1+U(1,2)˙y2)ifA1=(a1100a22).

    The modified Hamiltonian H[4]mod for the general case is given as

    H[4]mod(y,˙y)=H(y,˙y)+h2H2(y,˙y)+h4H4(y,˙y)

    with

    H2(y,˙y)=124((˙y1)2(2(43a11)U(2,0)6a12U(1,1))2˙y2˙y1(3a12U(0,2)+(3a11+3a228)U(1,1)+3a12U(2,0))6a22U(0,2)(˙y2)26a12U(1,1)(˙y2)2+(3a224)(U(0,1))2+3a11(U(1,0))2+6a12U(0,1)U(1,0)+8U(0,2)(˙y2)24(U(1,0))2).

    For further terms, we refer the reader to the Mathematica Notebooks of our accompanying source code. Hamilton's equations

    ˙z=(J[4]mod(z))1H[4]mod(z)

    for the modified Hamiltonian system are equivalent to the modified equation (2.5) truncating terms of order O(h6).

    Figure 1 shows a numerical experiment with a rotational invariant potential. The start values for the multistep formula were obtained using the fourth order modified equation. Trajectories computed with the multistep scheme look very regular. The quantities H[0]mod=H, H[2]mod, and H[4]mod evaluated along a trajectory show oscillatory energy error behaviour. Experiments with different values for the step-size h confirm the preservation of H[k]mod, k=0,2,4 up to truncation error. Initialising the multistep scheme with the fourth order modified equation, the effects of spurious solutions was minimised. However, as h is decreased, spurious solutions cause wriggles in the energy error H[k]modH(zinit), k=0,2,4 which eventually prevent further energy error decay.

    Figure 1.  Numerical experiment with the multistep scheme (4.1) with A2=diag(0.85,1.25), A1=I4A2, U(y)=exp(12(y21+y22)) and time-step h=0.1. The multistep formula was initialised at times 0, h, 2h, 3h by integrating the order 4 modified equation using very fine Euler-steps starting from the initial value (yinit,˙yinit)=((1,1),(0.1,0.2)). The figure at the top shows a trajectory up to time t=500, which looks like an orbit in a completely integrable system. The figures below show an evaluation along the trajectory of HH((yinit,˙yinit)) (blue) as well as H[2]modH[2]mod((yinit,˙yinit)) (orange) and H[4]modH[4]mod((yinit,˙yinit)) (green) up to time t=5000 and t=10, respectively. We see oscillatory energy error behaviour.

    If A2=αI with α14 and A1=I4A2, then (4.1) corresponds to a classical stable multistep scheme: the generating polynomial ρ to (4.1) is given as

    A multistep scheme for second order equations is stable if all roots of its generating polynomial lie in the closed unit disk and those on the circle are at most double zeros [3,XV.1.2].

    ρ(ξ)=α+(14α)ξ2(13α)ξ2+(14α)ξ3+αξ4.

    The polynomial ρ has a double root at 1 as well as the roots

    ξ3=112α+14α2α,ξ4=112α14α2α.

    Since α14, the roots ξ3 and ξ4 are complex conjugate to each other and lie on the unit circle such that the scheme is stable.

    Moreover, LΔ and LΔ are rotationally invariant because A1 and A2 commute with rotation matrices. An application of Noether's theorem to L[N]Δ yields the following modified angular momentum:

    I[N]=Mm=1m1k=0(1)kykL[N]Δ,dRy(m1k),dR=(0110)

    The integer M is the order of the highest derivative of y in L[N]Δ. For the truncation order N=4 we have M=5. Using (2.5) repeatedly, derivatives of y of order greater than two are replaced by terms in y, ˙y and O(h6)-terms. After truncation at order 4 in h we obtain the modified angular momentum I[4]mod(y,˙y). Figure 2 shows the conservation error of I=I[0]mod, I[2]mod, I[4]mod along a trajectory. We see oscillatory error behaviour. The oscillations shrink as the step-size h decreases until spurious solutions prevent further decrease.

    Figure 2.  The plots in the second row show the error II(yinit,˙yinit) (blue) in the angular momentum and I[2]modI[2]mod(yinit,˙yinit) (orange) and I[4]modI[4]mod(yinit,˙yinit) (green) along a trajectory. The same data as in Figure 1 was used apart from A2=diag(0.3,0.3) and h=0.15. The plot at the top shows a trajectory initialised as in Figure 1.

    For source code of our experiments refer to [11].

    Let us prepare the proof of 1.1. We will use the language of Differential Geometry.

    To motivate the following proposition and to fix notation recall the following fact: if (Z,ω,H) is a Hamiltonian system with Hamiltonian vector field XH, ZΨZ an embedding such that Ψ(Z) is a symplectic submanifold that is invariant under motions, i.e. XH(Ψ(z))TΨ(z)Ψ(Z) for all zZ, then the pull-back system (Z,ω,H) with ω=Ψω, H=HΨ is a Hamiltonian system with Hamiltonian vector field XH. Here Ψω denotes the pull-back of ω along Ψ. The Hamiltonian vector fields XH and XH relate by

    ΨXH=XHΨ,

    i.e. for a motion γ with ˙γ=XHγ the curve γ=Ψγ is a motion of (Z,ω,H), i.e. ˙γ=XHγ. Here Ψ denotes the push-forward map, i.e. (ΨXH)(z)=dΨ|z(XH(z)) for zZ.

    In the following, we will adapt this statement to a setting, where the definition of ω, H, and Ψ contain a formal variable h and where Ψ(Z) is left invariant only up to higher order terms in h.

    Definition 5.1 (Formal Hamiltonian system). Let NN be a truncation index, Z be a smooth manifold, h a formal variable, and H=Nk=0hkHk a formal polynomial whose coefficients are smooth maps ZR. Further, let ω=Nk=0hkωk be a formal polynomial whose coefficients are 2-forms on Z. The collection (Z,ω,H) is called a formal Hamiltonian system with truncation index N if the following conditions are satisfied.

    ● The formal symplectic form ω is closed, i.e. dω=Nk=0hkdωk is a formal polynomial whose coefficients are 3-forms that are zero.

    ● The 2-form ω0 is non-degenerate.

    The relation dH=ω(X,) defines a formal power series X=k=0hkXk, where Xk are vector fields on Z. The truncation XH:=Nk=0hkXk is called (formal) Hamiltonian vector field. The formal differential equation ˙γ=XHγ is called Hamilton's equation.

    Proposition 5.2. Let NN be a truncation index, Z and Z be manifolds, and let (Z,ω,H) be a formal Hamiltonian system with truncation index N. Consider a formal polynomial Ψ=Nk=0Ψkhk such that Ψk:ZZ are smooth and Ψ0:ZZ is an embedding. To zZ define formal tangent spaces

    TΨ(z)Ψ(Z):={Nk=0hkdΨk|z(v)|vTzZ}Nk=0hkTzkZ,

    where z:=Nk=0hkzk:=Ψ(z).

    Assume that the Hamiltonian vector field is tangential to Ψ(Z), i.e. trk(XH(Ψ(z)))TΨ(z)Ψ(Z) for all zZ. Here trk denotes the truncation of O(hN+1)-terms.

    Assume that Ψ(Z) is a symplectic submanifold of Z, i.e. for all zZ

    TΨ(z)Ψ(Z)TΨ(z)Ψ(Z)ω=0.

    Here Ψ(Z)ω denotes the symplectic complement

    TzΨ(Z)ω={vNk=0hkTzkZ|trk(ω(v,w))=0wTzΨ(Z)},

    where z:=Nk=0hkzk:=Ψ(z).

    Then the pull-back system (Z,ω,H) with ω=Ψω, H=HΨ constitutes a formal Hamiltonian system with truncation index N and

    trk(ΨXH)=XHΨ.

    Proof. Let ω=Nk=0hkωk Since Ψ0 is an embedding, it is an immersion and we conclude that ω0 is non-degenerate as ω0 is non-degenerate. Since pull-back and the differential d commute, ω is closed.

    In the following calculations we identify two formal polynomials P1 and P2 with coefficients of the same type (real numbers, n-forms, smooth functions, ) if and only if P1P2O(hN+1). For all zZ we have

    ωΨ(z)(dΨ|z(XH(z)),dΨ|z())=(Ψω)z(XH(z),)=ωz(XH(z),)=dH|z=d(HΨ)|z=dH|Ψ(z)dΨ|z=ωΨ(z)(XH(Ψ(z)),dΨ|z())
    ωΨ(z)(dΨ|z(XH(z))XH(Ψ(z)),dΨ|z())=0dΨ|z(XH(z))XH(Ψ(z))TΨ(z)Ψ(Z)ω.

    As the inclusion dΨ|z(XH(z))XH(Ψ(z))TΨ(z)Ψ(Z) holds as well and Ψ(Z) is a symplectic submanifold, it follows that

    dΨ|z(XH(z))=XH(Ψ(z)).

    We can now proceed to the proof of 1.1.

    Proof. Step 1. Construction and validity of the modified equation. The truncated power series L[N]Δ is a formal polynomial of the form

    L[N]Δ(y[M])=L(y,˙y)+Nk=1hkLkΔ(y[M]),

    where h is the formal variable. Since the Lagrangian L is regular, the Euler–Lagrange equations to L[N]Δ

    Mj=0(1)jddt(y(j)L[N]Δ(y[M]))=0

    are equivalent to an ordinary differential equation of the form

    ¨y=g0(y,˙y)+Nk=1hk˜gk(y[2M]) (5.1)

    when truncating terms of order O(h[N+1]). The ordinary differential equation is of order 2M because L[N]Δ(y[M]) is regular as well. Iterative replacements of derivatives of order j2 by derivatives of (5.1) yield the modified equation

    ¨y=g0(y,˙y)+Nk=1hkgk(y,˙y)+O(h[N+1]). (5.2)

    Solutions to

    ¨y=g0(y,˙y)+Nk=1hkgk(y[M]) (5.3)

    fulfil (5.1) up to O(h[N+1]) terms by construction. This proves the first part of 1.1.

    Step 2. Existence of Hamiltonian structure on a higher jet-space. Let Y denote the domain of y (smooth manifold), Z=Jet1(Y) the 1-jet space, and Z=Jet2M1(Y).

    The iterative substitution procedure gives rise to a formal polynomial Ψ of maps ZZ defined by y[1]y[2M1], where y(j) for j2 is expressed as a formal polynomial of functions depending on y[1]=(y,˙y). The expression for y(j) is obtained by deriving (5.1) j2 times, iteratively replacing derivatives of order greater than 2 by derivatives of (5.1) followed by a truncation of O(h[N+1]) terms.

    In the remainder of this part of the proof we suppress the fact that we operate on formal polynomials as the following steps can also be done when h is substituted with a sufficiently small real number h>0.

    In the following, we will show that there exists a Hamiltonian structure (Z,ω,H) on Z whose motions correspond to the motions induced by the order M-Lagrangian L[N]Δ(y[M]). To compute the Hamiltonian structure (Z,ω,H), we first construct a Hamiltonian system on TJetM1(Y) which we will then pull back to Z.

    As the Lagrangian L[N]Δ(y[M]) is regular, by Ostrogradsky's principle for high order Lagrangians [16], there exists a transformation χ:ZTJetM1(Y) between the jet variable y[2M1]Z and variables (q,p)TJetM1(Y) with

    q=(y,˙y,,y(M1))

    and

    pi=Mik=0(1)kdkdtky(k+i)L[N]Δ(y[M]),i=1,,M

    such that with

    Ω=Mi=1nj=1dpjidqij,andH[N](q,p)=Mk=1pk,qkL[N]Δ

    the motions of (5.1) are exactly mapped to the motions of the Hamiltonian vector field XH[N] on TJetM1(Y) with

    dH[N]=Ω(XH[N],) (5.4)

    by the transformation χ1. Let ω=χΩ and H=H[N]χ. Now (Z,ω,H) is a Hamiltonian system on Z whose motions are exactly the motions induced by the Lagrangian structure. This completes step 2 of the proof.

    Remark 5.3. In the setting of formal polynomials, Ω is a two form whose coefficients are formal polynomials or, alternatively, Ω is a formal polynomial whose coefficients are 2-forms. In this case, (5.4) defines a formal series XH[N] in h from which we truncate O(hN+1) terms. This is also done when defining XH through dH=ω(XH,). The differential equation defined by XH is then equivalent to the differential equation induced by the Lagrangian structure up to O(hN+1) terms.

    Step 3. Pull-back of (Z,ω,H) to Z. As the Lagrangian L is regular, the 2-form ω0=nj=1dyjdpj, where p is obtained by the Legendre transform for L, is a symplectic form. The pull-back form ω=Ψω is a formal polynomial in h, whose coefficients are 2-forms. The zeroth coefficient ω0 coincides with ω0. ω is, therefore, non-degenerate and Ψ(Z) is a symplectic submanifold of Z in the sense of 5.2. Moreover, by construction of H and Ψ the condition trk(XH(Ψ(z)))TΨ(z)Ψ(Z) for all zZ is fulfilled. Therefore, 5.2 completes the theorem's proof.

    Proof of 2.2. The construction method of the modified data J[N]mod and H[N]mod coincides with the construction method verified in the proof of 1.1.

    Proposition 6.1. All Hamiltonian systems

    ˙z(t)=ˉJ1(z(t))ˉH(z(t)),z=(y˙y)

    with a fixed symplectic structure represented by ˉJ can be formulated as variational problems of the form

    δˉS=0forˉS(y)=ˉL(y(t),˙y(t))dt (6.1)

    if and only if ˉJ is of the form

    ˉJ=(0n), (6.2)

    where 0n denotes an n×n-dimensional zero matrix with n the dimension of y, (i.e. if the distribution spanned by the vector fields ˙y1,,˙yn is Lagrangian).

    Proof. Denote the domain of the variable y by Y, the 1-jet space over Y by Jet1(Y), and the symplectic form represented by the matrix ˉJ by ˉω. Consider a Hamiltonian ˉH:Jet1(Y)R. If the distribution D spanned by ˙y1,,˙yn is Lagrangian w.r.t. ˉω, then there exists a primitive ˉλ of ˉω with kernel D [7,Cor.15.7]. The 1-form ˉλ is of the form

    ˉλ=nj=1(y,˙y)dyj.

    By Hamilton's principle, a curve γ:[tstart,tend]Y is a Hamiltonian motion if the action functional

    ˉS(γ)=γ(ˉλˉHdt)

    is stationary w.r.t. variations of γ through smooth curves fixing the endpoints. Thanks to the special structure of ˉλ (absence of d˙yj-terms), the pullback of ˉλˉHdt along a curve γ has the form

    ˉL(y(t),˙y(t))dt

    with y(t) describing γ(t) in the coordinate y. Therefore, in coordinates, the variational principle has the form (6.1).

    On the other hand, a variational principle of the form (6.1) with regular Lagrangian ˉL (i.e. invertible (2ˉL˙yi˙yj)i,j) can be formulated as a Hamiltonian system with Hamiltonian

    ˉH(y,˙y)=˙y,p(y,˙y)L(y,˙y),and withp(y,˙y)=ˉL˙y(y,˙y).

    The symplectic structure is given as

    ˉω=dˉλ,withˉλ=nj=1pj(y,˙y)dyj.

    As the distribution D spanned by ˙y1,,˙yn is in the kernel of a primitive of ˉω, the distribution is Lagrangian.

    Remark 6.2. The strength of 6.1 lies in the assertion that ˉL is a first-order Lagrangian in the original variable y, i.e. it depends on (y,˙y) only. If ˉJ is not of the required form, then, by Darboux's theorem, we can perform a change of variables on Jet1(Y) such that ˉω is the standard symplectic form dˉpidˉqi and ˉL is a 1st-order Lagrangian in q, i.e. depends on (q,˙q) but not on higher derivatives in q. However, as q and p each depend on (y,˙y), the variables have lost their dynamical meaning. This is because the required change of variables on Jet1(Y) is not fibred, i.e. the jet-space structure is not preserved. Expressed in the original variable y, the Lagrangian ˉL then depends on (y,˙y,¨y), i.e. describes a higher order variational structure.

    Remark 6.3. The computational example presented in 4 provides an example for which the modified symplectic structure is not of the form that is required in 6.1, unless A1 is of the form A1=αI, i.e. the method coincides with a multistep method with scalar coefficients. Indeed, if the considered method is a classical multistep method, i.e. all coefficients are scalar, then Lmod(y,˙y) exists by 2.3.

    We now proceed to the proof of 2.3. We exploit that linear multistep methods can be interpreted as 1-step methods on the original phase space [3,§XV.2]. Here and below we refer to the theory of linear multistep methods for 2nd order ODEs.

    Proof of 2.3. As proved in [1,§5], for the underlying 1-step method ϕ of a symmetric linear multistep method there exists a local diffeomorphism ψ such that ˜ϕ=ψϕψ1 is symplectic with respect to the original symplectic structure ω=jdpjdqj. The conjugacy ψ is given as a P-series applied to the original Hamiltonian vector field X0 given by the right hand side of

    ˙q=p˙p=U(q).

    The P-series ψ is a formal power series that is in general not convergent. Conjugacy, pull-back and push-forward operations are to be interpreted in a formal sense. The map ϕ is symplectic with respect to the modified symplectic structure ωmod=ψω and is the time-h-flow of a vector field X, for which the flow equations correspond to a first order formulation of the modified equation (2.5). The ω-symplectic map ˜ϕ is the time-h-flow of the ψ-related vector field ˜X=ψXψ1. By standard results on backward error analysis for symplectic integrators [3,§Ⅸ], ˜X is a Hamiltonian vector field w.r.t. the standard symplectic structure ω for a Hamiltonian H up to any order in the step-size h. Pulling back the Hamiltonian structure (ω,H) using ψ, we obtain a modified Hamiltonian system (ωmod,Hmod)=(ψω,Hψ) such that Hamilton's equations are equivalent to the modified equation (2.5).

    Since ψ is a P-series in X0, the distribution D spanned by the vertical vector fields p1,,pn is Lagrangian w.r.t. ωmod=ψω. By 6.1, for any order N in the step-size h there exists a first order Lagrangian Lmod(y,˙y) in the original variable such that the variational principle

    δ(Lmod(y(t),˙y(t))dt)=0

    recovers the modified equation up to higher order terms in h.

    Remark 6.4. The proof of 2.3 also shows that Lmod in 2.3 has the structure of an S-series applied to a P-series (see [1]). The modified Lagrangian Lmod can, thus, be computed with an ansatz as well. The modified data Hmod and Jmod can then be computed from Lmod by a Legendre transformation.

    Motivated by optimal truncation results for modified equations, it would be interesting to analyse the convergence properties of modified symplectic structures, modified Hamiltonians, and modified Lagrangians. Moreover, in view of 6.4, a systematic description of the combinatorial structure of the modified quantities appears feasible.

    We would like to thank Peter Hydon and Chris Budd for discussions on multisymplecticity and variational principles at the Isaac Newton Institute of the University of Cambridge in 2019. This research was supported by the Marsden Fund of the Royal Society Te Apˉarangi and the School of Fundamental Sciences, Massey University, Manawatˉu, New Zealand. The authors would like to thank the Isaac Newton Institute for Mathematical Sciences for support and hospitality during the programme Geometry, compatibility and structure preservation in computational differential equations (2019) when work on this paper was undertaken. RM would like to thank the Simons Foundation for their support. This work was supported by: EPSRC grant number EP/R014604/1. The authors acknowledge support from the European Union Horizon 2020 research and innovation programmes under the Marie Skodowska-Curie grant agreement No. 691070.


    Acknowledgments



    The authors would like to thank Dr. Chivukula for advice during model implementation. This publication was made possible with support from the Oregon Clinical and Translational Research Institute (OCTRI), grant number UL1TR000128 from the National Center for Advancing Translational Sciences (NCATS), a component of the National Institutes of Health (NIH), and NIH Roadmap for Medical Research; and grant NIH R01 HL094570. The content is solely the responsibility of the authors and does not necessarily represent the official views of grant giving bodies.

    Conflict of interest



    The authors declare no conflicts of interest.

    [1] Petrucci O, O'Brien SM, Jacobs ML, et al. (2011) Risk factors for mortality and morbidity after the neonatal Blalock-Taussig shunt procedure. Ann Thorac Surg 92: 642-652. doi: 10.1016/j.athoracsur.2011.02.030
    [2] Dirks V, Prêtre R, Knirsch W, et al. (2013) Modified Blalock Taussig shunt: a not-so-simple palliative procedure. Eur J Cardio-Thorac Sur 44: 1096-1102. doi: 10.1093/ejcts/ezt172
    [3] Li JS, Yow E, Berezny KY, et al. (2007) Clinical outcomes of palliative surgery including a systemic-to-pulmonary artery shunt in infants with cyanotic congenital heart disease: does aspirin make a difference? Circulation 116: 293-297. doi: 10.1161/CIRCULATIONAHA.106.652172
    [4] Pennati G, Corsini C, Hsia TY, et al. (2013) Computational fluid dynamics models and congenital heart diseases. Front Pediatr 1: 4. doi: 10.3389/fped.2013.00004
    [5] Marsden AL, Feinstein JA (2015) Computational modeling and engineering in pediatric and congenital heart disease. Curr Opin Pediatr 27: 587. doi: 10.1097/MOP.0000000000000269
    [6] Esmaily Moghadam M, Migliavacca F, Vignon-Clementel IE, et al. (2012) Optimization of shunt placement for the Norwood surgery using multi-domain modeling. J Biomech Eng 134: 051002. doi: 10.1115/1.4006814
    [7] Esmaily-Moghadam M, Murtuza B, Hsia TY, et al. (2015) Simulations reveal adverse hemodynamics in patients with multiple systemic to pulmonary shunts. J Biomech Eng 137: 031001. doi: 10.1115/1.4029429
    [8] Lagana K, Balossino R, Migliavacca F, et al. (2005) Multiscale modeling of the cardiovascular system: application to the study of pulmonary and coronary perfusions in the univentricular circulation. J Biomech 38: 1129-1141. doi: 10.1016/j.jbiomech.2004.05.027
    [9] Hsia TY, Cosentino D, Corsini C, et al. (2011) Use of mathematical modeling to compare and predict hemodynamic effects between hybrid and surgical Norwood palliations for hypoplastic left heart syndrome. Circulation 124: S204-S210. doi: 10.1161/CIRCULATIONAHA.110.010769
    [10] Arthurs CJ, Agarwal P, John AV, et al. (2017) Reproducing patient-specific hemodynamics in the Blalock–Taussig circulation using a flexible multi-domain simulation framework: applications for optimal shunt design. Front Pediatr 5: 78. doi: 10.3389/fped.2017.00078
    [11] Kowalski WJ, Teslovich NC, Dur O, et al. (2012) Computational hemodynamic optimization predicts dominant aortic arch selection is driven by embryonic outflow tract orientation in the chick embryo. Biomech Model Mechan 11: 1057-1073. doi: 10.1007/s10237-012-0373-z
    [12] Migliavacca F, Pennati G, Dubini G, et al. (2001) Modeling of the Norwood circulation: effects of shunt size, vascular resistances, and heart rate. Am J Physiol-Heart C Physiol 280: H2076-H2086. doi: 10.1152/ajpheart.2001.280.5.H2076
    [13] Piskin S, Altin HF, Yildiz O, et al. (2017) Hemodynamics of patient-specific aorta-pulmonary shunt configurations. J Biomech 50: 166-171. doi: 10.1016/j.jbiomech.2016.11.014
    [14] Zahorec M, Hrubsova Z, Skrak P, et al. (2011) A comparison of Blalock-Taussig shunts with and without closure of the ductus arteriosus in neonates with pulmonary atresia. Ann Thorac Surg 92: 653-658. doi: 10.1016/j.athoracsur.2011.04.008
    [15] Schiavazzi DE, Arbia G, Baker C, et al. (2016) Uncertainty quantification in virtual surgery hemodynamics predictions for single ventricle palliation. Int J Numer Meth Bio Eng 32: e02737. doi: 10.1002/cnm.2737
    [16] Marsden AL, Feinstein JA, Taylor CA (2008) A computational framework for derivative-free optimization of cardiovascular geometries. Comput Method Appl Mech Eng 197: 1890-1905. doi: 10.1016/j.cma.2007.12.009
    [17] Winberg P, Lundell BP, Gustafsson LE (1994) Effect of inhaled nitric oxide on raised pulmonary vascular resistance in children with congenital heart disease. Heart 71: 282-286. doi: 10.1136/hrt.71.3.282
    [18] Nelson DP, Schwartz SM, Chang AC (2004) Neonatal physiology of the functionally univentricular heart. Cardiol Young 14: 52-60. doi: 10.1017/S1047951104006304
    [19] Turitto VT, Hall CL (1998) Mechanical factors affecting hemostasis and thrombosis. Thromb Res 92: S25-S31. doi: 10.1016/S0049-3848(98)00157-1
    [20] Holme PA, Ørvim U, Hamers MJAG, et al. (1997) Shear-induced platelet activation and platelet microparticle formation at blood flow conditions as in arteries with a severe stenosis. Arterioscl, Throm, Vas Biol 17: 646-653. doi: 10.1161/01.ATV.17.4.646
    [21] Chowdhury UK, Venugopal P, Kothari SS, et al. (2006) Criterions for selection of patients for, and results of, a new technique for construction of the modified Blalock-Taussig shunt. Cardiol Young 16: 463-473. doi: 10.1017/S1047951106000631
    [22] Chowdhury UK, George N, Sankhyan LK, et al. (2019) A novel surgical technique of construction of the modified Blalock-Taussig shunt (UKC's modification): A video presentation. J Heart Cardiovas Med 2: 025-026.
  • This article has been cited by:

    1. Sina Ober-Blöbaum, Christian Offen, Variational learning of Euler–Lagrange dynamics from data, 2023, 421, 03770427, 114780, 10.1016/j.cam.2022.114780
  • Reader Comments
  • © 2020 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(5587) PDF downloads(305) Cited by(0)

Figures and Tables

Figures(8)  /  Tables(5)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog