Bifurcations in valveless pumping techniques from a coupled fluid-structure-electrophysiology model in heart development

We explore an embryonic heart model that couples electrophysiology and muscle-force generation to flow induced using a $2D$ fluid-structure interaction framework based on the immersed boundary method. The propagation of action potentials are coupled to muscular contraction and hence the overall pumping dynamics. In comparison to previous models, the electro-dynamical model does not use prescribed motion to initiate the pumping motion, but rather the pumping dynamics are fully coupled to an underlying electrophysiology model, governed by the FitzHugh-Nagumo equations. Perturbing the diffusion parameter in the FitzHugh-Nagumo model leads to a bifurcation in dynamics of action potential propagation. This bifurcation is able to capture a spectrum of different pumping regimes, with dynamic suction pumping and peristaltic-like pumping at the extremes. We find that more bulk flow is produced within the realm of peristaltic-like pumping.


I. INTRODUCTION
Various kinds of hearts are found throughout the animal kingdom [1], [2], [3]. In particular many invertebrates have valveless, tubular hearts from their infancy throughout adulthood [4]. These tubular hearts are similar to vertebrate heart morphologies during their first stage of vertebrate heart morphogenesis, e.g., the linear heart tube stage. We begin our discussion of heart tube morphologies by considering the evolution of hearts in the animal kingdom. Figure 1 shows the evolution of hearts from tunicates to humans. Tunicates have an open circulatory system from infancy through adulthood, in which blood is pushed through out the organism by a valveless-tubular  [5] illustrating the evolution of hearts from the valveless heart tubes in the open circulatory systems of tunicates to the adult multi-chambered-valvular of vertebrates. Tunicate, amphioxus, and lamprey images adapted from [6], [7], [8], respectively.
heart [1], [9], which is composed of only a single layer of myocardial cells. Next on the evolutionary chain is the amphioxus. The amphioxus heart is a rostrocaudally extended tube from its infancy through adulthood [10]. Similar to the tunicate heart, an amphioxus heart consists only of a monolayer of myocardial cells. Furthermore its heart has no chambers, valves, endocardium, epicardium, or other differentiated features of vertebrate hearts. Still, the amphioxus is regarded as the closest living invertebrate relative to vertebrates [11] and appears fish-like.
Furthermore, Figure 1 illustrates the bifurcation to multi-chambered hearts in a vertebrate -the lamprey. Cephalochordates, like amphioxus, have a series of four peristaltic vessels that serve as a pump, while tunicates have a single-chamber pump, which is composed of a single layer of myocardium (red) surrounded by stiff pericardial layer (pink). The earliest vertebrates, e.g., lampreys, have at least a two-chambered myocardium composed of a layer of cardiac myocardial cells (red), an endocardial cellular layer (yellow), valves that separate distinct chambers, and a surrounding pericardium (pink). Figure adapted from [16].
Lampreys are jawless fishes that are a very ancient lineage of vertebrates [12]. The lamprey is considered to have four heart chambers, which are the sinus venosus, atrium, ventricle, and conus arteriosus [13]. This is similar as to the zebrafish heart, which contains four chambers -the sinus venosus, atrium, ventricle, and bulbus arteriosus. Lamprey hearts also are valvular pumping systems, containing valve leaflets between chambers [14]. Moreover, lampreys are the first organism to develop an endocardial layer in addition to a myocardium, as well as, the first organism to develop cardiac valves [15]. An evolutionary depiction of heart morphology is illustrated in Figure 2, which was adapted from [16].
However, as discussed, the vertebrate embryonic heart is a valveless tube, similar to those in various invertebrates, such as urochordates and cephalochordates [17], [18], making invertebrates like sea squirts a possible model for heart development [19]. Historically, the pumping mechanism in these hearts has been described as peristalsis [17], [20], while more recently, dynamic suction pumping (DSP) has been proposed as a novel cardiac pumping mechanism for the vertebrate embryonic heart by Kenner et. al. in 2000 [21], and later declared the main pumping mechanism in vertebrate embryonic hearts by Forouhar et. al. in 2006 [22]. Debate over which is the actual pumping mechanism of the embryonic heart continues today, with the possibility that the mechanism may vary between species or may be some hybrid of both mechanisms [23], [24].
The Liebau pump, a dynamic suction pump, was first described in 1954 [25], and was studied as a novel way to pump water. It has not been until the past 20 years that scientists started looking at the pump as a valveless pumping mechanism in many biological systems and biomedical applications, including microelectromechanical systems (MEMs) and micro-fluidic devices. Direct applications of such pumps include tissue engineering, implantable micro electrodes, and drug delivery [26], [27], [28], [25].
With extensive industrial applications, dynamic suction pumping has proven to be a suitable means of transport for fluids and other materials in a valveless system, for scales of W o > 1 [29]. DSP can be most simply described by an isolated region of actuation, located asymmetrically along a flexible tube with stiffer ends. Flexibility of the tube is required to allow passive elastic traveling waves, which augment bulk transport throughout the system. The rigid ends of the tube cause the elastic waves to reflect and continue to propagate in the opposite direction, which when coupled with an asymmetric actuation point, can promote unidirectional flow. DSP is illustrated in Figure 3. Due to a coupling between the system's geometry, ma-terial properties of the tube wall, and pumping mechanics, there is a complex, nonlinear relationship between volumetric flow rate and pumping frequency [29], [30], [31]. Analytic models of DSP have been developed to address this relationship [32], [33], [34], [35], [30], [36]. Most models use simplifications, such as an inviscid assumption, long wave approximation, small contraction amplitude, and/or one-dimensional flow. Furthermore, no analytical model has described flow reversals, which can occur with changes in the pumping frequency. Relaxing many of these assumptions, physical experiments have been performed to better understand DSP [31], [37], [30], [25], as well as in silico investigations [38], [39], [40], [29], [41], [42]. Most of the joint experimental and computational studies focus on the 'high' W o regime (W o >> 1), besides studies by Baird et al. [41].  [43]. Spherical blood cells are seen within the tubular heart. The heart tube is roughly 5 blood cells thick in diameter.
In this paper we will investigate the pumping phenomena that occurs as a result of a coupled fluid-structureelectrophysiology model [44], and the bulk flow rates thereby induced. The electrophysiology model governs the propagation of action potentials, which then are coupled to muscular contraction, and hence the overall pumping dynamics. We then perturb the diffusion parameter in the electrophysiology model to investigate the bifurcation in pumping dynamics that occurs as a result of differing action potential propagation. This bifurcation is able to capture a spectrum of different pumping regimes, i.e., dynamic suction pumping to peristaltic-like waves of contraction. The electrophysiology model is governed by the FitzHugh-Nagumo equations [45], [44].
The power of this method is that it can be used to describe flow around complicated time-dependent geometries using a regular Cartesian discretization of the fluid domain. The elastic fibers describing the structure are discretized on a moving curvilinear mesh defined in the Lagrangian frame. The fluid and elastic fibers constitute a coupled system, in which the structure moves at the local fluid velocity and the structure applies a singular force of delta-layered thickness to the fluid.

