蔡尚炅 Shang-Gui Cai
email: shanggui.cai@univ-amu.fr

Abdellatif Ouahsine

Yannick Hoarau

Moving immersed boundary method for fluid-solid interaction

A strongly coupled algorithm is presented for the incompressible fluid-rigid body interaction using the moving immersed boundary method. The pressure and the boundary force are treated as Lagrange multipliers to enforce the incompressibility and no-slip wall constraints. To compute the two unknowns from the velocity field, we adopt the fractional step algorithm and successively apply the two constraints. A Poisson equation and a moving force equation are derived for the pressure and the boundary force respectively. As both coefficient matrices are formulated to be symmetric and positive-definite, the resulting linear systems are solved efficiently with the conjugate gradient solver. The stronglycoupled nonlinear fluid-solid system is achieved by a fixed point iteration. To improve the computational efficiency, we only iterate the moving force equation with the rigid body motion equation and the time-consuming pressure Poisson equation is solved once the sub-iteration is finished. The proposed method is validated with various benchmark tests and the results compare well with the literature.

I. INTRODUCTION

One primary challenge of fluid-solid interaction (FSI) simulation is dealing of complicated and time-dependent interface between two physical fields. Body-fitted mesh methods, such as finite volume or finite element discretization of fluid equations in arbitrary Lagrangian-Eulerian formulation, need mesh-deforming or re-meshing regularly to accommodate the time-varying solid boundaries. It is very costly for large scale problems and prone to mesh interpolation errors and divergence due to severe mesh distortion. The Chimera/overset grid method 1 circumvents these issues by employing multiple computational fluid domains covered by independent meshes. One background mesh, usually a Cartesian mesh, is utilized to discretize the global domain, while the others are body-fitted meshes surrounding the objects. The Chimera/overset grid method is capable to capture the boundary layer well at moderate and high Reynolds number, comparable to classical body-fitted mesh methods, since grids can be allocated efficiently when aligning to solid boundaries. However constructing a body-fitted mesh of high quality for complex geometries is still non-trivial. Simple discretization method like finite difference may require a coordinate transformation to work on a curvilinear mesh. In such circumstance finite volume or finite element method is usually used and the resulting system is solved with sophisticated solvers, such as preconditioned Krylov subspace methods (e.g. conjugate gradient) or algebraic multigrid method.

Moreover since the fluid equations are computed on different domains, it is difficult to ensure conservation.

Alternatively, the immersed boundary method (IBM) [START_REF] Peskin | Flow patterns around heart valves: A numerical method[END_REF][START_REF] Uhlmann | An immersed boundary method with direct forcing for the simulation of particulate flows[END_REF][START_REF] Kempe | An improved immersed boundary mehod with direct forcing for the simulation of particle laden flows[END_REF][START_REF] Cai | Computational fluid-structure interaction with the moving immersed boundary method[END_REF][START_REF] Cai | Improved implicit immersed boundary method via operator splitting[END_REF][START_REF] Cai | Implicit immersed boundary method for fluid-structure interaction[END_REF][START_REF] Cai | Moving immersed boundary method[END_REF] adopts a surface mesh for the solid boundary, thereby the fluid is only simulated on the background mesh. Using a simple fluid mesh is advantageous not only for the ease of mesh generation but also for the accuracy matter. The Cartesian mesh produces the best mesh quality and hence high order accuracy can be achieved. The matrix assembly is also simplified and the matrix condition is well preserved. Most importantly, it facilitates the use of many simple and efficient methods developed for Cartesian mesh. The finite difference discretization can be constructed into high order accuracy easily and highly efficient solvers, such as geometrical multigrid method and fast Poisson solver, can be used straightforward. Parallelization of program is relatively simple to a good efficiency. However as no underlying mesh will conform to the solid boundary, the imposition of boundary condition at the interface can be tricky. Interpolation methods are often selected to overcome this problem, either by locally reconstructing the velocity [START_REF] Cai | Coupling of turbulence wall models and immersed boundaries on Cartesian grids[END_REF][START_REF] Degrigny | Improved wall model treatment for aerodynamic flows in LBM[END_REF] or by adopting a source forcing function in the momentum equation [START_REF] Cai | Computational fluid-structure interaction with the moving immersed boundary method[END_REF][START_REF] Cai | Improved implicit immersed boundary method via operator splitting[END_REF][START_REF] Cai | Implicit immersed boundary method for fluid-structure interaction[END_REF][START_REF] Cai | Moving immersed boundary method[END_REF] . In fact, the two approaches are essentially the same as they both rectify the velocity in the vicinity of immersed boundary to produce solid effect [START_REF] Tseng | A ghost-cell immersed boundary method for flow in complex geometry[END_REF] .

For the source forcing approach, it is crucial to compute the boundary force accurately to ensure a physical boundary condition. Uhlmann 3 developed an efficient direct forcing IBM for simulating particle suspensions, where an explicit formulation is used for the boundary force. As a result, the no-slip wall condition is not satisfied and large errors occur at the interface. Kempe and Fröhlich 4 reduced the error through an additional forcing loop, which is solved with only a few number of iterations without convergence. Therefore the no-slip wall condition will not be fully verified. Full convergence however requires considerable computational efforts, thus it is often relaxed in practice [START_REF] Breugem | A second-order accurate immersed boundary method for fully resolved simulations of particle-laden flows[END_REF] . Actually, the boundary force can be viewed as a Lagrange multiplier for the no-slip constraint. an implicit immersed boundary projection method (IBPM) is proposed by Taira and Colonius 13 by integrating the boundary force and the pressure into a modified Poisson equation. The IBPM is very accurate since it enforces the divergence-free and the no-slip constraints simultaneously. In spite of the mathematical rigor and completeness, the IBPM is less efficient and stable than the explicit IBMs, as an enlarged system is formulated and evolves in time for moving bodies. Instead of solving the pressure and the boundary force in one single system, the moving immersed boundary method (MIBM) was proposed by Cai et al. [START_REF] Cai | Computational fluid-structure interaction with the moving immersed boundary method[END_REF][START_REF] Cai | Improved implicit immersed boundary method via operator splitting[END_REF][START_REF] Cai | Implicit immersed boundary method for fluid-structure interaction[END_REF][START_REF] Cai | Moving immersed boundary method[END_REF] to decouple the pressure and the boundary force into two independent systems. The pressure Poisson coefficient matrix is kept unchanged and hence preconditioners only need to be constructed in the initial step.

The boundary force is computed from a moving force equation after operator splitting. The size of the corresponding matrix is considerably reduced as it is only proportional to the number of Lagrangian markers defining the immersed boundary, which is much smaller than the fluid mesh numbers. The MIBM is found very efficient and accurate for moving boundary problems. In this work, we will extend the MIBM to fluid-solid interaction applications by incorporating the rigid body dynamics.

Generally, two approaches can be used for the fluid-solid interaction time coupling: the monolithic (direct) approach and the partitioned (segregated) approach. For the monolithic method, all the unknowns from the fluid and solid equations are formulated into a single large system along with interface conditions. While the partitioned method formulates the two subproblems into different systems and solve them separately. The monolithic approach is advantageous from the stability point of view, but it requires considerable efforts on code modifications. The partitioned approach can be further classified into two subgroups, namely the explicit coupling (weakly or loosely coupled) and the implicit coupling (strongly or tightly coupled). Weak coupling is very efficient since it solves the two subproblems in a sequential manner, but it leaves the interface conditions unfulfilled. As a result, spurious energy is generated at the interface, which leads to unstable solutions, which is often called the added mass effect [START_REF] Causin | Added-mass effect in the design of partitioned algorithms for fluid-structure problems[END_REF][START_REF] Förster | Artificial added mass instabilities in sequential staggered coupling of nonlinear structures and incompressible viscous flows[END_REF] . The stability criterion mainly depends on the density ratio between the solid and the fluid. Below a certain threshold the calculation diverges immediately for any chosen time step for incompressible fluids [START_REF] Kassiotis | Nonlinear fluid-structure interaction problem. Part I: implicit partitioned algorithm, nonlinear stability proof and validation examples[END_REF] . The strongly coupled method iterates the two subproblems until the interface conditions are satisfied, thus it is comparable to the monolithic method. In the present work, we will design a partitioned but strongly coupled scheme with the moving immersed boundary method.

The organization of this paper is as follows. First we present the general fluid-solid interaction problem in the immersed boundary formulation. In the following part the MIBM and the strongly coupled partitioned scheme are presented in detail. Numerical simulations are then performed with the proposed method and validated with benchmark tests. Conclusions are drawn in the final section. where the fluid and the rigid body occupy the physical domains Ω f and Ω s respectively. The fluidsolid interaction takes place at the their common boundary ∂ Ω i = Ω f ∩ Ω s . The whole system is subjected to the gravitational acceleration g.

II. METHODOLOGY

Ω f ∂ Ω f Ω s r ∂ Ω i g v s ω s O x y x s x c
FIG. 1: Schematic representation of the fluid-solid interaction problem.

The fluid subproblem is usually described in Eulerian coordinates, which is governed by the incompressible viscous Navier-Stokes equations

∂ v f ∂t + ∇ • (v f ⊗ v f ) = ∇ • σ f + g, (1a)
∇ • v f = 0, ( 1b
)
where v f represents the fluid velocity vector and the fluid stress tensor σ f (divided by ρ f ) is given by

σ f = - p ρ f I + ν ∇v f + (∇v f ) T , (1c)
where p designates the fluid pressure, ρ f the fluid density, ν the fluid kinematic viscosity. Appropriate initial and boundary conditions are assumed to the fluid Navier-Stokes equations to ensure that the problem is well posed.

Contrarily the solid subproblem is often treated in a Lagrangian framework. The motion of a rigid body is governed by the Newton-Euler equations, which can be expressed as

