Anisotropic shear stress patterns predict the orientation of convergent tissue movements in the embryonic heart

Myocardial contractility and blood flow provide essential mechanical cues for the morphogenesis of the heart. In general, endothelial cells change their migratory behavior in response to shear stress patterns, according to flow directionality. Here, we assessed the impact of shear stress patterns and flow directionality on the behavior of endocardial cells, the specialized endothelial cells of the heart. At the early stages of zebrafish heart valve formation, we show that endocardial cells are converging to the valve-forming area and that this behavior depends upon mechanical forces. Quantitative live imaging and mathematical modeling allow us to correlate this tissue convergence with the underlying flow forces. We predict that tissue convergence is associated with the direction of the mean wall shear stress and of the gradient of harmonic phase-averaged shear stresses, which surprisingly do not match the overall direction of the flow. This contrasts with the usual role of flow directionality in vascular development and suggests that the full spatial and temporal complexity of the wall shear stress should be taken into account when studying endothelial cell responses to flow in vivo.

AVC unfolding, quantification of the AVC shortening, and statistical analysis The threedimensional movement of the photoconverted cells within each heart was analyzed using a newly developed approach.
The following steps are implemented in Matlab scripts.
• Step 1: Points on the endocardium are identified by intensity thresholding of the acquired images and plotted in three dimensions. Points on the atrium and the ventricle are erased interactively using the erase command of Matlab to isolate points representing the narrowest and most straight part of the AVC. These points are selected and fitted with a cylinder as described previously (McMahon et al., 2008). A reference system can then be introduced with the axis of the cylinder as z-axis and the center of mass of the AVC points as origin. The endocardial points are reoriented such that the points on the ventricle have a positive z-(axial) position. The heart points are then rotated around the z-axis with this strategy: The AVC points are projected on the x-y plane and fitted with an ellipse. The endocardial points are then rotated around the z-axis to bring the major ellipse axis parallel to the y-axis and such that the internal side of the AVC (with respect to the embryo) points has positive y-values. At this point, the angular position of cells in the AVC is defined consistently for all the hearts.
• Step 2: The three-dimensional dataset is cut by planes obtained by rotating the half-plane {y=0, x>0} in the z-direction. The intensity of the neighboring pixels are projected on these planes and the endothelium traced with a line of points ordered with increasing arc-length from the atrium to the ventricular side of the AVC. These points are then interpolated by a three dimensional spline (spline toolbox, matlab). Each point on the resulting surface have a well defined angular position and the tissue length between two points at the same angular position is computed as the archlenght on the surface between the two points. This is different than the length of a linear segment between two points because the endocardium is curved.
• Step 3: The intensity of the neighboring pixels are projected on the parametric surface.
This gives a two-dimensional image with each column corresponding to the line on the surface at a given angle. This is equivalent to unfold the AVC into a two-dimensional image.
• Step 4: The unfolded AVC can then be segmented in two-dimensions with standard tools to find the edges of the non-photoconverted AVC and the AVC length L as a function of the angular position.
The same analysis was repeated at 36 and 48 hpf for each heart considered, and the AVC shrinking factor was computed as (L 48 − L 36 )/L 36 as a function of the angular position and averaged over the inferior, superior, internal and external regions of the AVC (Fig. 2), which are defined consistently based on the orientation of the elliptic cross-section of the AVC and its orientation within the embryo.
The statistical significance of the differences between the mean values calculated for control and perturbed flow conditions was determined by unpaired Student's t-tests and computed with the ttest2 function of the Matlab Statistical Toolbox, without assuming equal variances.
Shown standard deviations were computed as corrected sample standard deviations using the std Matlab function.