A. Equations of the IBM
Assume that the immersed boundary is described on a curvilinear, Lagrangian mesh, S, that is free to move. The fluid is described on a fixed Cartesian, Eulerian grid, Ω, that has periodic boundary conditions. Given the size of the domain and the localization of the flow to the tube, the boundary conditions do not significantly affect the fluid motion. The governing equations for the fluid, the Navier-Stokes equations, are given by Eqs. (1) and (2) are the Navier-Stokes equations written in Eulerian form, where Eq.(1) is the conservation of momentum for a fluid and Eq.(2) is the conservation of mass, i.e., incompressibility condition. The two constant parameters in these equations are the fluid density, ρ, and the dynamic viscosity of the fluid, µ. The fluid velocity, u(x, t), pressure, p(x, t), and body force, f (x, t), are unknown spatial time-dependent functions of the Eulerian coordinate, x, and time, t. The body force describes the transfer of momentum onto the fluid due to the restoring forces arising from deformations of the elastic structure. It is this term, f (x, t), that is unique to the particular model being studied. The material properties of the structure may be modeled to resist to bending, stretching, and displacement from a tethered position. Other forces that can have been modeled include the action of virtual muscles, electrostatic (contact) forces, molecular bonds, porosity, and other external forces [46], [57], [58], [59], [56], [48]. The immersed structure may deform due to bending forces and/or stretching and compression forces. In this paper, elastic forces are calculated as beams that may undergo large deformations and Hookean springs, i.e., Eq. (3) is the beam equation, which describes forces arising from bending of the elastic structure. Eq.(4) describes the force generated from stretching and compression of the structure. The parameters, k beam and k spring , are the stiffness coefficients of the beam and spring, respectively, and R L is the resting length of the Hookean spring. The variables X M and X S give the positions in Cartesian coordinates of the master and slave nodes in the spring formulations, respectively, X B (s) describes the deviation from the preferred curvature of the structure. In all simulations, X B (s) = 0 along the straight portion of the tube. A target point formulation can be used to tether the structure or subset thereof in place, holding the Lagrangian mesh in a preferred position that may be time dependent. An immersed boundary point with position X(s, t) that is tethered to a target point, with position Y(s, t) undergoes a penalty force that is proportional to the displacement between them. The force that results is given by the equation for a linear spring with zero resting length, where k target is the stiffness coefficient of the target point springs. k target can be varied to control the deviation allowed between the actual location of the boundary and its preferred position. The total deformation force that will be applied to the fluid is a sum of the above forces, A more detailed description of existing fiber models can be found in [48]. Once the total force from Eq. (6) has been calculated, it needs to be spread from the Lagrangian frame to the Eulerian grid. This is achieved through an integral transform with a delta function kernel, Similarly, to interpolate the local fluid velocity onto the Lagrangian mesh, the same delta function transform is used, Eqs. (7) and (8) describe the coupling between the immersed boundary and the fluid, e.g., the communication between the Lagrangian framework and Eulerian framework. The delta functions in these equations make up the heart of the IBM, as they are used to spread and interpolate dynamic quantities between the fluid grid and elastic structure, e.g., forces and velocity. The quantity X(s, t) gives the position in Cartesian coordinates of the elastic structure at local material point, s, and time t.
In approximating these integral transforms, a discretized and regularized delta function, δ h (x) [46], is used, where φ(r) is defined as