m s dv s dt = ρ f ∂ Ω i σ f • nds + m s (1 - ρ f ρ s )g, (2a)
I s dω s dt = ρ f ∂ Ω i r × σ f • n ds, (2b)
where m s , ρ s , I s represent the solid mass, density and the moment of inertia respectively. v s , ω s designate the translational velocity and the angular velocity of the solid. r = x s -x c is the position vector of any surface point with respect to the solid mass center x c . n represents the outwardpointing normal vector to the surface ∂ Ω i . The position of the rigid body can be obtained by integrating the kinematic equations

dx c dt = v s , (2c)
dθ c dt = ω s , (2d)
where θ c designates the rotation angle around the solid mass center.

On the fluid-solid interface ∂ Ω i the continuity of velocity

v f = v s + ω s × r, (3)
needs to be satisfied in order to take the fluid-solid interaction into account.

The immersed boundary method approximates the above fluid-solid interaction problem by replacing the solid domain with the surrounding fluid. To account for the presence of the immersed solid, a boundary force f is introduced and added into the fluid momentum equation. Therefore the fluid is simply simulated in a fixed domain Ω = Ω f (t) ∪ Ω s (t) irrespective to the movement of the immersed solid. The fluid-solid interaction problem in the immersed boundary formulation becomes

∂ v f ∂t + ∇ • (v f ⊗ v f ) = - 1 ρ f ∇p + ν∇ 2 v f + f in Ω, (4a)
∇ • v f = 0 in Ω, (4b)
m s dv s dt = -ρ f Ω s fdV + m s (1 - ρ f ρ s )g, (4c)
I s dω s dt = -ρ f Ω s r × fdV, (4d)
v f = v s + ω s × r on ∂ Ω i , (4e)
f(x,t) = Γ s F(s,t)δ (x -X(s,t))ds, (4f)
where F(s,t) represents the boundary force defined on the Lagrangian position X(s,t) and f(x,t)

on the Eulerian frame respectively. δ designates the Dirac delta function. Note that the effect of gravity in the fluid momentum equation is from now on incorporated into the pressure.

The fluid governing equations (4) are discretized in time as

v n+1 f -v n f ∆t + 3 2 N (v n f ) - 1 2 N (v n-1 f ) = - 1 ρ f G p n+1 + ν 2 L (v n+1 f + v n f ) + S F n+1 , (5a)
Dv n+1 f = 0, (5b)
T v n+1 f = v n+1 s + ω n+1 s × r n+1 on ∂ Ω n+1 i , (5c)
where second order semi-implicit time stepping schemes are utilized for the velocity, namely the explicit Adams-Bashforth scheme for the convection and the implicit Crank-Nicholson for the diffusion. Euler implicit scheme is applied to the two Lagrange multipliers, that is, the pressure and the boundary force. N , L , D, G represent the convection, Laplacian, divergence, and gradient operators respectively. Since the fluid grid is not coincident with the solid node, T and S are the interpolation and spreading operators that are responsible for the field exchange between two grids (see Figure 2). The discrete delta function is employed for the construction of T and S

δ h (x -X) = 1 h 2 φ ( x -X h )φ ( y -Y h ). (6)
In this work, we select the three-point kernel proposed in Ref. 17 aimed for staggered grid, whose one-dimensional form is With respect to the solid equation, we choose a simple Euler implicit scheme

φ (r) =              1 3 1 + -3r 2 + 1 , |r| < 0.5, 1 6 5 -3|r| --3(1 -|r|) 2 + 1 , 0.5 ⩽ |r| < 1.5, 0, otherwise. (7)
m s v n+1 s -v n s ∆t = -ρ f Ω s f n+1 dV + m s (1 - ρ f ρ s )g, (8a)
I s ω n+1 s -ω n s ∆t = -ρ f Ω s r × f n+1 dV, (8b)
x n+1 c -x n c ∆t = v n+1 s , (8c)
θ n+1 c -θ n c ∆t = ω n+1 s . (8d)
The total amount of force and torque are not changed during the transfer between Lagrangian and Eulerian locations, owing to the partition of unity of the discrete delta function. The solid equations can be solved as follows

m s v n+1 s -v n s ∆t = -ρ f n b ∑ l=1 F n+1 ∆V l + m s (1 - ρ f ρ s )g, (9a)
I s ω n+1 s -ω n s ∆t = -ρ f n b ∑ l=1 r × F n+1 ∆V l , (9b)
where ∆V l represents the area/volume associated with the Lagrangian marker. Therefore, provided the boundary force, the solid equations can be easily solved. For the fluid part, we employ an operator splitting scheme as

(1) Prediction step for vn+1

f vn+1 f -v n f ∆t + 3 2 N (v n f ) - 1 2 N (v n-1 f ) = - 1 ρ f G p n + ν 2 L (v n+1 f + v n f ). (10)
(2) Immersed boundary forcing step for imposing the no-slip wall condition at the interface

ṽn+1 f -vn+1 f ∆t = S F n+1 , (11a)
T ṽn+1 f = v n+1 s + ω n+1 s × r n+1 on ∂ Ω n+1 i . ( 11b
)
Applying the no-slip constraint, we obtain

T S F n+1 = v n+1 s + ω n+1 s × r n+1 -T vn+1 f ∆t , (12a) ṽn+1 f
= vn+1 f + ∆tS F n+1 . (12b)
Given M = T S the moving force coefficient matrix, then we obtain a compact form of the moving force equation

M F n+1 = F e , (13)
with

F e = v n+1 s + ω n+1 s × r n+1 -T ṽn+1 f ∆t , (14)
where F e is the explicit force used by Kempe and Fröhlich 4 . Given the solid velocity and position at time level n + 1, the moving force equation can be solved with various linear system solvers. In the present work the moving force coefficient matrix M is found to be symmetric and positivedefinite. Hence the moving force equation can be converged to the machine precision very quickly by using the conjugate gradient method, even without preconditioning. Moreover, the size of the moving force coefficient matrix is proportional to the number of nodes on the immersed boundary, which is in general much smaller than the size of the fluid matrix. Since the moving force coefficient matrix does not involve the fluid matrix, its update is computationally inexpensive in the case of moving boundaries.

(3) Projection step for obtaining a divergence free velocity

v n+1 f v n+1 f -ṽn+1 f ∆t = -G φ n+1 , (15a)
Dv n+1 f = 0. ( 15b
)
Applying the divergence operator to (15a) along with the divergence free condition (15b), the projection step is actually performed as

L φ n+1 = 1 ∆t D ṽn+1 f , (16a)
v n+1 f = ṽn+1 f -∆tG φ n+1 . ( 16b
)
The final pressure is advanced by

p n+1 = p n + φ n+1 - ν 2 D vn+1 f . (17)
Now we will present a novel coupling approach for the interaction of the two subproblems.

Mathematically speaking, the FSI problem of equation ( 4) is accomplished by the Lagrange multiplier method to joint two subdomains on a common boundary. The unknowns of the entire system are (v n+1 f , p n+1 ) for the fluid, (x n+1 s , ẋn+1 s ) for the solid, and F n+1 for the interface condition.

The conventional strongly coupled method iterates those variables at each time step in order to match the interface condition. Even if it leads to accurate and stable solutions, solving implicit coupling usually exhibits a prohibitive computational cost. This becomes more unaffordable in ALE formulation, since the mesh and the associated matrices have to be updated regularly.

Fernández et al. [START_REF] Fernández | A projection semi-implicit scheme for the coupling of an elastic structure with an incompressible fluid[END_REF] proposed an efficient coupling method by taking advantage of the projection method. The prediction (ALE-advection-viscous) step is ruled out from the FSI loop to reduce computational cost, so that the mesh and the associated matrices are computed once at each time step. Finally only the projection step is coupled implicitly to ensure stability. It has been shown that this semi-implicit scheme is stable for a reasonable range of the discretization parameters, compared to the explicit coupling approach.

We extend this idea to the moving immersed boundary method. However the projection step usually is the most time-consuming part in the projection method. In spite of various methods (e.g. Aitken relaxation) are available to accelerate the coupling procedure, the computational cost still remains high. We also notice that in equation ( 9) the solid is coupled to the fluid by the boundary force but not the pressure, and in equation ( 12) the boundary force is determined by the solid velocity and position. Therefore, we can move out the time-consuming projection step from the FSI iteration, and replace it by the immersed boundary forcing step. Following Ref. 18, the prediction step is not included in the FSI iteration either, as it is performed on a stationary combined domain in MIBM and the boundary force is not incorporated in this step for the solid effects. The moving force equation can be considered as an implicit equation for the no-slip boundary condition for the fluid at the interface. Therefore the implicit coupling of the immersed boundary forcing step with the solid equations features a strongly coupled FSI method.

Algorithm 1: Novel implicit coupling scheme and the pressure p n+1 using (16a), (16b), (17).