Computational Methods
Heart dynamics and anatomy We exploit the viscous regime of the zebrafish heart, where the Reynolds number is reported to vary from 0.017 at 26hpf to 0.342 at 4.5 days post fertilization (Santhanakrishnan and Miller, 2011). Therefore we neglect inertia effects, though they might alter details at 4.5dpf, as well as the overall curvature of the heart which will not alter flow substantively at so low of Reynolds numbers. In two dimensions, the dynamics of a point on the wall is simply defined as: y(x, t) = d 0 + A sin(t/f 1 2π + θ 1 ) for both the upper and lower walls, where we truncate the Fourier expansion of d(t) to the heart beat frequency f 1 . The time average radius d 0 and the oscillation amplitude A and the beating frequency f 1 were previously extracted from live imaging of 48hpf zebrafish hearts as described previously (Boselli and Vermot, 2016). Previous works showed how two-dimensional models can well predict typical shear stress in the heart of zebrafish embryos (Lee et al., 2013;Boselli and Vermot, 2016), which is the focus of this work. Therefore, despite the fact that an accurate description of some hemodynamic cues, including pressure, would require a three-dimensional model of the heart and of the hydrodynamic load represented by the vascular system, we will simplify the system and limit our analysis to two dimensions.
where S ij is the Stokeslet tensor, s is the arc-length position of a point x(s) on the membrane, and ∆σ j is the traction of the membrane on the fluid. For the cell membrane, the traction is given by the relation: with τ the membrane tension, and b the bending moment, with C the curvature and C 0 a constant initial curvature that cancels out in (2). The coordinate s 0 in Eqations (3) and (4) is the referential arch-length position of a point on the cell membrane in an hypothetic stress-free configuration. In practice, the stress-free shape of the RBCs does not need to be defined. It is sufficient to set the perimeter l 0 of the stress-free membrane and the area r 2 0 π of each of our two-dimensional RBCs such that l 0 = 2πr 0 . Starting from an arbitrary initial configuration (here a sphere of perimeter 2πr 0 ), the RBC will assume a configuration that minimize the elastic energy of the membrane. The time required by this process scales like the relaxation time of the red blood cell τ rlx = r 0 µ/T , and its ratio with the convective time scale τ U = r 0 /U gives a non dimensional index of the elasticity of the red blood cells (Freund, 2007). The main simulations were set-up such that τ rlx /τ U = 0.16, r 2 0 T /M = 12.5 and l 0 /(2πr 0 ) = 1.6π. Simulations presented in Fig. S3, where repeated for l 0 /(2πr 0 ) = 1.2π, τ rlx /τ U = 0.1, and τ rlx /τ U = 0.32, which, in the order, correspond to rounder, stiffer and softer red blood cells. The research code is available upon request.
Wall model The heart walls dynamics is obtained by imposing the wall traction which is equivalent to link each point on the wall to their desired location x w by a virtual spring of stiffness k w , which was kept constant and such that k w r 2 0 /T = 1.5 in the main simulations.
The flow velocity due to the wall traction can be expressed by the same boundary internal formulation (1) used for the cell membrane, but with ∆σ w i instead of ∆σ m i . As for the RBCs, the number of points necessary to describe the shape of the wall is smaller than that required to compute the integrals along the wall. Only these control points are carried on during time integration. The extra collocation points for the space integrals are then obtained by Fourier interpolation, which exploits the periodicity of the boundary.
Time integration A point x = {x i } on the cell membrane or on the wall moves then according to where u i now is the sum of the contribution of all the cell membranes and the walls in the model system. Time integration is performed by a second order Runge-Kutta scheme. The time step dt was set such that dt * f 1 = 5 · 10 −5 .
Numerical discretization The boundary integral (1)  where N a is the number of points on the cell membrane and wall used for the quadrature and N = N m + N w + N v is the number of evaluation points, which is the sum of N m collocation points used to represent the cell membrane, N w wall points, and N v points where to evaluate u i for post processing or simply for visualization.
The numerical discretization, including the PME method employed to speed up the evaluation of (7), is described in details in (Freund, 2007). Starting from the results in (Freund, 2007), the collocation points on the membrane of each RBCs was set to 128, and to 800 on the wall.
Computation of the shear stress The integral boundary methods extend the flow domain over the whole space such that flow is also computed outside the actual domain of interest.
Therefore, the wall shear stress at the endocardial interface is computed as where t w i is the unit tangent vector to the wall and τ out is the wall shear stress due to the flow outside the heart model. This is computed as to minimize σ w i t w i − τ out − τ in without RBCs. τ in is not computed directly to avoid the finite difference kernel to overlap with an RBC passing very close to the wall.The research code is available upon request.

Supplementary figures and movies
Supplements include one movie: Movie 1; there supplementary figures: Fig. S1, S2, 3; and Table S1.    Figure S2: Shear stress pattern during tissue convergence (higher harmonics) in complement to Fig. 3. A,C) Second (τ 2 ) and B,D) third (τ 3 ) harmonics. The results in A and B are for 106 red blood cells. The results on C and D are for the upper wall and for different numbers of red blood cells n p =1, 17, 54,106 (arrows point to larger values of n p ). The superscript * denotes that results are normalized by the space average of t f 1 |τ |dt. a: atrium; avc: atrioventricular canal; v: ventricle. Figure S3: Model sensitivity to red blood cell (RBC) shape and stiffness. Related to Fig. 3. A) Computational model for RBC flows at three different instants of the heart beat for round elastic red blood cells (cf. Fig. S1G-I). B) Flow direction in the AVC model: magenta, flow is from atrium to ventricle; black, flow is from ventricle to atrium; white, no flow or not visible. Solid arrows point in the direction of the flow. C) Fundamental harmonic τ 1 , D) time average τ 0 , and E) non-periodic oscillationτ of the shear stress. C-E) The blue solid line corresponds to the RBC shape of the test case of Fig. S1 and 3, n p =106, and τ rlx /τ U = 0.16. The three dashed lines correspond to the rounder RBC shape shown in A and are obtained using the same, stiffer (τ rlx /τ U = 0.1), and softer (τ rlx /τ U = 0.32) RBC membranes (dashed arrow points to larger value of stiffness). The superscript * denotes that results are normalized by the space average of t f 1 |τ |dt. a: atrium; avc: atrioventricular canal; v: ventricle.