B. Numerical Algorithm
As stated above, we impose periodic boundary conditions on the rectangular domain. To solve Eqs. (1), (2), (7) and (8) we need to update the velocity, pressure, position of the boundary, and force acting on the boundary at time n + 1 using data from time n. IBM does this in the following steps [46]: Step 1: Find the force density, F n on the immersed boundary, from the current boundary configuration, X n .
Step 2: Use Eq. (7) to spread this boundary force from the curvilinear mesh to nearby fluid lattice points.
Step 3: Solve the Navier-Stokes equations, Eqs. (1) and (2), on the Eulerian domain. In doing so, we are updating u n+1 and p n+1 from u n and f n . Note: because of the periodic boundary conditions on our computational domain, we can easily use the Fast Fourier Transform (FFT) [60], [61], to solve for these updates at an accelerated rate.
The above steps outline the process used by the IBM to update the positions and velocities of both the fluid and elastic structure. A more detailed discussion of the IBM can be found in [46].

C. Computational Model
We numerically model a 2D closed racetrack where the walls of the tube are modeled as 1D fibers. The closed tube is composed of two straight portions, of equal length, connected by two half circles, of equal inner and equal outer radii. The tube, or racetrack, has uniform diameter throughout. This is similar to the racetrack model geometry as in [44]. Furthermore, as in [44], we include the presence of an idealized stiff pericardium surrounding the flexible region of the heart tube.
The tunicate heart consists of a myocardium which is surrounded by a stiff pericardium [62], [63], which provides structural support to the myocardium. Muscle fibers spiral around the heart tube itself, and action potentials propagate to induce myocardial contraction. These action potentials have been previously measured [17]. Myocardial contractions may begin at either end of the heart tube, allowing the propagation of the action potential to occur in either direction [64]. However, we do not concern ourselves with flow reversals in this model. Although the tunicate heart tube has different material properties and physiological properties than the vertebrate embryonic heart, it still is an interesting model for vertebrate heart morphogenesis [19]. However, the conduction properties, e.g., velocities, of action potentials are much more uniform in tunicates than mammalian hearts [65]. The computational model we investigate is shown in Figure 5. Linear springs and beams connect adjacent Lagrangian points in the flexible region of the racetrack geometry. All other Lagrangian points of the boundary are modeled using target points, to hold the stiff portions of the racetrack and pericardium region nearly rigid. The parameters used in the model are found in Table I   Instead of prescribing contraction, we develop a model for the underlying electrophysiology of the heart, i.e., traveling action potentials arising from a single pacemaker region, to couple to myocardial contraction and hence intracardiac fluid flow. The model of action potential propagation is given by the FitzHugh-Nagumo equations [45], [44] below, where v(x, t) is the membrane potential, w(x, t) is the blocking mechanism, D is the diffusion rate of the membrane potential, v a is the threshold potential, γ is the resetting rate, is the blocking strength parameter, and I(t) is an applied current, e.g., an initial stimulus potentially from pacemaker signal activation. Note that v is the action potential and that w can be thought to model a sodium blocking channel. We note that the FitzHugh-Nagumo equations (11)- (12) are a reduced order model of the Hodgkin-Huxley equations, which were the first quantitative model to describe the propagation of an electrical signal across excitable cells [66]. The parameters used in the electrophysiology model are found in Table  II.
Next we need to interpolate the information from the electrophysiology model to the fluid-structure interaction solver, i.e., immersed boundary method. Time is scaled in order to match the dynamics of the generated action potentials to the desired active wave of contraction and is given by:  where dt is the time-step associated with the fluid solver, F is a non-dimensional scaling parameter, and T is the desired pumping period. The spatial location, x, in (11)- (12) is also scaled to match the dynamics of the active wave of contraction on the tube. When the propagating action potential reaches one of the muscles along the tube, the associated spring stiffness of said muscle model is given by The simplified muscle model is given by a dynamic spring stiffness coefficient, given by k e (x, t), which is a non-linear function of the traveling action potential, v(x, t). This idea was adapted from Baird et al. [41], [44]. The force generated by the springs that connect the bottom and top of the elastic tube can then be computed. These forces represents muscular contraction. The value of k m is tuned to produce the amount of contraction observed in Ciona hearts, as in [41], [44]. The idea for electro-dynamic pumping can be seen in Figure 6, which is a schematic for electro-dynamical pumping behavior. First the tube is at rest until a pacemaker initiates a potential signal, which contracts the tube in one singular region. Next the action potential propagates along the tube inducing contraction. Once the action potential passes outside a region on the tube, that location no longer has active contraction, but can return to its resting position. Furthermore the main electrophysiology idea behind the model is illustrated in Figure 7. In diagram 1 the flexible tube is at rest. Next 2 depicts a pacemaker initiating an input signal (current). Then that voltage (action potential) travels down the tube, while the input signal dissipates. Once the action potential reaches a muscle fiber, the tube contracts based on a non-linear relationship between spring stiffness and the magnitude of the action potential (voltage).