1 Given: v n f , p n , x n s , ẋn s ; 2 (
Note that the moving force equation is non-linear, because the interpolation, spreading op-erators and the boundary force are functions of solid position x n+1 s , namely M n+1 F n+1 = T (x n+1 s )S (x n+1 s )F(x n+1 s ). We can linearize this equation by treating the moving force coefficient matrix explicitly M (x n s ), but this will decrease the overall accuracy due to the time lag, as indicated in Ref. 19. In order to preserve a high accuracy, fully implicit implementation of the moving force equation is considered in the Algorithm 1.

III. NUMERICAL RESULTS

A. Freely falling cylinder in a confined channel

We first investigate the motion of a circular cylinder falling freely between two parallel walls in a quiescent fluid. The width and height of the computational domain are chosen to be 2 cm and 6 cm. The circular cylinder with a diameter of D = 0.25 cm is released initially from the position (1 cm, 4 cm) and falls down because of gravity. The density of the cylinder and the surrounding fluid are ρ s = 1.25 g/cm 3 and ρ f = 1.0 g/cm 3 respectively. The fluid dynamic viscosity µ is set to 0.1 g/(cm•s).

The calculations are carried out on two different uniform meshes, i.e. h = 1/48 cm and h = 1/96 cm, to check the mesh sensitivity. The immersed cylinder surface is discretized evenly with a mesh size approximated to the size of surrounding fluid cells, due to the inherent limitation of the discrete delta function. Therefore the resulting Lagrangian marker numbers are 38 and 76 respectively. The time step is chosen to be ∆t = 0.001 s and the resulting maximum CFL number is reported to be 0.46. The vorticity around the falling cylinder is shown in Figure 3

B. Freely falling and rising cylinder in an open domain

Next we consider an object freely falling and rising in an open domain as another test. This phenomenon happens frequently in nature and a large amount of work can be found in the literature. Here we compare to the numerical results of Namkoong et al. [START_REF] Namkoong | Numerical analysis of two-dimensional motion of a freely falling circular cylinder in an infinite fluid[END_REF] and Lacis et al. 19 , using the body-fitted ALE formulation and the immersed boundary projection method respectively.

Two density ratios are considered in this study, i.e. ρ s /ρ f = 1.01 for the falling case and ρ s /ρ f = 0.99 for the rising simulation. A large computational domain is taken as [-5D, 5D] × [-70D, 70D] with free-slip boundary conditions applied at all the exterior boundaries, where D = 0.5 cm is the cylinder diameter. A uniform mesh is employed to cover the computational domain, and the mesh resolution is kept to 0.04D in order to compare with Lacis et al. 19 . Initially the cylinder is located at ±65D, depending on the situation (65D for the falling case, -65D for the rising case). The Reynolds number Re = V t D/ν f here is 156, where V t is the terminal velocity.

Note that the Reynolds number depends on the Galileo number G = (|ρ s /ρ f -1|gD 3 ) 1/2 /ν f (the gravity force divided by the viscous force, G = 138) and the density ratio ρ s /ρ f . The instantaneous vorticity fields are presented in Figure 5 for the falling cylinder case. Initially symmetric vortex pair forms behind the cylinder in the beginning of falling. After that the numerical error accumulates and breaks the symmetry. At around tV t /D = 40, the flow becomes unsteady and periodic vortex shedding occurs. The time histories of the velocity of the cylinder are plotted in Figure 6. Table I shows the Strouhal number St = f D/V t ( f is the shedding frequency)

and the coefficients of drag and lift. Present results are compared to those of Namkoong et al. [START_REF] Namkoong | Numerical analysis of two-dimensional motion of a freely falling circular cylinder in an infinite fluid[END_REF] and Lacis et al. 19 . Good agreements have been obtained.

To illustrate the efficiency of the proposed coupling method, different strong FSI coupling of the FSI coupling time within nearly the same number of sub-iteration.

C. Elliptical particle sedimentation in a confined channel

In this example, we consider the sedimentation of an elliptical particle in a narrow channel, to demonstrate the ability of current FSI algorithm for non-circular object. This example was stud-modes of sedimentation have been found ranging from oscillating, tumbling along the wall, vertical sedimentation, horizontal sedimentation to an inclined mode. In their work, a multi-block lattice Boltzmann method is used and compared to the traditional ALE formulation. To compare with Xia et al. [START_REF] Xia | Flow patterns in the sedimentation of an elliptical particle[END_REF] , the computational domain is selected to be [0, L] × [0, 7L] with L = 0.4 cm, as shown in Figure 8. The aspect ratio of the ellipse is α = a/b = 2, where a and b are the major and minor axes respectively. The blockage ratio is defined as β = L/a = 4. The density ratio is ρ s /ρ f = 1.1. The kinematic viscosity of fluid is set to ν = 0.01 cm 2 /s. The particle starts falling in a quiescent fluid from the centroid at (0.5L, 6L) with an initial angle of π/4 to break the symmetry.

No-slip boundary conditions are applied at four boundaries. A uniform mesh is employed with a gird resolution of 0.0027 cm. The time step is chosen such that the CFL condition is satisfied.

Figure 9 shows the vorticity fields at different times at t = 0.1 s, 0.3 s, 0.5 s, 1.0 s, 1.5 s, 2.0 s. The trajectory and orientations are compared to the results of Xia et al. [START_REF] Xia | Flow patterns in the sedimentation of an elliptical particle[END_REF] in Figure 10. Good agreements have been obtained.

D. Flow around a rotating NACA0012 airfoil

The rigid objects simulated in this section so far have been the circular or elliptical particles.

In this example, we consider the incompressible viscous flow over a NACA0012 airfoil to demon-strate the capacity of current MIBM for handling sharp geometries. The shape of the NACA0012 airfoil is given by

Y = ±0.6 • (0.2969 √ X -0.1260X -0.3516X 2 + 0.2843X 3 -0.1015X 4 ), (18)
where X ∈ (0, 1.009) cm. The characteristic length, i.e. the airfoil length, is 1.009 cm. Following Glowinski et al. [START_REF] Glowinski | A fictitious domain approach to the direct numerical simulation of incompressible viscous flow past moving rigid bodies: application to particulate flow[END_REF] and Wan and Turek 20 , we select a computational domain of [-4 cm,16 cm] ×

[-2 cm, 2cm] with the airfoil centered at (0.42 cm,0), as shown in Figure 11. The airfoil is fixed at its mass center and is free to rotate due to hydrodynamic forces. The density of fluid is taken as ρ f = 1.0 g/cm 3 and the density of solid is ρ s = 1.1 g/cm 3 in this simulation. The viscosity of the fluid is ν f = 0.01 cm 2 /s. Initial angular velocity and incident angle of the airfoil are set to zero. The boundary conditions of flow are given as u = (0, 0) at y=-2 cm, 2 cm and u = (1, 0) cm/s at x = -4 cm, 16 cm. Those boundary conditions are used in Glowinski et al. [START_REF] Glowinski | A fictitious domain approach to the direct numerical simulation of incompressible viscous flow past moving rigid bodies: application to particulate flow[END_REF] and are adopted here in order to compare the results of two methods. The Reynolds number is around 101 based on airfoil length and the maximum inflow speed.

ω s u = 1 v = 0 u = 1 v = 0 u = 0, v = 0 u = 0, v = 0
This flow is quite challenging as the leading edge of the airfoil has very small radius of curvature. To resolve the flow near the leading edge, a good resolution of the Cartesian mesh is required.

Two sets of grids are chosen here to test the grid sensitivity, namely h = 1/96 cm and h = 1/64 and 0.25 respectively. The flow fields are shown in Figure 12. The airfoil keeps a stable position with its broadside perpendicular to the in-flow direction in the beginning and finally reaches a periodic motion of oscillation. The time histories of the rotational angle and the angular velocity are plotted in Figure 13. The results are pretty smooth. For comparison, we also include the results of Glowinski et al. [START_REF] Glowinski | A fictitious domain approach to the direct numerical simulation of incompressible viscous flow past moving rigid bodies: application to particulate flow[END_REF] in Figure 13 obtained by the distributed Lagrange multiplier method. Presents results match well those of Glowinski et al. [START_REF] Glowinski | A fictitious domain approach to the direct numerical simulation of incompressible viscous flow past moving rigid bodies: application to particulate flow[END_REF] . We observe that the direction of oscillation depends on the numerical errors. The oscillating direction changes over different calculation parameters. An opposite oscillating case can be found in Wan and Turek [START_REF] Wan | Direct numerical simulation of particulate flow via multigrid fem technique and the fictitious boundary method[END_REF] . Present direction of oscillation is in accord with Glowinski et al. [START_REF] Glowinski | A fictitious domain approach to the direct numerical simulation of incompressible viscous flow past moving rigid bodies: application to particulate flow[END_REF] . respectively. The results of Glowinski et al. [START_REF] Glowinski | A fictitious domain approach to the direct numerical simulation of incompressible viscous flow past moving rigid bodies: application to particulate flow[END_REF] are marked with circles.

A revised algorithm for the fluid-solid interaction was presented in this paper based on the moving immersed boundary method. Operator splitting is successively employed to achieve computational efficiency and modularity. First the pressure-velocity coupling is handled by a second order fractional step method, breaking the Navier-Stokes equations into a Helmholtz-type equation for the velocity and a Poisson equation for the pressure; Followed by the immersed boundary force splitting, a moving force equation is derived to account for the solid effect; At last a partitioned approach separates the fluid field from the solid field. By doing this, all the variables are decoupled and computed effectively with their own favored solvers. To achieve a more stable FSI coupling, a strongly coupled scheme is employed through a fixed point iteration. For better performance, we have removed the time-consuming pressure Poisson solver out from the FSI iteration and the computational time has been considerably reduced. Various cases have been tested, ranging from circular to non-circular objects with large displacements and rotations. The numerical results are compared to the benchmark, which demonstrates a good accuracy of the proposed method.

Figure 1

 1 Figure 1 illustrates a general description of a rigid body motion in an incompressible fluid,

FIG. 2 :

 2 FIG.2: Illustration of interpolation and spreading procedures with a discrete delta function.

3 5 (; 6 (s )ṽ n+1 f ; 7 ( 8 (, 9 if 16 (

 5678916 Fluid) Predict the velocity ṽn+1 f using(10); Initialize values: (•) k=0,n+1 = (•) n , where (•) includes x s , ẋs , F; 4 for k = 0 to k max do Fluid) Construct or update the interpolation operator matrix T (x k,n+1 s ) and the moving force coefficient matrix M (x k,n+1 s )Fluid) Interpolate the fluid velocity T (x k,n+1 Interface) Solve the moving force equation (12a) for F k+1,n+1 with ẋk,n+1 s ; Solid) Compute the solid equations for (x k+1,n+1 s ||(•) k+1,n+1 -(•) k,n+1 ||/||(•) k+1,n+1 || < tolerance then 10 (•) n+1 = (•) k+1,n+1 ; Fluid) Correct the fluid velocity to vn+1 f with F n+1 using (12b); 17 (Fluid) Solve the pressure Poisson equation and compute the final velocity v n+1 f

  at different times t = 0.2 s, 0.4 s, 0.6 s, 0.8 s. The cylinder quickly reaches a uniform falling velocity until it hits the bottom of the channel. We plot the flow quantities as a function of time in Figure 4, including the longitude position y c of the cylinder center, the vertical velocity v c , the Reynolds number Re and the translational kinetic energy E T . Here Re and E T are defined as Re = (ρ s D u 2 c + v 2 c )/µ and E T = 0.5m s (u 2 c + v 2 c ) respectively, where u c is the horizontal velocity component. For comparison, the results of Ref. 20 are included in Figure 4, taken from h = 1/96 cm. Good agreements have been obtained.

FIG. 3 :

 3 FIG. 3: Vorticity fields at different times t = 0.2 s, 0.4 s, 0.6 s, 0.8 s for the freely falling cylinder in a confined channel problem. The contour levels are set from -15 (blue) to 15 (red) with an increment of 1.

FIG. 4 :

 4 FIG. 4: Time evolution of longitude position y c , velocity v c of the cylinder center, the Reynolds number Re and the translational kinetic energy E T for the freely falling cylinder in a confined channel problem. "•", result of Ref. 20; "-", present result with h = 1/96 cm; "----", present result with h = 1/48 cm.

FIG. 6 :

 6 FIG. 5: Snapshots of vorticity fields for a freely falling cylinder in an open domain . The contour level is set from -6 (blue) to 6 (red) with an increment of 0.4.

FIG. 7 :

 7 FIG. 7: Comparison of different strong FSI coupling strategies. (a) The number of sub-iterations in the strong FSI coupling; (b) Time consumption ratio of various strong coupling methods compared to the weak coupling method. "•", strong coupling with pressure Poisson equation in the sub-iteration; "+", the proposed strong coupling.

FIG. 8 :

 8 FIG. 8: Computational domain of the elliptical particle sedimentation problem.

FIG. 9 :FIG. 10 :

 910 FIG. 9: Vorticity fields at different times: (a) t = 0.1 s ; (b) t = 0.3 s; (c) t = 0.5 s; (d) t = 1.0 s; (e) t = 1.5 s; (f) t = 2.0 s. The contour levels are set from -15 to 15.

FIG. 11 :

 11 FIG. 11: Computational domain of the flow past a rotating NACA0012 airfoil. The dashed lines represent the initial position of the airfoil.

FIG. 12 :

 12 FIG. 12: Instantaneous vorticity (a) and velocity (b) of the flow over a rotating NACA0012 airfoil.

FIG. 13 :

 13 FIG. 13: Time histories of the angle (a) and the angular velocity (b) of the rotating NACA0012 airfoil. h = 1/96 cm, solid line; h = 1/64 cm, dashed line. The measures are in rad and rad/s

TABLE I :

 I The drag, lift coefficients and the Strouhal number for the freely falling and rising circular cylinder in an open domain.

ACKNOWLEDGMENTS

This work was part of the PhD study of the first author, where the financial support of the China Scholarship Council is acknowledged. The computation was performed on the platform PILCAM2 at the Université de Technologie de Compiègne and the HPC at the Université de Strasbourg.


Xiaochen Huang

Boxin Huang

Yong He

Liang Feng

Jian Shi

Li Wang

Juan Peng
email: juan.wang@ens.psl.eu

Yong Chen
email: yong.chen@ens.psl.eu

Sars-Cov-2 spike protein-induced damage of hiPSC-derived cardiomyocytes

Keywords: Sars-Cov-2, Spike protein, hiPSC, Cardiomyocytes, ACE2, Ang II

Sars-Cov-2 may trigger molecular and functional alterations of cardiomyocytes of the heart due to the presence of receptor angiotensin-converting enzyme 2 (ACE2) of the host cells. While the endocytic itinerary of the virus via cleavage of the spike protein of Sars-Cov-2 is well understood, the role of the remaining part of the spike protein subunit and ACE2 complex is still elusive. Herein, we investigate the possible effects of this complex by using synthetic spike proteins of Sars-Cov-2, human induced pluripotent stem cells (hiPSC), and a culture device made of an arrayed monolayer of crosslinked nanofibers. hiPSCs were firstly differentiated into cardiomyocytes (CM) which formed cardiac tissue-like constructs with regular beating and expression of both ACE2 and gap junction protein Connexin 43. When incubated with the spike proteins, the hiPSC-CMs underwent a rhythmic fluctuation with overstretched sarcomere structures and dispersed gap junction proteins. When incubated with the spike proteins and supplementary Angiotensin II, the damage of the spike protein on hiPSC-CMs was enhanced due to downregulated ACE2, chromatin margination, altered Connexin 43 expression, sarcomere disruption, and beating break. This discovery may imply latent effects of the spike proteins on the heart.

Introduction

The symptoms of COVID-19 could be associated with injury and dysfunction of multi-organs, including lung, heart, blood vessels, kidney, gut, and brain [START_REF] Wiersinga | Pathophysiology, Transmission, Diagnosis, and Treatment of Coronavirus Disease 2019 (COVID-19): A Review[END_REF][START_REF] Maghool | Gastrointestinal and renal complications in Sars-Cov-2-infected patients: Role of immune system[END_REF][3][START_REF] Sato | Cardiovascular disease in Sars-Cov-2 infection[END_REF][START_REF] Bearse | Factors associated with myocardial Sars-Cov-2 infection, myocarditis, and cardiac inflammation in patients with COVID-19[END_REF][START_REF] Dekosky | Blood Biomarkers for Detection of Brain Injury in COVID-19 Patients[END_REF] for which all studies attribute the infection to the binding of the spike protein of Sars-Cov-2 and the receptor Angiotensin-converting enzyme 2 (ACE2) of the host cells [START_REF] Knyazev | Endocytosis and Transcytosis of Sars-Cov-2 Across the Intestinal Epithelium and Other Tissue Barriers[END_REF] . While the endocytic itinerary of the virus after cleavage of the spike protein is well understood, the role of the remaining part of the spike protein subunit 1 and ACE2 complex (spike-ACE2 complex) is still elusive. Indeed, in vitro models have shown systematically an overactive immune response of the cells to Sars-Cov-2 [START_REF] Chen | Studying Sars-Cov-2 infectivity and therapeutic responses with complex organoids[END_REF] but less has been devoted to the effect of the spike-ACE2 complex.

The endocytosis of the virus is initiated with a two-step cleavage on the cell membrane with help of Furin and TMPRSS2 surface proteins [START_REF] Zhang | Molecular mechanism of interaction between Sars-Cov-2 and host cells and interventional therapy[END_REF] . After a membrane fusion, the virus is endocytosed and replicated in the cytoplasm, while the spike-ACE2 complex remains on the cell membrane (Figure 1A). Intuitively, the spike-ACE2 complexes could have three possible ends: (1) Endocytosed into the cytoplasm; (2) Total or partial shed from cell membrane; (3) Stay on the cell membrane. Previous studies have shown that the spike protein could be endocytosed by clathrin mediation due to its ACE2-dependence [START_REF] Bayati | Sars-Cov-2 uses clathrinmediated endocytosis to gain access into cells[END_REF] . The spike protein-induced downregulation of ACE2, activation of several signaling pathways, and inflammation could also be observed [11] [12] . The objective of the present study is to investigate whether the spike protein alone, outside the context of infectious virus, can trigger molecular and functional alterations in cardiac complications which occur in 20% to 44% of hospitalized patients, and constitute an independent risk factor for COVID-19 mortality [START_REF] Sharma | Human iPSC-derived cardiomyocytes are susceptible to Sars-Cov-2 infection[END_REF][START_REF] Navaratnarajah Chanakha | Highly Efficient Sars-Cov-2 Infection of Human Cardiomyocytes: Spike Protein-Mediated Cell Fusion and Its Inhibition[END_REF][START_REF] Perez-Bermejo Juan | Sars-Cov-2 infection of human iPSCderived cardiac cells reflects cytopathic features in hearts of patients with COVID-19[END_REF][START_REF] Bailey | Sars-Cov-2 Infects Human Engineered Heart Tissues and Models COVID-19 Myocarditis[END_REF][START_REF] Xie | Long-term cardiovascular outcomes of COVID-19[END_REF] . Human-induced pluripotent stem cells (hiPSC) are used for differentiation toward hiPSC-CMs with an arrayed nanofiber membrane (culture patch) [START_REF] Tang | Induction and differentiation of human induced pluripotent stem cells into functional cardiomyocytes on a compartmented monolayer of gelatin nanofibers[END_REF] . The derived hiPSC-CMs (cardiac patch) are then incubated with synthetic spike proteins in the presence or absence of supplementary Angiotensin II (Ang II).

The hiPSC-CMs are chosen due to their modeling relevance and easy observation of the beating behavior which should be sensitive to Sars-Cov-2. Indeed, the sarcomere structures of hiPSC-CMs which are tightly correlated to the beating can be quantitatively analyzed by fluorescence imaging of α-actinin and troponin T marked cells. Moreover, the morphology and the beating performance of hiPSC-CMs are determined by the gap junctions which can also be analyzed by fluorescence imaging of actin and connexin 43 marked cells. Finally, both ACE2 and Ang II are involved in the regulation of the renin-angiotensin system (RAS) since ACE2 can hydrolyze Ang II to yield Angiotensin 1-7 (Ang 1-7) which has opposite effects of Ang II. While Ang 1-7 promotes vasodilation and anti-inflammation, Ang II may cause vasoconstriction and inflammation [19]. Therefore, the presence of Ang II is meaningful for the investigation of the spike effect and the downregulation of ACE2 is a risk factor of RAS abnormally [20]. Accordingly, the changes of the cardiac beating, sarcomere, and nuclear structure, and the expression of connexin 43 and ACE2 are analyzed to illustrate the structural and functional damage induced by the spike proteins.

Results

Beating and structure of hiPSC-CMs

The culture devices (patch) used in this work for hiPSC-CMs differentiation were produced by lithography, electrospinning, and chemical crosslinking techniques [START_REF] Tang | Induction and differentiation of human induced pluripotent stem cells into functional cardiomyocytes on a compartmented monolayer of gelatin nanofibers[END_REF] . The characteristics of the nanofiber patches are summarized in Figure S1, including the porosity, the average diameter of the nanofibers, permeability, and stiffness of the nanofiber membrane, together with SEM pictures to illustrate the morphology of the nanofiber membrane. Statistically, the size of 98% of pores was smaller than 7 µm, on which cells can stay without falling down. By measuring the flow rate of DI water across the membrane as a function of pressure, a membrane permeability of 3.44×10 3 nm² could be deduced, which is much larger than that of the commonly used etch ion tracking membrane [START_REF] Liang | A microfluidic system for rapid determination of Darcy permeability of filtration membranes: from patterned micro-holes to artificial basement membranes[END_REF] . By measuring the deflection as a function of pressure, an effective Young's modulus of the membrane in the order of 4 MPa could be obtained [22] . Accordingly, an effective in-plane Young's modulus in the order of 9.5 kPa could be estimated by using a nanofiber diameter of 0.4 µm, pore size of 6 µm, and curve fitting based on a theory of cellular solids [START_REF] Gibson | Cellular Solids: Structure and Properties[END_REF] . Such an effective Young's modulus is well-suited to support hiPSC-CM differentiation [START_REF] Wang | Substrate elasticity dependent colony formation and cardiac differentiation of human induced pluripotent stem cells[END_REF] . Finally, the cell adhesion strength was evaluated with a spinning technique, showing a better attachment of the cells on nanofibers than on glass substrates [START_REF] Huang | Generation of Interconnected Neural Clusters in Multiscale Scaffolds from Human-Induced Pluripotent Stem Cells[END_REF] .

To generate hiPSC-CMs, the standard differentiation protocol was used [START_REF] Tang | Induction and differentiation of human induced pluripotent stem cells into functional cardiomyocytes on a compartmented monolayer of gelatin nanofibers[END_REF] (Figure 2A). hiPSCs were firstly seeded on the nanofiber patch to form dome-like hiPSC colonies and then differentiated to cardiomyocytes in each of the membrane areas (Figure S2). During the first 7 days, the hiPSC colonies became flat due to cardiac mesoderm induction. Starting from Day 9, 1-5% of hexagonal areas started beating slowly. The percentage of beating areas increased to 20% on day 15 and to 80% on day 64 (Supplementary Movie 1). Meanwhile, the beating rate of the hiPSC-CMs increased from 65±13 bpm on day 15 to 153±4 bpm on day 64 (Figure 2B) by analyzing the timelapse images with a MATLAB program [START_REF] Maddah | A Non-invasive Platform for Functional Characterization of Stem-Cell-Derived Cardiomyocytes with Applications in Cardiotoxicity Testing[END_REF] . By dividing the time-lapse images into 30 areas (Figure 2C), three types of beating waveforms could be identified as plotted in Figure 2D. Statistically, 26 over 30 areas (marked in green) showed regular and high magnitude beatings (plotted in green). The two areas with low magnitude beatings (marked in yellow) fell mostly in the frame area, whereas the fluctuation in beating amplitude of the two red marked areas can be attributed to the inhomogeneity of the tissue construct. Nevertheless, the beating frequency of all areas was the same.

To characterize the fabricated hiPSC-CMs, α-smooth muscle actin (SMA), a marker for smooth muscle cells and myofibroblasts [START_REF] Shinde | The role of α-smooth muscle actin in fibroblast-mediated matrix contraction and remodeling[END_REF] , and troponin T (TnT), a structural marker of cardiomyocytes, were used for immune-stained analysis. In Figure 2E, immunofluorescence images showed the expression of SMA and troponin T on the hiPSC-CMs. Clearly, the derived hiPSC-CM tissue constructs consisted of at least three types of cells (TnT+, TnT-/SMA+, and TnT-/SMA-), in consistent with the previous observations [28] [29] . Moreover, the cross-section views of the tissue construct in Figure 2E showed that the cells in the hexagonal area were stacked to form a multilayer structure where the cells on the frame or nearby are mostly spindle-shaped TnT-SMA+ cells. Interestingly, the TnT-SMA+ cells formed a thin layer at the perimeter (top layer), TnT+ cells were localized below the TnT-SMA+ cell layer (middle layers), and nostained cells were also found at the bottom. Schematically, the hiPSC-CM tissue construct can be divided into three layers due to tissue-like self-organization (Figure 2F). Consequently, the percentage of each cell type in three layers could be estimated (Figure 2G). Apparently, the TnT+ cells were quite homogeneously distributed in all layers, while a high percentage of TnT-SMA+ cells (17.5±4.58%) were in the top layer compared to that in the middle (2.52±1.82%) and the bottom (4.92±3.09%) layers. In contrast, more TnT-SMA-cells (46.72±3.13%) were in the bottom layer compared to that in the top (12.5±4.05%) and the middle (22.01±2.54%) layers.

To evaluate the suitability of the hiPSC-CMs to the studies of the effects of spike protein, the expression levels of ACE2 and connexin 43 were analyzed. As illustrated in Figure 2H, 74.84±3.05% of the cells expressed connexin43, and 76.55±3.88% of the cells expressed ACE2, some of them around the nucleus in agreement with the previous studies [START_REF] Ahmad Mulyadi Lai | Expression of Endogenous Angiotensin-Converting Enzyme 2 in Human Induced Pluripotent Stem Cell-Derived Retinal Organoids[END_REF] . The sarcomere structure of the hiPSC-CMs can be seen from α-actinin and troponin T immunostained images (Figure 2I). Here, α-actinin is specific for the attachment of actin filaments to the Z-lines of the hiPSC-CMs, and troponin T is associated with the coupling between tropomyosin and actin filament. The observed period arrangement of both α-actinin and troponin T were characteristic of cardiac sarcomeres which confirmed the maturation of our hiPSC-CMs. Statically, 80.29±5.33% of the troponin T were aligned in the troponin T expressed areas. Likewise, 66.11±2.29% of the α-actinin formed stripy patterns in the α-actinin expressed areas.

Functional and structural changes of hiPSC-CMs caused by spike proteins

Effect on the rhythmic behavior

The effect of the spike protein on the rhythmic activity of cardiomyocytes was studied by co-culturing the spike protein with the hiPSC-CMs. A high dosage of synthetic biotinylated spike protein (1 µg/mL) was added to the culture medium for incubation 24 h. To track the trajectory of the spike proteins, the cells were fixed and stained with streptavidin-594 for spike proteins and phalloidin-FITC for actin. The results showed striated actin structures, confirming the cardiac cell type (indicated with white arrows in Figure 3A). Spike proteins were found inside 88.1±4.56% of the cells and they were mostly localized around the nucleus of hiPSC-CMs (Figure 3A). This means that the spike proteins partially diffused and that the striated organization of actin filaments could be affected, in agreement with the previous studies [START_REF] Schneider | Sars-Cov-2 direct cardiac damage through spike-mediated cardiomyocyte fusion[END_REF] .

After introducing the spike proteins, the beating frequency of the hiPSC-CMs varied. Figure 3B shows the change of normalized beating frequency over 24 h.

Compared to the samples without spike protein treatment, the beating frequency of the samples with spike proteins increased during the first 3 h, then decreased until 18 h.

During the last 6 h, the beating frequency of the samples with and without treatment was comparable. This result suggested that the spike protein has a limited impact on the beating frequency of hiPSC-CMs and that apparently the cells remained alive (Supplementary Movies 2, 3). We also noticed a remarkable change of the beating waveform after 18 h incubation and then a recovery at 24 h, which can be attributed to the spike protein-induced perturbation and the cellular regulation.

Disturbed intracellular structure

The spike protein-induced beating fluctuation of hiPSC-CMs may be correlated to the disturbance of sarcomere structures. The expression and distribution of α-actinin and troponin T of the samples with and without spike protein treatment were then compared (Figure 3C). Without spike protein, striated organizations of α-actinin and troponin T were observed in the immunostaining images (periodic narrow stripes in red and double stripes in green, respectively). Here, the structure of troponin T marked double stripes which represent the relaxation phase of sarcomeres was intentionally selected for comparison. With spike protein, the α-actinin associated stripes cannot be seen and the distance between the troponin T marked double stripes remarkably increased. In fact, the α-actinin proteins aggregated or diffused irregularly (Figure 3C, 3F), keeping the averaged intensity of α-actinin unchanged (Figure 3E). Normally, the distance between the troponin T double stripes varies with the contraction phase of the sarcomere, which is maximum when the sarcomere is completely relaxed and minimum when it is contracted. The sarcomere length (SL) of hiPSC-CMs with and without spike protein treatment could then be quantified (Figure 3D). Statistically, it varied between 0.70 µm to 2.17 µm for the cells without the treatment. After the spike protein treatment, 25.56±7.18% of the sarcomeres have an SL showed a maximum value of 2.68 µm (Figure 3G). Since the sarcomeres are organized in the form of cross and slidable bridges of actin and myosin filaments, the excess in SL may cause of weakening of contraction ability [START_REF] Rassier | Sarcomere mechanics in striated muscles: from molecules to sarcomeres to cells[END_REF] . Together, the results suggest that the spike protein impacted indeed the sarcomere structures by disrupting the organization of α-actinin and increasing the distance between the double stripes of troponin T marked structure.

Remodeled gap junction

Gap junctions are composed of connexins in the form of membrane channels, which facilitate communication between neighboring cells and thus coordinate rge electrical and contractile activity of the heart [START_REF] Kurtenbach | Gap junction modulation and its implications for heart function[END_REF] . The spike protein-induced rhythmic change might be associated with the remodeling of the gap junctions. Therefore, the samples were immunostained with connexin 43 and actin for outlining the cell periphery. Without spike protein treatment, connexin 43 appeared spotty pattern on the intercellular space between cells. After the spike protein treatment, connexin 43 was aggregated around the nucleus (Figure 3H and Figure 3I). Statistically, 74.19±3.45% of the cells expressed connexin 43 in the perinuclear areas after the spike protein treatment (Figure 3J), suggesting an induced remodeling of the gap junction due to alternation of the connexin 43 expression caused by endocytosis [START_REF] Totland | Regulation of gap junction intercellular communication by connexin ubiquitination: physiological and pathophysiological implications[END_REF] .

Damages caused by spike protein in the presence of supplementary Ang II

Beating arrest of hiPSC-CMs

The changes of hiPSC-CMs caused by the spike protein may be related to the alternation of the available ACE2. By coculture of hiPSC-CMs and spike proteins with supplementary Ang II, we observed a significant change of cell vitality. To have a striking effect, a high dosage of Ang II (100 ng/mL) was added to the culture medium after incubation for 3h with or without the spike proteins. After incubation for 48 h, the medium was refreshed without spike protein and Ang II.

Figure 4A shows a comparison of the normalized beating frequency as a function of time over 54 h of three samples with and without the spike protein and Ang II treatment. Apparently, the hiPSC-CMs without the spike protein treatment remain alive, even a relatively high dose of Ang II was used. The beating frequency of the hiPSC-CMs treated with Ang II increased in the first 6 h and then decreased smoothly in the following 18 h. Afterward, the beating frequency of the samples with and without Ang II treatment varied in a similar manner, suggesting that the damaging effect of Ang II on the beating activity of the hiPSC-CMs is limited and that the cells remain alive after moving the Ang II. In contrast, the spike protein treatment with supplementary Ang II caused a fatal change in the rhythmic activity of the hiPSC-CMs. As can be seen, the variation of the beating frequency during the first 29 h is comparable to that of the sample treated with Ang II. Then, the beating frequency gradually decreased and the cardiac beating completely stopped in the end even after renewing the medium (Supplementary Movies 4,5, 6). This can be considered as a direct evidence of the damaging effect of the spike proteins on hiPSC-CMs in the presence of supplementary Ang II.

Altered nuclear morphology

The beating arrest of hiPSC-CMs caused by spike protein with supplementary Ang II may partially be due to disturbance of the ACE2 expression. Figure 4B shows immunofluorescent staining images of ACE2 under different culture conditions, indicating a variation in ACE2 expression. In the presence of spike protein and Ang II, the intensity of ACE2 marked fluorescence was significantly lower than that without the spike protein (Figure 4C). during both necrosis and apoptosis, while chromatin margination is a typical symbol of apoptosis [START_REF] Cummings | Measurement of cell death in mammalian cells[END_REF] . This result suggested spike protein may cause cell death through apoptosis consistent with the finding reported before [START_REF] Sharma | Human iPSC-derived cardiomyocytes are susceptible to Sars-Cov-2 infection[END_REF] .

Disrupted sarcomere and impaired expression of connexin 43

The spike protein with Ang II-stopped beating of hiPSC-CMs may be related to the disruptions of sarcomere structures and intracellular connection. Figure 5A indicates that the immunostaining of troponin T and α-actinin was varied depending on the treatment. Compared to the samples without treatment, consistently these proteins kept assembling regularly as striated structures in most cells incubated with Ang II.

While with the treatment of spike protein and Ang II, the specific organization of troponin T and α-actinin to represent sarcomere was changed to random distribution in most of the cells. Figure 5B shows the statistical plots of the mean intensity of troponin T and α-actinin fluorescence of the samples. With the treatment of spike protein in the presence of Ang II, the expression level of both troponin T and α-actinin was lower than that of the two other samples. Likewise, Figure 5C shows the fluorescent area of troponin T and α-actinin of the sample after spike protein and Ang II treatment, indicating a remarkable variation compared with the two other samples. In other words, the regular striated distribution of the troponin T and α-actinin became a sparse random distribution. This suggests that the supplementary Ang II aggravates the effect of the spike protein on sarcomere structure from overstretched troponin T bands to randomly distributed ones. Accordingly, the expression level of connexin 43 in the hiPSC-CMs incubated with spike protein and Ang II was significantly lower than that without spike protein due to the altered actin cytoskeleton and the absence of sarcomere structures (Figure 5D,5E).

Discussions

COVID-19 pandemic has already lasted for more than two years with the emergence of several Sars-Cov-2 variants due to spike protein mutations. Expanded research is necessary to develop new models and to improve our understanding of Sars-Cov-2 infection mechanisms. In this regard, in-vitro modeling of the virus-host interaction by using hiPSC-derived tissue-like constructs could be significant. Amongst other tissue types, hiPSC-CM based modeling is specific due not only to its availability [36-40] but also to the simplicity of its damage observation, as shown in this work.

Previously, we have shown that the response of hiPSC-CMs to external stimuli is substrate stiffness dependent and that hiPSC-CMs on a substrate of stiffness of ~9 kPa are more sensitive to electric stimulation [24]. The present study was based on a culture device with an effective in-plane Young's modulus of the same order. We also found a tissue-like cell organization of hiPSC-CMs with 77.47% of CMs (Troponin T+ cells), 15.77% of myofibroblasts (SMA+ cells), and 6.76% of troponin T-and SMA-cells on our nanofiber patch. While the majority of myofibroblasts were localized on the top of the tissue-like structure, the troponin T-and SMA-cells were mostly on the bottom of the structure. These cells are mostly stromal cells but some of them could be undifferentiated cells. Moreover, a few pace making cells may be enough to induce beating, with the faster cells over-riding quiescent cells. This discovery shows the tissue-like features in the presence of myofibroblasts could be improved in a similar way to that observed in hiPSC-derived cardiac spheroids [START_REF] Beauchamp | 3D Co-culture of hiPSC-Derived Cardiomyocytes With Cardiac Fibroblasts Improves Tissue-Like Features of Cardiac Spheroids[END_REF] . Compared to other engineering techniques [START_REF] Joanne | Nanofibrous clinical-grade collagen scaffolds seeded with human cardiomyocytes induces cardiac remodeling in dilated cardiomyopathy[END_REF][START_REF] Orlova | Electrospun nanofibers as a tool for architecture control in engineered cardiac tissue[END_REF][START_REF] Li | Human Pluripotent Stem Cell-Derived Cardiac Tissue-like Constructs for Repairing the Infarcted Myocardium[END_REF][START_REF] Beauchamp | 3D Co-culture of hiPSC-Derived Cardiomyocytes With Cardiac Fibroblasts Improves Tissue-Like Features of Cardiac Spheroids[END_REF][START_REF] Hoang | Engineering spatial-organized cardiac organoids for developmental toxicity testing[END_REF] , the suspended culture of hiPSC-CMs has additional advantages of easy stiffness regulation and a higher degree of mini tissue-like selforganization.

By using hiPSC-CMs with tissue-like features, we were able to demonstrate the structural and functional damage effects of synthetic Sars-Cov-2 spike proteins which are schematically summarized in Figure 6, including i) sarcomere structure, ii) connexin 43 distribution, and iii) chromatin phenotype. The functional consequence of these structural changes was the variation of the rhythmic activities, including desynchronization and beating frequency fluctuation without supplementary Ang II and beating stop with supplementary Ang II.

