## 2 Materials and Methods
Based on the stretches in the fibre and sheet directions, and the amount of shear in the material point of the continuum and the normal and shear stresses are calculated directly from the curves presented in Figure 2. Then, the tangential moduli \(\overline{E}_{T1}\) and \(\overline{E}_{T2}\) are calculated based on the stress increment with respect to the strain increment within the time step, and finally, the constitutive matrix C is calculated as \(C_{ii}=\overline{E}_{Ti},i=1,2,3;C_{ii}=\overline{G}_{Ti},i=4,5,\ldots,9\) where \(\overline{G}_{Ti}\) are the tangent shear moduli. Since, using the previous procedure, loading in one material direction would not produce deformations in the other two coordinate directions, coupling of the local stresses is then performed using the penalty method [28] and a novel formulation for compressible material behaviour [29]. The cardiac muscle is activated during systole, creating mechanical tension in the muscle cells. To determine the active stresses, we use the most common relationship between active stress and calcium concentration provided by [PERSON] and colleagues [32]:  
\[\sigma_{a}=\frac{\left[Ca^{2+}\right]^{n}}{\left[Ca^{2+}\right]^{n}+C_{50}^{n }}\sigma_{\max}[1+\eta(\lambda-1)] \tag{4}\]  
where details about the coefficients used in Equation (4) are given in [32], while \(Ca^{2+}\) is the calcium concentration which can be used directly from Figure 1 for different drugs. The main advantage of [PERSON]'s model is that it is computationally efficient, which makes it more suitable for practical use in comparison to biophysical models. While it can take several hours to run a finite element simulation with a biophysical muscle model at the micro level, it only takes a few minutes to run the same simulation with [PERSON]'s muscle model. Solid deformations are calculated using the force balance equation [28], and the final form of the incremental-iterative equation for the force balance of the finite element within the solid domain is then written as [27]  
\[\left(\frac{1}{\Delta t}\mathbf{M}+\overline{\mathbf{K}}^{(i-1)}\right)\Delta \mathbf{V}^{(i)}=\mathbf{F}^{ext(i)}-\mathbf{F}^{int(i-1)}-\frac{1}{\Delta t} \mathbf{M}\Big{(}\mathbf{V}^{(i-1)}-\mathbf{V}^{t}\Big{)} \tag{5}\]  
where, in addition to the explanation provided for Equation (3), \(\mathbf{F}^{ext(i)}\) is a vector of the external nodal forces acting on the element at the iteration \"\(i\)\", and \(\mathbf{V}^{(i-1)}\) and \(\mathbf{V}^{t}\) are the nodal velocities at previous iteration \"\(i-1\)\" and at the start of the time step \(t\), respectively. The vector of the internal nodal forces \(\mathbf{F}^{int(i-1)}\) is a function of the stresses calculated from the constitutive curves presented in Figure 2, while the matrix \(\overline{\mathbf{K}}^{(i-1)}\) takes into account the tangent constitutive matrix \(\mathbf{C}_{T}^{(i-1)}\) based on the curves from Figure 2. The details regarding the calculation of passive stresses based on the constitutive curves and computational procedure for solving the coupled motion of the heart are given in [29; 34; 35]. The loose-coupling FSI interaction procedure described in [33] is used for the coupling of fluid flow and solid mechanics. Fluid and solid fields are modeled separately using two different solvers: computational fluid dynamics (CFD) and computational solid dynamics (CSD); see Figure 2c. In this procedure, the same velocity values at the solid-fluid interface are assumed for both components. Within one time step of the numerical simulation, we first solve the equations for the fluid field, and the resulting nodal forces, for an element \(E\), at the fluid-solid boundary are calculated as  
\[\mathbf{F}^{E}=\Big{[}\frac{1}{\Delta t}\mathbf{M}+\mathbf{K}_{vv}\quad \mathbf{K}_{vp}\Big{]}\bigg{\{}\mathbf{V}\bigg{\}} \tag{6}\]  
where matrices \(\mathbf{M}\), \(\mathbf{K}_{vv}\), and \(\mathbf{K}_{vp}\) are the same as from Equation (3).