III. RESULTS
In this study, we conducted numerical experiments of the electro-dynamic pumping model, which encompassed fully coupled electrophysiology to pumping behavior for a heart tube, modeled as a closed racetrack geometry. We investigated various diffusivities, D, which give rise to different pumping regimes, e.g., either a 'dynamic suction pumping-esque' or 'peristalticlike' pumping regime. Furthermore, we explored these regimes for over 3 orders of magnitude in W o.

A. Results of the FitzHugh-Nagumo Model
Here we present the varying action potential dynamics given via the FitzHugh-Nagumo equation, which models the electrophysiology. We explored this model for a variety of diffusive coefficients, D = {0.1, 1.0, 10.0, 100.0}.  Figure 8 illustrates the kinds of traveling action potentials that arise out of the electrophysiology model. It is clear that the D = 0.1 case resembles a signal that could be reminiscent of that of dynamic suction pumping, where as D gives rise to a propagating action potential that could model a more peristaltic-like contraction. It is clear that as diffusivity increases, the waves propagate outwards, and with greater wave-speed. Furthermore, the wave-form itself gets wider.

B. Results of the electro-dynamical heart tube model
In this section we present the results describing how bulk flow rates are affected by varying the diffusivity, to capture different pumping behaviors for a variety of W o.  Furthermore, the wave-form in the D = 100 case is different between the W o = 0.1 and W o = 10 cases. There is a single peak for the case when W o = 10 and a dual peaks for W o = 0.1 for the forward flow; however, in the backflow, the situation is reversed, where a dualpeak is observed for W o = 10 and a single peak for W o = 0.1.
In attempt to maximize bulk flow for the dynamic suction pumping-esque regime, the stretching-stiffness and bending stiffness coefficients of the tube were varied. The results are shown in Figure 12. It is clear that as the stiffness is varied there is a non-linear relationships between flow speed (spatially-and temporally-averaged non-dimensional velocity across a cross-section of the racetrack) and stiffness. However, not a considerable amount of more bulk flow was produced from increasing these stiffness coefficients.
Lastly we compared the spatially-and temporallyaveraged non-dimensional velocities across a crosssection of the racetrack against W o for a variety of D.
The results are shown in Figure 13. It is clear there is a non-linear relationship in average velocity and scale arising from this model of pumping in every pumping regime, given by D. Furthermore, the highest bulk flow rates were seen in the case of D = 100 for W o ∼ 0.8, which correspond to the W o around that of tunicate tubular hearts [29], [41].