Synthesized biotinylated spike proteins of 1 µg/mL (37nM) concentration were used. The number of spike proteins per mL at this dosage is ~ 210 [START_REF] Sharma | Human iPSC-derived cardiomyocytes are susceptible to Sars-Cov-2 infection[END_REF] , which is significantly larger than the available ACE2 receptors for a total cell number of 10 6 .

Experimentally, this dosage was sufficiently high for the observation of the effect of spike protein-mediated cell signaling [START_REF] Suzuki | Sars-Cov-2 spike protein-mediated cell signaling in lung vascular cells[END_REF] . Nevertheless, cells in culture could still be alive when the spike protein was used alone. Covid-19 Clinically, the free spike protein could be detected from the serum samples of patients in the range of 2.5-17.5 µg/mL [43] , which is much higher than the concentration we used. Therefore, the testing concentration of the spike proteins would mimic the minimal effect of free spike protein on cells. On the other hand, the dosage of Ang II we used was about 100 times higher than that of the physiological and pathological conditions [START_REF] Osman | Expression of ACE2, Soluble ACE2, Angiotensin I, Angiotensin II and Angiotensin-(1-7) Is Modulated in COVID-19[END_REF] . However, our data suggested that the Ang II alone at this high dosage has a much less injury effect than that of the combined use with the spike proteins.

The spike protein interacts with ACE2 through its RBD site. As illustrated in Figure 1, the spike can enter into the cell, shedding or staying on the cell membrane.

Whether the entry of spike proteins is through endocytosis is unclear. However, the presence of the spike proteins in the intracellular areas was demonstrated in this work (Figure 2) and others [START_REF] Kim | Spike Proteins of Sars-Cov-2 Induce Pathological Changes in Molecular Delivery and Metabolic Function in the Brain Endothelial[END_REF] . It seems that the spike proteins were mostly located near the nuclei. Although the structural changes of the cells mentioned above could be attributed to the spike proteins, the mechanisms underlying remains to elucidate. Previous studies have discussed that the actin cytoskeleton rearrangement occurs due to coronavirus internalization [START_REF] Wen | Cytoskeleton-a crucial key in host cell for coronavirus infection[END_REF] and the spike protein could bind with nonmuscle myosin in pulmonary epithelial cells [START_REF] Chen | Nonmuscle myosin heavy chain IIA facilitates Sars-Cov-2 infection in human pulmonary cells[END_REF] . In the case of hiPSC-CMs, such an attachment may cause dysfunction of the sarcomere.

The overstretched sarcomere caused by spike proteins may decrease the beating frequency and dysfunction of actin/myosin filaments or their bridging structure [START_REF] Prondzynski | Disease modeling of a mutation in α-actinin 2 guides clinical therapy in hypertrophic cardiomyopathy[END_REF] . Furthermore, the altered localization of connexin 43 remodeled the intercellular connection, since connexin 43 is a component of gap junction localized on the cell membrane to play a crucial role in the synchronized contraction of the hiPSC-CMs. In the presence of spike proteins, connexin 43 was dispersed in the cytoplasm close to the nucleus. Consequently, the synchronized beating might be disrupted.