IV. DISCUSSION AND CONCLUSION
This 2D model coupled the propagation of action potentials, given via the FitzHugh-Nagumo equations, to the force generation and myocardial contraction, given through a non-linear spring-like muscle model, to induce pumping behavior in a flexible tube, where the fully coupled fluid-structure interaction model was solved using the immersed boundary method. This model was first described in [44]. We explored the effect of perturbing a diffusive coefficient in the electrophysiology model to capture different pumping regimes.
It was clear that by varying this diffusive term, D, the model was able to recreate a spectrum of pumping mechanisms, ranging from one that in which the action potential remained localized and did not diffusive, i.e., a dynamic suction pumping-esque behavior, and one where the action potential diffused along the heart tube in as a more traveling wave, e.g., peristaltic-like active wave of contraction. Our model showed that when D was in the more peristaltic-like regime, i.e., D ∼ 100, more bulk flow was produced in the racetrack geometry, as compared to more negligible amounts from the dynamic suction pumping-esque regime, D ∼ 0.1. This result was consistent for the range of W o considered.
Moreover, in all cases considered, there was a nonlinear relationship between average flow rate, scale (W o), and diffusivity (pumping behavior). More bulk  flow was produced on average (both spatially and temporally), with a maximum around W o ∼ 0.8 than for higher W o, up to W o = 30, in the peristalic-like regime. Similar behavior, in that peristalsis produces more bulk flow than DSP, has been observed before when using prescribed pumping behavior, as in [29], [42].
However, perturbing the material properties of the tube could potentially affect bulk flow rates across all pumping regimes, given by D. Our focus was limited to perturbing the stretching and bending stiffnesses of the tube specifically within the dynamic suction pumpingesque regime, D ∼ [0.1, 1]. Furthermore, our study only considered increasing the stiffnesses and not decreasing them. For the regime and material properties considered, we found a non-linear relationship between flow rates and stiffness.
As blood flow and the resulting hemodynamic forces are essential for proper heart development [67], it is important that the pumping model capture as much biology as possible. Each pumping regime considered here, given by the diffusivity of action potential propagation, will give rise to a different force distribution along the endothelial lining of the heart and hence impact the epigenetic signals that are transmitted through mechanotransduction [68], [69]. Furthermore the flow profiles resulting from each pumping mechanism would be different. These differences in the flow patterns itself could impact the way morphogens advect and diffuse during embryogenesis [70], [71], opening the realm to a lot more interesting biological questions to explore.