The decrease of ACE2 expression level due to spike protein or Sars-Cov-2 infection was previously observed [START_REF] Lei | Sars-Cov-2 Spike Protein Impairs Endothelial Function via Downregulation of ACE 2[END_REF] . The related effects on immune response, or cytokine production of the cell due to the imbalance of RAS system has also been discussed [START_REF] Aleksova | Effects of Sars-Cov-2 on Cardiovascular System: The Dual Role of Angiotensin-Converting Enzyme 2 (ACE2) as the Virus Receptor and Homeostasis Regulator-Review[END_REF] . A pathological amount of Ang II (100 ng/mL) was added, which could affect cardiac beating rate and Ca 2+ signaling pathway [START_REF] Sedan | 1, 4, 5-Inositol trisphosphate-operated intracellular Ca2+ stores and angiotensin-II/endothelin-1 signaling pathway are functional in human embryonic stem cell-derived cardiomyocytes[END_REF] even without spike protein.

This concentration has been reported to affect cardiac activities [START_REF] Lagerqvist | Endothelin-1 and angiotensin II modulate rate and contraction amplitude in a subpopulation of mouse embryonic stem cell-derived cardiomyocyte-containing bodies[END_REF] but is not fatal to cardiomyocytes [START_REF] Gao | Angiotensin II induces apoptosis in human induced pluripotent stem cell-derived cardiomyocytes[END_REF] . Indeed, more cells underwent apoptosis and necrosis and more cells showed abnormal nuclear morphology. After a short period of incubation, however, most of the cells were alive. The change of the cell state can be attributed to the unbalanced signalization of pro-and anti-inflammation pathways due to excess of Ang II. Indeed, Ang II can be converted by ACE2 to Ang 1-7 to activate Mas signaling pathways while unconverted Ang II mainly activates AT1R. In the case of a high dosage of Ang II, excess of Ang II occurs since the conversion rate of ACE2 is limited. The signaling pathways mediated by AT1R can then be over-activated, giving rise to inflammation-like behaviors of the cells. In the presence of both spike proteins and Ang II, ACE2 was partially consumed so that the conversion from Ang II to Ang 1-7 becomes less efficient so that the signaling pathways mediated by AT1R can be even more over-activated, giving rise to more pronounced damage effects.

The presence of the spike protein in the intracellular area was previously observed when cells were infected by the coronavirus. Then, cytokine production, sarcomere disassembly, and cell death could be expected [START_REF] Bailey | Sars-Cov-2 Infects Human Engineered Heart Tissues and Models COVID-19 Myocarditis[END_REF] . Our data provide new insights on the damage of hiPSC-CMs caused by the spike-ACE2 complex, independent from the mechanism of virus replication. Recent studies have shown that the full-length spike protein may act as a ligand to induce non-infective cellular stress by activation of the extracellular signal-regulated kinase 1/2 (ERK1/2) through the CD147 receptor in human cardiac pericytes [START_REF] Avolio | The Sars-Cov-2 Spike protein disrupts human cardiac pericytes function through CD147 receptor-mediated signalling: a potential non-infective mechanism of COVID-19 microvascular disease[END_REF] . Together, the effect of the spike protein might be significant, suggesting that adjuvant care might be necessary. More systematic investigations are expected by considering different dosages of spike protein and different concentrations of Ang II as well as antagonist drugs such as ACE2 inhibitors and AT1R-blockers. Finally, the effects of the spike-ACE2 complex released after endocytosis of Sars-Cov-2 under pathological conditions must be studied to elucidate the joint effect of the virus replication and the complex.

Conclusions

We reported a study of the effect of Sars-Cov-2 spike protein on hiPSC-CMs. A cardiac tissue-like construct was generated by the differentiation of hiPSCs into CMs on an arrayed monolayer of crosslinked nanofibers. Our results showed that the spike

Experimental Section

Patch fabrication

The patch was fabricated as described previously [START_REF] Tang | Induction and differentiation of human induced pluripotent stem cells into functional cardiomyocytes on a compartmented monolayer of gelatin nanofibers[END_REF] . Firstly, a hexagonal microframe with an outer ring was patterned by using two-level photolithography and vacuum-assisted UV molding techniques. The pattern was defined in a chromium mask with a micropattern generator (μPG 101, Heidelberg Instruments, Germany) and chemical etching. A photoresist layer (AZ 40XT, MicroChem, France) of 50 µm thickness was spin-coated on the mask and exposed from the backside with UV light. Afterward, a ring pattern was added by spin-coating another 50 µm photoresist layer and frontside UV-exposing. After development, the resist pattern (mold) was antisticking treated in trimethylchlorosilane (TMCS) vapor for 1 min. Then, a solution of polydimethylsiloxane (PDMS) was prepared by mixing PDMS pre-polymer and its cross-linker (GE RTV 615, GE, France) at a weight ratio of 10:1 and poured onto the mold. After bubble removal, PDMS was cured at 75 °C for 4 h and peeled off from the mold. Next, the PDMS replica was placed on a glass slide (pattern side in touch of the glass) for vacuum-assisted UV molding. After gas removal (5 min), a low viscous photoresist, OrmoStamp (Micro Resist Technology), was introduced to fill the cavity between PDMS and glass. After UV exposure for 120 s at 21.6 mW cm -2 , the solidified OrmoStamp structure (frame) was separated from the PDMS mold and the glass slide.

A monolayer of gelatin nanofibers was deposited on the frame by electrospinning and chemical crosslinking. To collect the nanofibers more uniformly, the frame was coated with Au with a sputter (K6785X, Emitech) and placed on an aluminum-covering dish (collector). Gelatin solution was prepared with 12 wt% gelatin powder (G2625, Sigma) and a solvent mixture of DI water, ethyl acetate, and acetic acid at a volume ratio of 10: 14: 21. The gelatin solution was ejected at a flowrate of 0.2 mL/h at a bias voltage of 11 kV and a distance of 10 cm between the stainless-steel needle (23-gauge) of the syringe and the collector. After electrospinning, the sample was dried overnight in a desiccator to remove the residual solvent and the nanofibers were cross-linked in an ethanol solution containing 0.2 M N-hydroxy succinimide (NHS, Sigma) and 0.2 M 1-ethyl-3-(3-dimethylaminopropyl) carbodiimide hydrochloride (EDC, Sigma) for 4 h.

Finally, the sample (patch) was washed three times in ethanol and dried in a vacuum overnight to eliminate the residual solvent.

SEM characterization

The patch was soaked in 30% ethanol solution (in DI water) for 30 min, and dehydrated by subsequent immersion in 50%, 70%, 80%, 90%, 95%, and 100% ethanol in DI water, each for 10 min and dried with a vacuum pump. Afterward, a gold layer of 2 nm thickness was deposited with a sputter and the samples were observed by using a scanning electron microscope (SEM, Hitachi S-800, 10 kV).

Differentiation of hiPSC-CMs

hiPSCs (Human episomal line, Thermofisher Scientific) were cultured with complete E8 flex medium in a dish coated with 1:100 diluted vitronectin (A14700, Thermo Fisher) and placed in an incubator of 37 °C in 5% CO2. The medium was refreshed every two days until the cells reached ~80% confluence. Cells were then collected by dissociation with a 0.5 mM EDTA DPBS solution. To promote the adhesion of hiPSCs on gelatin nanofibers, the patch was coated with 1:500 diluted vitronectin in DPBS and then fixed in a culture dish with a PDMS hollow cylinder. 2 × 10 5 hiPSCs in 500 µL E8 flex medium containing 10 µM Y-27632 (ROCK inhibitor, Y0503, Sigma) were dropped on the culture patch through the hollow area of the cylinder and incubated for 2 h. Afterward, the hollow cylinder was removed and a fresh E8 flex medium of 2 mL without 10 µM Y-27632 was gently added to the dish. After 48 h, hemisphere hiPSC colonies were formed for cardiac differentiation using the protocol of Lian et al [START_REF] Shinde | The role of α-smooth muscle actin in fibroblast-mediated matrix contraction and remodeling[END_REF] .

Briefly, RPMI 1640 medium containing B27 minus insulin supplement was used as a base medium during the cardiac differentiation. On day 0, the E8 flex culture medium was replaced by the basal medium with 12 µM of CHIR99021 (SML1046, Sigma). It is a GSK3 inhibitor to induce mesoderm differentiation. After 24 h incubation, CHIR99021 was removed. Then cells were cultured with the basal medium for 48 h. On day 3, 5 µM of IWP2 (I0536, Sigma), the WNT activator was added into the basal medium for cardiac mesoderm induction. On day 5, the medium was replaced with the basal medium. Then, the culture medium was replaced by RPMI 1640 containing B27 every two days. Generally, cells started beating around days 8 and 10 and synchronous contraction was observed between days 21 and 30.

Incubation with spike protein of Sars-Cov-2 or Ang II

The recombinant biotinylated Sars-Cov-2 spike protein (bs-46003P-Biotin) and Ang II (A9525) were purchased from Bioss and Sigma. The hiPSC-CM samples of at least 60 days of differentiation were used for spike protein or Ang II treatment.

For the experiment for spike protein treating, 1 µg/mL spike protein was added to cardiac patch for 24h, another group was treated with the normal medium as control.

For the experiment of spike protein and Ang II treatment, the cardiac patch was treated with 1µg/mL spike protein and 100 ng/mL Ang II under three different conditions: (1) Spike + Ang II: Spike protein treatment for 48h, together with Ang II added after 3h of spike protein treatment; (2) Ang II only: Ang II treatment for 48h; (3) No treatment: normal medium culture for 48h.

Time-lapse imaging and beating analyses

The beating videos of the cardiac patch under different conditions were captured as sequential phase-contrast images by a live cell imaging system-Muvicyte (Perkin Elmer, France) with a 10× objective at a frame rate of 20 fps, a resolution of 1396 × 1052 pixels. At least 3 positions on the cardiac patch were selected for further monitoring. Before any treatment, a 20s video was captured per position as the baseline.

After treatment, 20s video was captured per 3h for 24h or 54h depending on the experiment. All captured videos were analyzed by a MATLAB program to acquire beating frequency and beating waveform.

To evaluate the beating frequency as a function of incubation time, the recorded time-lapse images were spatially segmented into 30 rectangle areas. The beating frequency of each area was determined and averaged over the whole image. Data obtained with an untreated hiPSC-CM sample were used as a baseline to normalize the beating frequency of treated samples for each time point. To ensure a statistical relevance, data were collected from 3 or 5 randomly selected positions of each sample and then averaged before reporting.

Immunostaining characterization

The cardiac patches were rinsed with PBS to remove detached and dead cells. Then, the patches with cells were fixed in a solution of 4% paraformaldehyde for 15 min and permeabilized in 0.5% Triton X-100 for 10min at room temperature. For non-specific binding blocking, cells were treated with a PBS solution containing 0.1% Tween-20, 0.1% Triton X-100, and 3% bovine serum albumin for 2h at room temperature.

Afterward, the samples were treated with primary antibodies, i.e., anti-α-actinin were used to stain actin when necessary. Finally, the stained samples can be observed using an LSM Zeiss 710 confocal microscope. For the fluorescent intensity analysis in each experiment, all data were analyzed by Image J software. For the nuclei morphology analysis, the nuclei of 120-150 cells from 10 images per sample were counted and analyzed. For the sarcomere length analysis, 252 sarcomeres from 12 images per sample were analyzed.

Statistic analysis

Data are presented as mean ± standard error. The statistical differences between groups were performed using Student's t-test. Differences were considered statistically significant at P value<0.05 ( * P < 0.05; * * P < 0.01; * * * P < 0.001). The P values are indicated in the figure captions where necessary. The Sars-Cov-2 spike protein can affect the behavior of host cells. By using cardiomyocytes derived from human induced pluripotent stem cells, this study shows the clear effects of the spike protein on cardiac beating, sarcomere structure, and distribution of gap junction proteins. In the presence of supplementary angiotensin II, the effects become more evident, showing even a beating arrest of the cardiomyocytes.

  The effect of the spike protein on hiPSC-CMs in the presence of Ang II was further analyzed by observing the nuclear morphology. Typically, the DAPI staining images show three types of nuclear morphology as illustrated in Figure4D: (1) chromatin diffusion homogenously in the nucleus; (2) the nucleus is split into multiple small pieces (nuclear fragmentation); (3) chromatin condenses as clusters at the edge of the nuclear membrane (chromatin margination). Compared to the samples without treatment, the presence of Ang II or spike protein with Ang II increased the percentage of nuclei having nuclear fragmentation comparably (Ang II/No treatment: ~2.83 times; Spike protein + Ang II/No treatment: ~3.43 times). While unlike the samples with no treatment or treated with Ang II, half of the nuclei (53.70±4.0%) treated with spike protein and Ang II appeared chromatin margination. Nuclear fragmentation appears

  proteins alone were able to provoke a damaging effect in terms of rhythmic fluctuation overstretching of sarcomere structures, and dispersion of gap junction proteins. With supplementary Ang II, the structure and function of the hiPSC-CMs underwent drastic changes due to downregulated ACE2, chromatin margination, Connexin 43 mislocalization, sarcomere disruption, and inevitably beating break. Further studies on the effect of the spike proteins shed from Sars-Cov-2 should enrich our understanding of the COVID-19 pandemic.

(

  Sarcomeric) antibody (1:200 diluted, Sigma, Japan), anti-Troponin T-C (1:200 diluted, Santa Cruz), anti-connexin 43 (1:100 diluted, Sigma), and anti-ACE2 (1:100 diluted, Santa Cruz) in blocking buffer overnight at 4 °C. Afterward, the samples were incubated with the corresponding secondary antibody for 2 h at room temperature. For spike protein staining, the samples were also incubated with streptavidin, Alexa Fluor™ 594 conjugate (1:200 diluted, Thermo Fisher) to visualize the localization of biotinylated spike protein. Then, cell nuclei were stained with 300 ng/mL-1 4′-6-diamidino-2phenylindole (DAPI) (1:100 diluted, Sigma) for 15 min, and then rinsed with PBS. Phalloidin-FITC (1:200 diluted, Sigma) or Alexa Fluor™ 594 (2 drops per 1mL, Sigma)

Figure 1

 1 Figure 1 Overview of possible ends of shedding spike protein and evaluation process. (A) After endocytosis of Sars-Cov-2, the spike-ACE2 complex may: ① Enter into the cytoplasm; ② Totally or partially shed from the cell membrane; ③ Stay on the cell membrane. (B) Workflow of cardiac patch fabrication and spike protein testing. hiPSCs are first differentiated to cardiomyocytes on a culture patch made of arrayed nanofiber membrane. They are then incubated with synthetic spike protein with or without supplementary Ang II. The beating behavior, cellular structure, and nuclear morphology of the cells are analyzed to evaluate the effect of the spike protein of Sars-Cov-2.

Figure 2

 2 Figure 2 Fabrication and characterization of the cardiac patch. (A) Timeline of hiPSC differentiation toward cardiomyocytes (CM). (B) Beating frequency of the hiPSC-CMs as a function of incubation time. (C) Spatial segmentation of time-lapse images of mature hiPSC-CMs on day 64. (D) Beating waveforms of three types of segmented areas of (C). Over 30 segmented areas, 26 exhibited regular and high magnitude beating (green). (E) Immunostaining images of a selected area showing tissue-like organization with multilayer and multicellular types. The dash lines indicate the borders of the nanofiber membrane. Green: SMA (α-smooth muscle actin), Red: TnT (troponin T), Blue: Nuclei. Scale bar: 20 µm. (F) Schematic of the cell organization. (G) Percentages of the cells expressing TnT and/or SMA. (H) Immunofluorescent images of ACE2 (red) and Connexin 43 (Con43, green) marked hiPSC-CMs. Blue: Nuclei. Scale bar: 5 µm. (I) Maximum intensity projection of Z-stack, fluorescence confocal image of TnT (green) and α-actinin (red) marked hiPSC-CMs. Blue: Nuclei. Scale bar: 5 µm. Zoomed areas display more clearly the sarcomere structures. Scale bar: 5 µm.

Figure 3

 3 Figure 3 Effect of spike protein on hiPSC-CMs. (A) Immunofluorescent images of actin (green) and spike protein (red) marked hiPSC-CMs. Blue: Nuclei. Scale bar: 15 µm. Two zoomed areas show the distribution of spike protein around the nuclei (up) and in the cytoplasm (bottom). Scale bar: 5 µm. (B) Beating frequency variation as a function of incubation time during the first 24 h. Red: sample treated with spike protein,

Figure 4

 4 Figure 4 Beating arrest and nuclear distortion of hiPSC-CMs. (A) Comparison of the rhythmic activity of hiPSC-CMs incubated for 54 h, with or without spike protein and Ang II. 1µg/mL spike protein was added at -3 h, 100 ng/mL Ang II was added after another 3 h. (B) Immunofluorescent images of hiPSC-CMs, showing expression of ACE2 (Green), Actin (Red) and Nucleus (Gray or Blue) under different conditions. Scale bar: 5 µm. (C) The bar graph shows the mean intensity of ACE2 fluorescence under different conditions. (D) Schematic figures to show different nuclear morphology and data statistical analyses, showing quantitatively the changes in chromatin distribution.

Figure 5

 5 Figure 5 Degradation of the cellular structures of hiPSC-CMs. (A) Immunofluorescent images of hiPSC-CMs, showing the expression of TnT (Green), Actinin (Red) under different conditions. Scale bar: 5 µm. (B-C) Bar graphs at left show

Figure 6

 6 Figure 6 Schematic summaries of the effects of spike protein on hiPSC-CMs. (A) Sarcomere structures with and without spike protein treatment. Without spike protein, myosin slides along actin to contract the hiPSC-CMs. With spike protein, the distribution of α-actinin proteins might be disturbed and the overlap between myosin and actin became shortened. (B) Damage effects of the spike protein on cells. Without spike protein, the cell-cell contact with the help of gap junction proteins (connexin 43) enables coordinated contraction of the cardiac tissue. With spike protein, the gap junction disappeared due to mislocalization of connexin 43. When incubated with supplementary Ang II, the spike protein-induced changes became fatal in terms of ACE2 downregulation, chromatin margination, Connexin 43 expression, and sarcomere disruption, leading to cell death.

Figure 3 Figure 4 33

 34 Figure 3

Figure 5 Figure 6

 56 Figure 5

Acknowledgment X.H has been supported by a fellowship from China Scholarship Council. This work was supported by the Agence Nationale de la Recherche under contract No ANR-17-CE09-0017 (AlveolusMimics), DIM ELICIT program from Région Ile-de-France (AutoStemLab), and PSL-Valorization through Pre-maturation project, and IPGG Carnot.

Author contributions

X.H. performed the device preparation, hiPSC differentiation, and data acquisition, J.P and Y.C. were involved in planning and data analyses. All authors discussed the results and commented on the manuscript.

Declarations

Conflicts of Interest. No conflict of interest exists for the authors.