Fluid mechanics of sarcomeres as porous media

Abstract

Muscle contraction, both in skeletal and cardiac tissue, is driven by sarcomeres, the microscopic units inside muscle cells where thick myosin and thin actin filaments slide past each other. During contraction and relaxation, the sarcomere's volume changes, causing sarcoplasm (intra-sarcomeric fluid) to flow out during contraction and back in as the sarcomere relaxes. We present a quantitative model of this sarcoplasmic flow, treating the sarcomere as an anisotropic porous medium with regions defined by the presence and absence of thick and thin filaments. Using semi-analytic methods, we solve for axial and lateral fluid flow within the filament lattice, calculating the permeabilities of this porous structure. We then apply these permeabilities within a Darcy model to determine the flow field generated during contraction. The predictions of our continuum model show excellent agreement with finite element simulations, reducing computational time by several orders of magnitude while maintaining accuracy in modelling the biophysical flow dynamics.

---

1 Introduction

Many biological processes, ranging from locomotion of microorganisms to large-scale muscle movement, rely on molecular motors, nano-scale machines that generate motion.¹ One such molecular motor is the actin-myosin motor, which converts energy stored in ATP molecules into movement when myosin filaments repeatedly bind (also known as cross-bridging) and unbind to parallel actin filaments. This process is observed as a characteristic “walking” motion of the myosin along the actin, and the two filaments slide relative to one another.^2,3 One notable example of biological components that utilise the actin-myosin motor are sarcomeres, the fundamental contractile units of so-called striated muscle (both skeletal and cardiac), composed of periodic hexagonal arrays of axially-aligned, interdigitating (thin) actin and (thick) myosin filaments.^4,5 Sarcomeres are approximately cylindrical, with a typical diameter of one micron and a typical length of a few microns.^6–9 Many micro-scale sarcomeres arranged end-to-end form a myofibril, and it is the collective periodic contraction and relaxation of many myofibrils that produce large-scale muscle movement (see illustrations in Fig. 1A and B). In this paper, we address the fluid mechanics of sarcomere contraction; we first present in the Introduction a broad overview of the biophysics of sarcomere contraction before detailing the motivation for our modelling study.

Fig. 1 Physical and mathematical modelling of sarcomere. (A) A muscle fibre containing many myofibrils.¹⁰ (B) A single myofibril, containing many aligned sarcomeres. Each Z disk defines the boundary from one sarcomere to the next.¹⁰ (C) Electron micrograph of a single sarcomere, showing the defining banding patterns of striated muscle;¹¹ regions of different intensities correspond to distinct regions in the sarcomere, which are labelled as in D. (D) Simplified structure of a single sarcomere, with thin actin filaments (blue), thick myosin filaments (red), and Z and M disks; the three regions are defined by the presence of the thick and/or thin filaments: regions 1 (thick only), 2 (overlap) and 3 (thin only). (E) Darcy model representation of the sarcomere as an anisotropic porous medium (with the three regions highlighted in different colours), showing radius R and overall sarcomere length 2L. (F) Blow-up of top-right quarter of the Darcy model representation, with non-dimensionalisations giving radius and half-length of 1, and boundary conditions as indicated.

There is great variation in the precise structure of sarcomeres across different species, and even a single organism can contain a range of different sarcomeres.^9,12–15 However, all sarcomeres of striated muscles have the same characteristic banding pattern (see electron micrograph in Fig. 1C): an interdigitating array of thick and thin filaments attached to Z disks (Fig. 1D). The Z disks are thick and dense structures of interconnecting proteins that separate adjacent sarcomeres.^16–18 Each Z disk anchors an array of thin filaments, composed primarily of actin.^19,20 Each sarcomere also contains an array of thick filaments, composed mostly of myosin.^21–25 The midpoints of the thick filaments are anchored to the M disk, which is a thick, multilayered, dense structure of proteins that separates the two halves of the sarcomere^26–29 (this M disk is absent from some sarcomeres, mostly in invertebrates^12,14,15); as we will show below, the presence of M disks in fact has little impact on the fluid flow. We also note that thick filaments continue unbroken through the M disk (if present), with ‘half-filaments' on either side of the M disk; in what follows, we use the term thick half-filament to avoid confusion. Meanwhile, the thin filaments from one sarcomere do not continue unbroken through the Z disk, into the next sarcomere, and so the Z disk has entire thin filaments on either side. An example of the periodic hexagonal arrangement of filaments is illustrated in Fig. 2A, which shows a cross-section of the overlap region of the sarcomere. Upon activation of the muscle, calcium ions are released into the sarcomere from the sarcoplasmic reticulum, an intricate tubular structure that surrounds myofibrils,³⁰ triggering myosin-actin cross-bridging^31–33 and causing the thin filaments (and the Z disks to which they are attached) to be pulled inwards with some speed V, on the order of 1 μm s⁻¹.^34–36 This shortens the sarcomere, seen as contraction of the muscle.^24,37–41

Fig. 2 (A) Sarcomere cross-section in the overlap region (i.e. region 2 in Fig. 1D), here illustrated for a sarcomere of hypothetical radius 200 nm, with a 30° sector of symmetry indicated. (B) 30° sector of cross-section, as indicated in (A). (C)–(E) Hexagonal arrangement of filaments in the overlap region. Packing ratios (number of thin filaments per thick half-filament) of 2 (C), 3 (D) and 5 (E) are shown, with 2 being typical for vertebrates and 3 and 5 often being found in invertebrates.^9,42,43 Fundamental triangles, each being 1/12 of a hexagon, are shown at the bottom. (F)–(H) Fundamental triangles, highlighted in (C)–(E), used to compute the permeabilities and the traction parameter β.

Despite the many protein filaments (i.e. myofilaments) within the sarcomere, the majority of the space between the M disk and Z disk is filled with the liquid cytoplasm of muscle cells, termed sarcoplasm in this context.⁴⁴ Studies of sarcomeres from various species have shown that inter-filament spacing and sarcomere radius can change slightly over the course of a contraction. However, radial strain proves insufficient to maintain a constant volume.³⁶ Consequently, assuming incompressibility, this implies a flow of sarcoplasm out of the sarcomere during contraction, and back in during relaxation. It has recently been demonstrated⁴⁵ that this fluid flow could augment the transport of substrates such as calcium ions and ATP necessary for sarcomere function; see also our preliminary results in Appendix F.

The many myofilaments of the sarcomere can be modelled as long, slender cylinders, and provide complicated but periodic obstructions to fluid flow. The explicit flow past these filaments can be calculated numerically, such as by the finite element method.⁴⁶ However, this not only involves designing an appropriate explicit model, but also requires specialised computational modelling or software, and significant computational power and time to produce accurate results. Furthermore, the process must be repeated whenever considering a different type of sarcomere, or the same sarcomere at a different stage of contraction.

In this paper, we propose an alternative modelling approach and derive a theoretical model of fluid flow in and out of the simplified sarcomere of Fig. 1D that circumvents the need to consider the myofilaments individually and the concomitant small-scale fluid flow. The model instead considers the averaged large-scale fluid behaviour over many filaments, treating the sarcomere as an anisotropic porous medium, as summarised in Fig. 1E. The results of this model lead to excellent agreement with full numerical computations that resolve all the sarcomere length scales, but with vastly reduced computational cost, time and complexity.

This paper is split into two main sections, one devoted to the physical and mathematical modelling (Section 2), and a second one discussing the results for flows inside sarcomeres and comparing our model with full numerical simulations (Section 3). In Section 2, we justify some of our physical assumptions, and then set up the geometry of the model and non-dimensionalise the problem. We next derive and solve the Darcy equations for the pressure and fluid flow in an anisotropic porous medium in each of the three regions of the sarcomere, subject to the relevant boundary and interface conditions. These solutions involve numerical physical parameters (namely, the permeabilities and the traction parameter), whose values we determine by calculating the small-scale fluid flow through the periodic cells of the filament lattice, via semi-analytic methods. We next compute in Section 3 the large-scale fluid flow of the Darcy model, and compare with results of full numerical computations, finding excellent agreement in several orders of magnitude less time. We then compare results of the Darcy model between different sarcomeres, identifying common flow properties, as well as differences. We conclude with a discussion (including potential improvements) of the modelling approach, as well as potential applications, in Section 4.

2 Modelling flow within a sarcomere

2.1 Physical assumptions and simplifications

Our model is derived using two simplifying physical assumptions. We first model the sarcoplasm as a Newtonian fluid. Previous observations indicate that the cytoplasm of muscle fibres is approximately Newtonian for all but the fastest and highest frequency contractions.^36,47–49 However, since the non-Newtonian behaviour of cytoplasm is primarily due to the cytoskeleton,^50,51 which is absent from the interior of the sarcomere, the behaviour of sarcoplasm likely remains Newtonian even at such extremes. Nano-scale objects in the sarcoplasm, such as enzymes,⁵² are likely insufficient, both in abundance⁵³ and individual size, to cause non-Newtonian behaviour.^54–56

Secondly, in the model we take the sarcomere radius (and the inter-filament spacing) to be constant. While it is well established that sarcomere radius (and consequently filament spacing) is not constant throughout contraction, radial strains (i.e. fractional changes in cross-sectional area) are generally smaller than axial strains (i.e. fractional changes in sarcomere length).³⁶ The quadratic dependence of cross-sectional area on radius, and the fact that sarcomere radius is typically smaller than length (see Table 1^6,57,58), indicate that absolute changes in radius are much smaller than absolute changes in length. Furthermore, the radial strain profile over the course of a contraction varies between species, and the Poisson ratio (i.e. the ratio of radial expansion to axial compression) can vary between positive and negative values,³⁶ making precise and general modelling difficult. In terms of fluid flow (see Section 2.3), radial dilation of a matrix of fibres cannot drive any net fluid flux due to overall volume conservation, and we found that the changes to the permeability values (and especially their ratios) and the traction parameter are negligible, at most a few percent, for biologically observed radial strains.³⁶ Therefore, it is appropriate to assume that the sarcomere radius, and inter-filament spacing, remain constant. In other words, we model the filaments as being rigid, not expecting this simplification to have any significant effect on fluid flow.

Table 1 Physical parameters (length scales) used across all sarcomere models, as found in ref. 4, 8, 9, 19, 46, 59 and 60; the thick filament spacing is given by the D₁₀ spacing (horizontal spacing between two adjacent columns of thick filaments in Fig. 2) multiplied by . All values are given to 1 or 2 significant figures

Parameter	Symbol	Value (nm)	Ref.
Thick filament radius	R thick	7	4, 46 and 59
Thin filament radius	R thin	5	19 and 46
D 10 spacing	D 10	45	9, 46 and 60
Thick filament spacing		52
Thick half-filament length	L thick	800	9 and 46
Thin filament length	L thin	1000	9 and 46
Sarcomere radius	R	500	8

---

2.2 Setup and non-dimensionalisation

Following these assumptions, we aim to solve for the flow inside the model sarcomere illustrated in Fig. 1E. Our analysis, unless stated otherwise, will assume a dimensionless sarcomere described by cylindrical polar coordinates (r,z,θ), denoting the radial, axial and azimuthal coordinates respectively, where axial (z) lengths are scaled with the sarcomere half-length L, radial (r) lengths with the sarcomere radius R and axial velocities with the contraction speed V. As a result, the dimensionless sarcomere half-length, radius and contraction speed are all equal to 1. The standard incompressibility equation ∇·u = 0 implies a scaling for lateral velocities of RV/L. Further, pressure is scaled with the viscous scaling μV/L, where μ is the dynamic viscosity of the fluid. We may finally exploit z → −z symmetry to solve the problem on a half-sarcomere. In this dimensionless problem, region 1 (thick filaments only) is located in 0 < z < L₁, region 2 (overlap region) is located in L₁ < z < L₂, and region 3 (thin filaments only) is located in L₂ < z < 1.

We will show in our analysis that the porous medium exhibits lateral isotropy (all flow in a cross-section has the same permeability, regardless of direction) which, coupled with the cylindrical geometry of the sarcomere, leads to azimuthal symmetry, and thus we consider the problem in the (r,z) polar plane. The resulting dimensionless setup is shown in Fig. 1F (which also shows boundary conditions, that we discuss below). Note that the three regions of the sarcomere (separated by interfaces at z = L₁ and z = L₂) are considered separately: we solve for the pressure and fluid flow in each region individually, and apply interface conditions to relate them.

2.3 Darcy flow modelling

2.3.1 Introduction to Darcy flow. Within sarcomeres, there is a clear separation of length scales between, on one side, the thick and thin filaments (typical width 10 nm) and the spacing between the various filaments (typical size 10's of nm) and, on the other, the sarcomere itself (typical width and length of 1 μm). Motivated by this observation, we model in our paper the sarcomere as a porous medium, i.e. a medium that is not simply free fluid, but rather consists of fluid (occupying a volume fraction Φ, termed the porosity) that must navigate fixed physical obstructions (here, the myofilaments, occupying a volume fraction 1 − Φ). The fluid/obstruction mixture is considered a continuum, which is valid assuming that the length scale of flow is much larger than the pore scale (the characteristic size and separation of the obstructions), as justified above. The resultant flow problem is known as Darcy flow, a statistical average of fluid flow over many pore scales.^61,62

In our sarcomere model, we deviate from traditional Darcy flow in two ways. First, Darcy flow typically considers isotropic porous media, but here, the alignment of the myofilaments along the long axis of the sarcomere makes the medium anisotropic, with axial flow typically easier than lateral flow. Second, most porous media consist of pores that exhibit randomness, with their averaged properties determined statistically. Here, instead, the porous medium consists of regular periodic cells, and we perform averages over individual cells rather than many pore scales. Darcy flow is concerned not with the fluid velocity, but rather the fluid flux (often called the Darcy flux or Darcy velocity), which is the rate of fluid flux per unit area through the porous medium. Pressure also exists as a spatial average. To avoid confusion, we thus reserve throughout the terms u and p for conventional, interstitial (fluid) velocity and pressure, respectively, which undergo sharp changes within a periodic cell, and U and P for the Darcy flux and Darcy pressure, which eliminate these sharp changes via spatial averaging, with U = Φ 〈u〉 and P = 〈p〉.

2.3.2 Darcy fluxes in the anisotropic porous medium. In order to identify the large-scale flow behaviour in the sarcomere, we first consider the small-scale flow between the filaments and properly derive the Darcy model. To do so, we focus on the individual periodic hexagonal cells that comprise the larger periodic lattice, as shown in Fig. 2C–E. By neglecting variation in the flow between adjacent periodic cells, we exploit various symmetries to calculate semi-analytic expressions for the fluid flow within a single periodic cell, and obtain explicit linear relationships between fluid flux and pressure gradient. Details of the derivations, including necessary assumptions, are given in Appendix A (axial permeabilities), B (traction parameter β), C (lateral permeabilities), but essentially we need only assume that the length scales (both axial and lateral) over which ∇P varies are much larger than the size of each individual hexagonal cell of myofilaments.

By solving for the axial fluid flow in a hexagonal cell, we calculate an overall axial fluid flux per unit area of

for some constant axial permeability k_‖ that is a function only of the geometry of the hexagonal lattice and which of the three regions we are in (thin filaments vs. thick filaments vs. overlap region); the various axial permeabilities are calculated in Appendix A.

We now turn our attention to the lateral fluid flow, that is, fluid flow within a cross-section of the cylindrical sarcomere, perpendicular to the z axis. In general, the resultant flux can be expressed in terms of a 2 × 2 lateral permeability matrix k_⊥, where ⊥ indicates the plane perpendicular to the z axis:

U_⊥ = −k_⊥∇_⊥P, (2)

where ∇_⊥P denotes the 2-dimensional gradient of P in the (r,θ) plane, perpendicular to the z axis. Once again, we calculate k_⊥ by deriving the explicit lateral fluid flow within a hexagonal cell. Note, however, that the hexagonal cells, in Fig. 2C–E, exhibit a six-fold rotational symmetry. In particular, this means that a lateral pressure gradient ∇_⊥P (of fixed magnitude) in the direction normal to any given side of the hexagons will yield the same magnitude of lateral fluid flux per unit area, i.e. Darcy flux, and this flux will be directed in the same direction as the applied pressure gradient. Importantly, two such (independent) directions form a basis of the plane. We further note that the linearity of the Stokes equations allows us to reconstruct any flow (including the corresponding pressure field) as a linear combination of said basis. We consequently deduce that the lateral Darcy flux is equal in magnitude, and parallel to ∇_⊥P, regardless of the direction of the applied pressure gradient. We therefore calculate a lateral flux per unit area U_⊥ = −k_⊥∇_⊥P for some constant scalar lateral permeability k_⊥ (here also only a function of local lattice geometry and which of the three regions of interest is being considered), which we calculate in Appendix C.

We conclude that the Darcy model of the sarcomere has cross-sectional isotropy, at which point the cylindrical geometry of the entire porous sarcomere implies an overall Darcy fluid flow that is independent of the azimuthal angle θ and has no component in the θ direction, but is instead purely axial and radial. We therefore obtain the Darcy equations for flow in a sarcomere as a porous medium as

where we have used superscript j to refer to the three regions of the sarcomere, and no Einstein summation is implied.

The result in eqn (3) is the fluid flux that arises solely from pressure changes, but it does not account for the fluid flow induced directly by the thin filaments moving relative to the thick ones during contraction or relaxation. This movement has no effect on any of the lateral fluid fluxes U_r, nor on the axial flow in region 1, U_z⁽¹⁾. In the thin region, the absence of stationary thick filaments means that the moving thin filaments pull the entirety of the fluid along with them, and hence U_z⁽³⁾ is decreased by a value equal to the porosity Φ⁽³⁾. In the overlap region, the stationary thick filaments partially resist the flow induced by the moving thin filaments, and U_z⁽²⁾ decreases by a value βΦ⁽²⁾, where the traction parameter β is a dimensionless constant between 0 and 1 that is determined entirely by the geometry of the cell, and calculated similarly to the axial permeabilities, in Appendix B. Together with eqn (3), this gives us modified, final Darcy equations

where γ⁽¹⁾ = 0, γ⁽²⁾ = β, γ⁽³⁾ = 1 capture the effects of the moving thin filaments in each region.

2.3.3 Pressure equation. Since the fluid flow is incompressible, it follows that the same is true for the fluid flux and thus ∇·U = 0. Substituting the two Darcy equations from eqn (4) into the incompressibility condition leads to the governing equation for the Darcy pressure field P as

This is accompanied by relevant boundary conditions, as discussed below. Note that although the interstitial pressure p is harmonic, as required by Stokes flow, the Darcy pressure P is not exactly harmonic; however, by a suitable rescaling of z, we would recover Laplace's equation and so, with appropriate boundary conditions, we deduce that P has a unique solution.

2.3.4 Boundary conditions for the full sarcomere model. We first identify appropriate boundary conditions for the full sarcomere model of Fig. 1D, so that we can approximate them in the Darcy model, and later perform full numerical computations. We may exploit axial symmetry to consider only the right half-sarcomere z ≥ 0. The surfaces of the individual filaments are subject to no-slip conditions, i.e. the fluid velocity u is equal to the velocity of the filament. The Z disk and M disk (if present) are protein structures that are much more densely packed than the three filament regions, so it is legitimate to approximate a no-slip condition on the disks as well. Alternatively, when the M disk is absent, we should apply a symmetry condition at z = 0: normal velocity is zero, u_z = 0, and normal gradient of tangential velocity is also zero, ∂u_r/∂z = ∂u_θ/∂z = 0. We model the sarcomere as being immersed in free fluid, so in our simulations we also apply a zero-stress far-field condition away from the sarcomere.

2.3.5 Boundary and interface conditions for the Darcy model. We now determine the appropriate boundary conditions, based on the above, for the Darcy model of the sarcomere. Being a rescaled Laplace equation, eqn (5) requires a single boundary condition at each boundary. We set P = 0 at r = 1 to simulate the stress-free far-field condition. In particular, this can be understood by noting that the dense internal structure of the sarcomere is expected to generate much larger pressure gradients inside the sarcomere compared to outside, differing by a factor of (R/l)², where R ([scr O, script letter O] (1 μm)) is the sarcomere radius and l ([scr O, script letter O] (10 nm)) is the pore scale, according to scalings of the Stokes equations. This assumption is confirmed by the full numerical simulations discussed above. We apply no-penetration, ∂P/∂z = 0, at the disks. This may produce a non-zero tangential slip velocity, in contrast to the full boundary conditions; however, we will see that the overall effect of this turns out to be very small on the scale of the full sarcomere. Meanwhile, if the M disk is absent, the two symmetry conditions are u_z = 0 and ∂u_r/∂z = 0, both of which can be achieved by setting ∂P/∂z = 0. In other words, the Darcy flow is the same regardless of whether the M disk is present or not. Finally, there is an additional regularity condition ∂P/∂r = 0 at r = 0 that arises from exploiting azimuthal symmetry in the cylindrical geometry.

Writing these explicitly, we therefore need to solve eqn (5) subject to

P⁽¹⁾ = P⁽²⁾ = P⁽³⁾ = 0 at r = 1, (7)

There remains the task of identifying the interface conditions between adjacent regions of the sarcomere. Assuming that the Darcy equations hold everywhere, we identify the change in pressure across an interface by performing a force balance. We integrate the axial flux equation in a small region about the interface at z = L_k, giving

Taking ε → 0, we see that the change in pressure across the interface is zero; continuity of pressure is a standard condition for interfaces in porous media.^62,63 Meanwhile, volume conservation gives continuity of axial flux at the interface between regions 2 and 3 (z = L₂). In contrast, the interface between regions 1 and 2 (z = L₁) is more subtle. As the sarcomere contracts, the interface between regions 1 and 2 moves leftwards. This moving boundary produces an effective flux source, and the resulting increase in axial fluid velocity when moving from region 1 to region 2 is precisely (1 − Φ⁽³⁾), which is the volume fraction occupied by the thin filaments. The interface conditions are therefore

P⁽¹⁾ − P⁽²⁾ = 0 at z = L₁, (13)

P⁽²⁾ − P⁽³⁾ = 0 at z = L₂. (14)

2.3.6 Solution for pressure. The solution to eqn (5) in region j (with j = 1, 2, 3) that satisfies the boundary conditions in eqn (6) and (7) is obtained classically as the serieswhere J₀ is the zero order Bessel function of the first kind, λ_n is the nth zero of J₀, and . This solution satisfies the boundary conditions eqn (6) and (7) analytically; the two remaining boundary conditions eqn (8) and (9), as well as the four interface conditions eqn (11)–(14), must be satisfied by appropriate choice of coefficients A^(j)_n and B^(j)_n. Since there are three regions, this amounts to 6 coefficients to be determined for each value of n. By using the orthogonality condition between Bessel functions, we may apply the boundary and interface conditions eqn (8)–(14) to obtain an appropriate 6 × 6 matrix equation for the six coefficients above, for each value of n independently (see Appendix D). Once we calculate the permeabilities and the traction parameter β, we easily solve these matrix equations numerically, and hence obtain all the coefficients. Of course, in order to compute the solution, eqn (15), we must truncate the series to some finite number of terms N,where, to reiterate, the permeabilities are incorporated via, with λ_n being the nth zero of the Bessel function J₀. In all results below, we set N = 50 (see Appendix D for justification).

2.4 Biophysical parameter values in biological sarcomeres

In order to validate the Darcy model, we must compare it with precise numerical computations. To this end, we identify here suitable values of the various biophysical parameters. We note that there is a wide range of biological parameters, observed in a variety of sarcomeres; below we will test the Darcy model in a minimal suite of configurations, for which we expect the worst agreement, and will next discuss the extension to more ideal configurations.

2.4.1 Fixed value parameters: filament radii and spacing. A review of the existing literature, examining both intact sarcomeres and isolated actin and myosin filaments, reveals that certain biophysical parameters are approximately constant across all sarcomeres.^{4,9,19,46,59,60} The radii of both thick and thin filaments, as well as the D₁₀ filament spacing, D₁₀, which is defined as the horizontal distance between two adjacent columns of thick filaments in Fig. 2, will be assumed to be constant. Consequently, the distance between any two adjacent thick filaments, given by , is also constant. These values are listed in Table 1.

Geometrically, the thick and thin filament radii are dictated by the fixed size of the individual myosin and actin monomers, respectively, whilst the filament spacing is dictated by the distance over which myosin molecules can reach during the cross-bridge cycle; we can therefore also consider these parameters to be constant. Structural variation between sarcomeres, such as different packing ratios, can affect these values to an extent; we selected values, consistent with the larger body of literature, to match those used in previous numerical work.⁴⁶

In what follows, we will also set the dimensional muscle contraction speed V to 1000 nm s⁻¹, which is the observed order of magnitude for sarcomere contraction.^34–36 Obviously, due to the linearity of the (viscous) physical system, all flow will instantaneously scale linearly with the chosen value of V.

2.4.2 Extreme value parameters: filament lengths and sarcomere length; sarcomere radius. The remaining parameters could be selected from a wide range of biologically observed values. Filament lengths have been reported to vary through a factor of at least around 5, though there is a systematic linear relationship between thick and thin filament lengths.⁹ Sarcomere radius can also vary significantly.⁸ The Darcy model is expected to be most valid when there are large relative length scales of variation, and in order to validate it we will focus on the worst-case scenario, i.e. situations of shortest observed length scales. We thus set sarcomere radius and filament lengths (and consequently overall sarcomere length) to take their smallest observed values; these have been added to Table 1.

2.4.3 Variable parameters: contractions; M disk; packing ratio. With all biophysical parameters set, we vary in our investigation the level of contraction, α ≡ 1 − L₁/L₂, between states of extreme contraction (α large) and extreme relaxation (α small), thereby considering the full range of possible values for the lengths of the three regions. We expect states of extreme contraction or relaxation, both of which lead to the presence of short regions, to produce the least accurate agreement with the full numerical computations. In reality, sarcomeres do not typically contract over the full range, although contractions of hundreds of nanometres are observed even in small sarcomeres.^34,35 We also consider the presence or absence of the M disk, and cases of differing packing ratios.

2.4.4 Permeabilities and traction parameter β. Thanks to the aforementioned constant filament radii and spacing, the traction parameter β and the six (dimensional) Darcy permeabilities are constant for each packing ratio. While β is a truly dimensionless parameter, the dimensional permeabilities each have dimensions of length squared, and must be appropriately non-dimensionalised to be used in eqn (16). Specifically, each dimensional k^(j)_‖ (being a property of axial flow) should be divided by L², and each dimensional k^(j)_⊥ (being a property of lateral flow) by R², for dimensional sarcomere length L and radius R. Therefore, despite each dimensional permeability being independent of L and R, their non-dimensional counterparts in eqn (16) are not. As such, in order to present the permeabilities in a succinct, dimensionless form, that does not depend on L or R, we introduce new non-dimensionalisations specifically for the triangular cells of Fig. 2F–H within which the permeabilities are calculated. All lengths are scaled with the triangle height, , so that the triangles now have a dimensionless height of 1. The resultant dimensionless permeabilities have been computed in Appendix A and C, with β being calculated in Appendix B, and these are listed in Table 2 for packing ratios of 2 (left column), 3 (middle) and 5 (right). To restore the permeabilities to their dimensional values, they must be multiplied by D₁₀²/3. In order to use the permeabilities in eqn (16), each k^(j)_‖ must then be divided by L² and each k^(j)_⊥ must be divided by R². Therefore overall, to use the permeabilities in Table 2 within eqn (16), each k^(j)_‖ must be multiplied by D₁₀²/3L², and each k^(j)_⊥ must be multiplied by D₁₀²/3R². In contrast, β is a truly dimensionless quantity that requires no such treatment.

Table 2 Computed dimensionless axial and lateral permeabilities in the three regions, as well as the traction parameter β, for packing ratios of 2, 3 and 5, using a triangle with height 1, see Fig. 2 and Appendix A, B, C. To use these values in eqn (16), each stated value of k^(j)_‖ must be multiplied by D₁₀²/3L², and each k^(j)_⊥ must be multiplied by D₁₀²/3R², where L and R are the dimensional sarcomere length and radius, respectively (see main text and Appendix A, C). All geometrical parameters are taken from Table 1

Variable	Packing ratio
	2	3	5
k ‖ ^(1)	0.3729	0.3729	0.3729
k ‖ ^(2)	0.0728	0.0433	0.0247
k ‖ ^(3)	0.2038	0.1071	0.0759
β	0.6175	0.7127	0.7238
k ⊥ ^(1)	0.1864	0.1864	0.1864
k ⊥ ^(2)	0.0362	0.0212	0.00618
k ⊥ ^(3)	0.0938	0.0435	0.00737

---

3 Fluid flux in sarcomeres

In this section, we apply our modelling approach to compute the fluid flux in and out of contracting sarcomeres. Specifically, we compute the predictions of our Darcy model and compare them with results of full numerical computations obtained using the finite element method (FEM) within COMSOL Multiphysics,⁶⁴ wherein the entire geometry of the sarcomere is fully specified and the flow past the array of filaments is explicitly calculated (see Appendix E for details).

It should be emphasised that the two approaches (Darcy vs. FEM) are fundamentally different, with the latter considering conventional fluid properties, and the former considering only the bulk fluid properties; as such, we will need to average the FEM data in order to compare with the Darcy model. Note that the comparisons will be made using dimensional variables, to help with application to biological systems. As we will see, the Darcy model is able to accurately reproduce the bulk fluid flow of the precise FEM computations in an averaged sense, serving as a highly accurate, and far more time-efficient and less resource-intensive, method of calculating the flow within the sarcomere.

3.1 Illustrative example

3.1.1 Fluid flow. We start our comparisons between the two models by focusing on the particular case of a sarcomere with a packing ratio of 5 in the absence of an M disk (typical of invertebrates^9,12,14,15), at 50% contraction, i.e. α = 1 − L₁/L₂ = 0.5. To carry out a detailed comparison between the Darcy and FEM approaches, we azimuthally average the FEM data and plot the flux magnitude as a heat map, with arrows indicating direction of flux. The results are shown in Fig. 3 with Darcy predictions in Fig. 3A, and FEM numerics in Fig. 3B.

Fig. 3 Comparisons of fluid flow between the Darcy and FEM models, for a sarcomere at 50% contraction with a packing ratio of 5 in the absence of an M disk. (A) and (B) Heat map of fluid flux with arrows indicating direction for Darcy (A) and FEM (B) models. Note that in A, the maximum flux magnitude (found at the radial boundary, at the interface between regions 1 and 2) is approximately 1600 nm s⁻¹, whereas the colour scale has a maximum value of 1100 nm s⁻¹; this is to improve readability of A and B, and less than 0.1% of the domain surpasses 1100 nm s⁻¹. (C) Plot of radial flux against z for the Darcy (red) and FEM (blue) models, with r = 100, 200, 300, 400, 500 nm from bottom to top. Radial fluxes have been individually rescaled to have the same maximum value of radial Darcy flux, to improve readability.

We observe excellent agreement between the two models, in terms of the magnitude, direction and spatial distribution of the flow. The striations in the FEM model, most visible in the overlap region, correspond to the presence and absence of filaments; these filaments cause fluctuation in the axial fluid flux (but not the radial fluid flux) which are effectively averaged out by the Darcy model.

3.1.2 Radial flux. Since the heat maps in Fig. 3A and B alone do not offer a proper quantitative comparison between the two models, we next compute the rescaled radial flux against z at a number of fixed values of r and compare with the average of the computational results in Fig. 3C. The Darcy model accurately captures the radial fluxes across nearly all values of r.

Some small quantitative differences can be observed. For very small r, the model deteriorates somewhat; this represents only a few percent of the entire sarcomere volume and is not of practical concern. Some discrepancies are also seen for very large r, as the periodic pattern begins to terminate. Specifically, whilst the Darcy model considers a porous medium in a perfect cylindrical geometry, the FEM model explicitly considers a hexagonal lattice that must be fitted inside said cylinder. As such, the lattice must be terminated as it approaches the outer radial boundary. As a result, within the FEM model, there are small regions near the outer radial boundary but still within the cylinder for which the hexagonal lattice of filaments is absent.

Considering variation of z, the most noticeable discrepancies between the porous media approach and the full simulations appear near the interfaces between the regions of the sarcomere, where the assumptions necessary for Darcy flow break down, and near the no-slip condition at z = L, which cannot be precisely satisfied by the Darcy model. These create short regions of disagreement, across thicknesses that are about the size of the hexagonal cells.

3.1.3 Axial flux. We next consider the axial fluxes through the sarcomere. Given the very good agreement of radial flux illustrated in Fig. 3C, it follows from mass conservation that the axial fluid fluxes should also agree in an averaged sense. This is confirmed directly in Fig. 4, where the azimuthally averaged FEM axial flux and the Darcy axial flux are plotted as functions of r for select values of z. Overall we also find excellent agreement at almost all values of r and z, with notable exceptions as discussed above. In particular, note that the fluctuations of the FEM axial flux caused by the presence of the filaments are averaged out by the Darcy model.

Fig. 4 Comparisons of fluid flow between the Darcy and FEM models, for a sarcomere at 50% contraction with a packing ratio of 5 in the absence of an M disk. Plot of axial flux against r for the Darcy (red) and FEM (blue) models, with z = 100, 300, 500, 700, 900, 1100, 1300 nm from bottom to top. Axial fluxes have been individually rescaled to have the same maximum value of axial Darcy flux, to improve readability.

3.1.4 Lagrangian transport. Given this excellent agreement, we can now use the Darcy model to compute and consider the Lagrangian deformation of fluid particles within the sarcomere. We numerically integrate (via a first-order forward Euler scheme) the interstitial fluid velocity, u = U/Φ, of the Darcy model to determine the deformation of marked fluid elements over the course of a contraction. The results are illustrated in Fig. 5 for sarcomere contractions of α = 5% (A), 25% (B), 50% (C), 75% (D) and 95% (E).

Fig. 5 Lagrangian deformation of marked fluid particles throughout a sarcomere contraction, shown for contractions of α = 5% (A), 25% (B), 50% (C), 75% (D) and 95% (E). Dashed black lines indicate interfaces between the three regions.

We observe that fluids in the thin and overlap regions typically undergo very little deformation. Meanwhile the fluid in the thick region undergoes large displacements and deformations in both the axial and radial directions, with significant axial compression and radial stretching of the fluid parcels, and large radial efflux out of the sarcomere. Due to the reversibility of Stokes flow, this fluid will be drawn back (deterministically) into the sarcomere as it relaxes. This may help to draw in useful substrates (such as ATP) during relaxation.⁴⁵ Furthermore, various metabolic processes occur within the sarcoplasm, involving a variety of substrates with uses and effects both inside and outside the sarcomere.⁵² Fluid flow may beneficially redistribute these substrates, and in particular may help rid the sarcomere of waste products (such as lactate resulting from anaerobic respiration, and excess ADP resulting from hydrolysis of ATP over prolonged activity) during contraction. Interestingly, the axial compression of the fluid is insufficient to keep the fluid parcels, that begin in the thick region, in the thick region; the growing overlap region absorbs much of this fluid. Similarly, much of the fluid that begins in the thin region is pushed into the growing overlap region. Both of these are likely to help advect useful substrates, such as ATP and calcium ions needed for cross-bridging,^31–33 into the overlap region. Importantly, advection of ATP is likely to be highly impactful to sarcomere function, since the internal structure of the sarcomere has been shown to impede ATP diffusion by several orders of magnitude compared to free cytosol.^65,66

3.2 Analysis and comparisons of a variety of sarcomeres

Having established agreement between the Darcy model and full numerical simulations in the above illustrative example, we investigate in this subsection the impact of varying the nature of the sarcomere. Specifically, we show a similar comparison in the case of a sarcomere with a packing ratio of 3 in the absence of an M disk at various stages of contraction in Fig. 6; this is the case characteristic of invertebrate organisms.^9,12,14,15 We also show a similar result for a sarcomere with a packing ratio of 2 in the presence of an M disk in Fig. 7, here a situation relevant to vertebrate sarcomeres.^{9,12,26–29,42,43}

Fig. 6 Radial fluid flux in a sarcomere with a packing ratio of 3 (cross-section of overlap region illustrated on top left) in the absence of an M disk at various stages of contraction, α = 5% (A), 25% (B), 50% (C), 75% (D), 95% (E), showing the Darcy model (red) against the FEM model (blue). Each figure shows five pairs of data sets, at r = 100, 200, 300, 400, 500 nm from bottom to top.

Fig. 7 Radial fluid flux in a sarcomere with a packing ratio of 2 (cross-section of overlap region illustrated on top left) in the presence of an M disk at various stages of contraction, α = 5% (A), 25% (B), 50% (C), 75% (D), 95% (E), showing the Darcy model (red) against the FEM model (blue). Each figure shows five pairs of data sets, at r = 100, 200, 300, 400, 500 nm from bottom to top.

In both cases, we again observe overall excellent agreement between the two models. When the M disk is present (Fig. 7), some discrepancies are seen in its vicinity, all the more present if the thick region is very small (limit of large contractions). Despite this, the agreement between the Darcy approach and the averaged numerics is essentially perfect at most values of r and z.

We are now able to compare the properties of the flow between the different sarcomeres. By examining Fig. 3, 6 and 7, we see that the radial efflux is low within the overlap region, owing to the low permeability caused by the filaments, but is typically greater in the thick and thin regions, a consequence of the higher permeabilities. Generally speaking, the thick region contributes the majority of radial flow out of the sarcomere – the sarcomere contracts, and the thin filaments invade the thick region, thus pushing fluid out of the sarcomere. A similar argument of thick filaments invading the thin region (after an appropriate change of frame) explains the significant radial outflow occurring in the thin region. The radial outflow in the thick region is most significant for higher packing ratios; indeed, as the packing ratio increases, the relative importance of the thin filaments increases. As such, the ‘plunger’ effect of the thin filaments invading the thick region is effectively increased, leading to a more significant radial efflux in the thick region. Simultaneously, this increase in packing ratio reduces the permeability in the thin region, thereby reducing the radial efflux there. Predictably, this ‘plunger’ effect is most evident near the interfaces with the overlap region.

Overall, we see that the Darcy model accurately predicts the radial flux, and hence all fluid properties, for all levels of contractions and sarcomere geometries. This agreement also allows us to confirm, a posteriori, that the methods by which the permeabilities were calculated (values in Table 2) were correct, and hence, again in an averaged sense, the Darcy model accurately determines the pressure and the interstitial flow profiles between the filaments.

4 Discussion

4.1 Summary

The fluid flow past the many interdigitating filaments of the sarcomere requires significant computational resources to calculate accurately. In this paper, we demonstrated how to derive a Darcy model that considers the sarcomere as an anisotropic porous medium, and we compare its predictions with those of full numerical (FEM) simulations.

The model gives a simple governing equation for the bulk pressure P (eqn (5)) with further equations to obtain the bulk fluid flux from this pressure (eqn (4)). The boundary conditions valid in a full computational (FEM) model can be approximated in the Darcy model, and force and flux balances allow for appropriate interface conditions between the regions of the sarcomere. The final solution for P is a truncated series (eqn (16)) whose coefficients are easily calculated from these boundary and interface conditions, by exploiting the orthogonality condition between Bessel functions. In order to do this, we calculate six permeabilities and the traction parameter β (Table 2), which is done by adapting classical semi-analytic schemes to triangles representing 1/12 of a hexagonal cell. The calculated values for the permeabilities and β are then fed into the equations for the coefficients (Appendix D). Having determined these coefficients, we can evaluate P and U at any given point in the sarcomere practically instantly.

Overall, we obtain excellent agreement between the Darcy and FEM models, with minor exceptions near the radial centre and exterior radial boundary of the sarcomere, as well as near interfaces and disks. The Darcy model therefore allows us to calculate, with surprising accuracy, the fluid flow within the sarcomere, in a vastly reduced amount of time. Indeed, the FEM model required a multi-processor machine with specialised software, taking around 24 hours to compute a particular solution. Meanwhile, once the permeabilities are calculated, the Darcy model can be computed using basic software in less than a second on a simple laptop computer; the permeabilities themselves can be calculated in less than a minute on the same device.

4.2 Potential applications

In particular, this allows us to calculate the fluid flow easily in a continuous range of contractions, enabling us to accurately determine the time-evolution of the fluid within the sarcomere. This allows us to study, both quantitatively and qualitatively, the behaviour of the fluid flow, and consider how this may affect the function of the sarcomere, such as by the transport of substrates beyond that achieved by diffusion alone. A demonstration of this phenomenon is illustrated in Appendix F, where preliminary results on substrate advection, reaction and diffusion are shown.

Another potential application is in the calculation of viscous drag or viscous energy dissipation as a result of fluid flow. Importantly, these must be calculated at the pore scale, analysing the semi-analytic flow profile to obtain a linear (though anisotropic) relationship between pressure gradient (or fluid flux) and viscous drag, and a quadratic relationship between pressure gradient (or fluid flux) and viscous energy dissipation. These can then be integrated over the sarcomere to obtain an expression for the total drag force and total rate of doing work. Alternatively, order-of-magnitude estimates can be obtained using scaling arguments, with viscous drag per unit filament surface area on the order of μu/l, and viscous dissipation rate per unit fluid volume on the order of μu²/l². Here u is a characteristic fluid velocity (10⁻⁶ m s⁻¹ (ref. 34–36)) and l is the relevant length scale (i.e. the pore scale, 10⁻⁸ m). Even using an upper estimate for sarcoplasm viscosity (being a few tens of times that of water, [scr O, script letter O] (10⁻²) Pa s^67,68) we find that the viscous drag experienced by each filament is at most [scr O, script letter O] (10⁻¹⁴) N. This figure was also classically obtained by Huxley, who additionally noted that the total viscous drag experienced by the filaments in the sarcomeres is several orders of magnitude smaller than the forces muscles generate, and so is unlikely to be of consequence.⁶⁹ Viscous energy dissipation rate per unit volume can be estimated either using the aforementioned direct scaling argument, or simply by multiplying the viscous drag force experienced by each filament by the contraction speed, and accounting for the number of filaments per unit volume, and is found to be at most [scr O, script letter O] (10²) W m⁻³. This is several orders of magnitude smaller than the rate of energy metabolism from ATP hydrolysis, being around 10⁵ W m⁻³,^52,70 and so is also unlikely to be of consequence.

4.3 Versatility of semi-analytic solutions for flow within periodic cells

The semi-analytic methods of calculating the fluid flow, and hence the fluid flux and permeability, through the periodic cells are highly adaptable. For example, whilst only three possible packing ratios are shown in Fig. 2, packing ratios of 4, 6 and even 12 have been reported in biological sarcomeres.⁹ It is a simple matter to adapt the above methods to these more exotic packing ratios, simply by placing thin filaments appropriately along the lower edge of the triangle (Fig. 2F–H). Furthermore, the general method can be adapted to other arrangements of filaments, such as a square arrangement (e.g. as considered by Sangani & Acrivos⁷¹) by splitting the traditional periodic cell into appropriate right-angled triangles. In the case of a square arrangement, these would be isosceles right-angled triangles. Since only one side of the triangle is exposed to the boundary of the larger periodic cell (the other two are simply lines of symmetry within the cell) only that boundary can contain auxiliary filaments, making it a theoretically easy (though sometimes numerically cumbersome) matter to alter the numerical boundary conditions as appropriate in the presence or absence of auxiliary filaments. Lastly, there is nothing in our analyses that demands the thin filaments are of equal size, or that all of the thin filaments are moving when calculating the traction parameter β. This is not relevant for sarcomeres, but could prove useful in other situations requiring calculation of the flow past a complex periodic array of (active) cylinders.

4.4 Limitations of the model sarcomere

The model studied in this paper (Fig. 1D) is a mathematical and physical simplification of reality that allows us to demonstrate the effectiveness of the Darcy approach to muscle flows, and a number of potential improvements could be proposed.

Whilst myosin and actin are the two most abundant proteins of the sarcomere, composing approximately 50% and 20% of myofibrillar protein mass respectively,^72,73 we did not include the giant protein titin found in vertebrates, which makes up around 10% of the myofibrillar protein mass.^12,74 Titin has a diameter of around 4 nm, and each titin molecule runs from the M disk to the Z disk, via the thick filaments.^75–78 Titin can also bind to actin in certain places, though some of these are uncertain,^79–86 and it is unclear what precise shape the titin molecules take in the thin region. Adding this to the fact that there are believed to be six titin molecules per myosin half-filament,⁸⁷ we see that the presence of titin may have a non-trivial effect on fluid flow within the thin region, and this could be the basis of future work.

A further complication in the specific case of vertebrate sarcomeres is that the thin region is often not a neat continuation of the hexagonal geometry in the overlap region. The thin filaments can become disorganised within the thin region, and appear to ultimately bind with the Z disk in a square arrangement.^12,88,89 The resultant uncertainty in the permeabilities within the thin region motivates a demonstration of how the model can be applied to a continuous range of permeabilities, which could be measured experimentally, without the need to resolve the precise physical structure within the thin region; we refer the reader to Appendix G for further details.

We finally note that in this work we do not focus on cardiomyocytes (i.e. heart muscle cells) that contain branching myofibrils whose radii can be as small as 200 nm,^90–93 significantly smaller than the 500 nm considered in Table 1. This smaller length scale could impact the validity of the Darcy model, and accounting for this should be the focus of future work.

We hope that the simplicity of our Darcy approach, and its agreement with full simulations, will encourage the further development of porous media models for muscle hydrodynamics.

Data availability

The MATLAB code and processed COMSOL data needed to compute and verify the Darcy flow, and produce Fig. 3–7, have been made freely available on GitHub.⁹⁴

Conflicts of interest

There are no conflicts of interest to declare.

Appendices

A Calculating axial permeabilities

In this first Appendix, we apply and adapt a semi-analytic method to calculate the axial flow induced by a constant axial pressure gradient.⁹⁵ The cross-section of the sarcomere seen in Fig. 2A, of hypothetical radius 200 nm, consists of a hexagonal filament lattice, as displayed in Fig. 2C–E, which shows packing ratios of thin actin filaments to thick myosin half-filaments of 2, 3 and 5 respectively; note that this illustrates the lattice within the overlap region, where both thick and thin filaments are present. The lattices in the thick and thin regions are obtained simply by removing the appropriate filaments. Therefore, each region exhibits a hexagonal lattice. These hexagonal cells can themselves be split into 12 right-angled triangles, seen in Fig. 2F–H respectively, and reproduced for a packing ratio of 5 in Fig. 8. Note that, for all the packing ratios presented, and any region of the sarcomere, the filaments present as sectors occupying the vertices of these triangles.

Fig. 8 Triangular domain, 1/12 of a periodic hexagonal cell, on which the axial and lateral flow profiles are calculated, shown here for the overlap region with a packing ratio of 5 (see Fig. 2). The triangle's original height of is rescaled (using a different non-dimensionalisation to that used for the Darcy model) to 1. Polar coordinates (r,θ), distinct from those used for the Darcy model, are shown, with rescaled thick and thin filament radii R_thick and R_thin indicated.

These triangles within the lattice have dimensional height , but to ease these computations, we rescale the triangle to have a height of 1. We also introduce temporary, local polar coordinates (r,θ), which notably are distinct from the cylindrical polar coordinates used for the full sarcomere, and these have been added to Fig. 8. Assuming that the pressure gradient is entirely axial (i.e. along z), and varies over length scales much greater than the size of the triangle, we calculate the flow w driven by a uniform axial dimensionless pressure gradient of 1 by solving the resultant dimensionless Stokes equation in the triangle

We next apply a substitution w* = w − r²/4 and solve the resultant Laplace equation in plane polar coordinates

We now assume that the length scale on which the applied pressure gradient ∂P/∂z changes is much larger than the size of each triangle. Consequently, the axial flow profile in each triangle and the three adjacent triangles is identical, after appropriate reflections. Therefore, symmetry conditions manifest along each edge of the triangle, ∂w/∂n = 0. Note also the no-slip condition w = 0 on the surfaces of the filaments. Thus we obtain a boundary condition along each boundary, and the problem is fully determined with a unique solution.

The solution to eqn (17) for the axial flow that satisfies the symmetry conditions at θ = 0 and θ = π/6 is given by the truncated series

Once again, this value of N is temporary and distinct from the value of N seen in eqn (16). We use R_thick to denote the rescaled radius of the thick myosin filament. Within the thick and overlap regions, applying w = 0 on the thick filament gives

whilst in the thin region, regularity at r = 0 gives simplyEach of these solutions involves N + 1 coefficients, which we determine numerically by applying the unused symmetry condition on the lower edge, and the no-slip conditions on the thin blue filaments (if present) via a least squares method. Upon computing these coefficients, we have thus obtained a semi-analytic solution for w.

Noting that the height of the triangle was rescaled to 1, giving an area of , and recalling that permeability is defined as the ratio between fluid flux per unit area and pressure gradient, we then calculate a dimensionless permeability in each region , where is the total fluid flux through the triangle. The resulting values for the dimensionless axial permeabilities of real sarcomeres (k^(j)_‖, j = 1, 2, 3) are given in Table 2 in the main text. These permeabilities must be re-dimensionalised to have the appropriate dimensions of length squared by multiplying k^(j)_‖ by H², where is the original, dimensional height of the triangle. To then incorporate these dimensional permeabilities into eqn (16), they must be non-dimensionalised under the same scheme as the sarcomere. Since each k^(j)_‖ refers to axial flow, the appropriate scaling involves a division by L², where L is the dimensional sarcomere length. Therefore, each k^(j)_‖ in Table 2 should overall by multiplied by D₁₀²/3L² to be used in eqn (16).

It should be noted that the least squares method is typically more numerically robust if we rescale the terms of the semi-analytic solution to not increase rapidly with n. For example, since the maximum possible value of r is , we may rewrite eqn (21) as

Here, the value of N required to produce accurate results, as given in Table 2, depends on the packing ratio and varies between the regions of the sarcomere, as well as the particular implementation of the least squares matching. In general, retaining a few dozen terms allows for excellent accuracy, such that the values given in Table 2 are accurate to the given number of digits.

B Calculating the traction parameter β

With the axial permeabilities known, there remains the task of calculating the value of the traction parameter β. This may be done by adapting the axial flow solution from Appendix A in the overlap region (since β exists only in the overlap region), by removing the pressure gradient, and replacing the boundary conditions on the thin blue filament(s) with w = 1. The semi-analytic solution that satisfies w = 0 on the thick filament r = R_thick and the symmetry conditions ∂w/∂θ = 0 on the vertical edge θ = 0 and the hypotenuse θ = π/3 is given by

Again, the coefficients are determined numerically by a least squares method, matching the remaining symmetry condition along the lower edge, as well as the no-slip condition w = 1 on the thin blue filament(s). We then calculate the value of β from the flux, , as

Note that β requires no re-dimensionalisation. The resulting values for β of real sarcomeres are given in Table 2 in the main text.

C Lateral permeabilities

The lateral permeabilities are more complicated to determine, but we can adapt a semi-analytic scheme, originally used by Sangani & Acrivos,⁷¹ who considered a hexagonal array of cylinders of equal radii (i.e. our thick region) as a series of repeating rectangles with quarter-circles in a pair of opposite vertices. We could immediately apply this method to the thick region, but in the other regions, the thin filaments would render the method ineffective. However, by considering additional symmetries not originally noted in ref. 71, we can significantly reduce the size of the domain and the number of numerical boundary conditions, as well as generalise the problem to the presence or absence of thick and/or thin filaments.

To do so, we simultaneously solve for the two solutions resulting from two orthogonal pressure gradients on the same triangles illustrated in Fig. 2F–H, and reproduced for a packing ratio of 5 in Fig. 8, thereby reducing the size of the domain on which the solution is solved by a factor of six compared to the original study.⁷¹ Within the thick region, there is only one boundary condition that we must match numerically, compared to three originally.⁷¹ Finally, we can include thin filaments in any of the three packing ratios of Fig. 2C–H, and even remove the thick filament, and can still solve the problem semi-analytically.

As with the axial permeability calculations, we rescale the triangle from a dimensional height of to a dimensionless height of 1, and we employ a temporary, local polar coordinate system, (r,θ), as in Fig. 8. We begin by identifying the biharmonic equation for the streamfunction ψ, relating to the vorticity ω

∇⁴ψ = −∇²ω = 0, (25)

which produces a velocity

On the surfaces of filaments, we have no slip, giving zero normal derivative, ∂ψ/∂n = 0. We also have zero tangential derivative, ∂ψ/∂t = 0, which we implement as ψ = const, where the value of the constant is set by a flux condition. It will be most convenient to universally specify the Darcy flux (i.e. flux per unit area) to be 1, so the actual flux through the triangle will depend on the direction of the pressure gradient. Specifically, we set the flux across the triangle to be 1 for a horizontal pressure gradient, and for a vertical pressure gradient. We then consider four different directions for the pressure gradient – one horizontal, one vertical, one normal to the hypotenuse, and one parallel to the hypotenuse, with streamfunctions ψ_h, ψ_v, ψ_⊥ and ψ_‖ respectively. As with the axial case, we assume that changes in ∇P occur over length scales much larger than the size of each triangle, giving symmetry conditions along the edges of the triangle, dependent on the direction of the applied pressure gradient. Whenever the pressure gradient is parallel to a side of the triangle, we have ∂ψ/∂t = ω = 0 at that side boundary. When the pressure gradient is normal to a side of the triangle, we have ∂ψ/∂n = ∂ω/∂n = 0 there. We relate the four solutions by

We apply boundary conditions on the hypotenuse for ψ_‖ and ψ_⊥, as described above. Therefore, we have four boundary conditions along the hypotenuse for ψ_v and ψ_h together, leading to four conditions on each and every boundary for the coupled system (ψ_v,ψ_h), which is hence fully determined.

Complete expressions for the streamfunctions, for the three regions of the sarcomere, with corresponding expressions for vorticity and pressure, are provided in the following subsections. Note that, for any given pressure gradient, giving streamfunction ψ in the triangles of Fig. 2F–H, we can use appropriate superpositions of ψ_h and ψ_v (along with additive constants) to determine the extension of the streamfunction to the other triangles, and hence the entire hexagonal cell. Somewhat remarkably, we find that the solution in the entire hexagonal cell is simply ψ extended to a domain of all θ between 0 and 2π.

Finally, the expressions for pressure allow us to determine the dimensionless pressure change Δp over a length l through the hexagonal cell, giving each permeability as

As mentioned in the main text, the value of k^(j)_⊥ is independent of the chosen direction for the applied pressure gradient. The resulting values for the lateral permeabilities of real sarcomeres (k^(j)_⊥, j = 1, 2, 3) are listed in Table 2 in the main text. As with the axial permeabilities, each k^(j)_⊥ must be multiplied by D₁₀²/3 to be re-dimensionalised. Since each k^(j)_⊥ refers to lateral flow, they must then be divided by R², where R is the dimensional sarcomere radius, in order to have the same non-dimensionalisation used for the sarcomere. This means the k^(j)_⊥ in Table 2 must overall be multiplied by D₁₀²/3R² to be used in eqn (16).

C.1 Solutions in the thick and overlap regions. We obtain solutions for the streamfunctions for horizontal and vertical flow:and

These give vorticities

and

We calculate dimensionless pressure analytically by noting that ∂p/∂θ = r∂ω/∂r:

and

C.2 Solutions in the thin region. In the thin region, the solutions are different due to the regularity condition at r = 0. We haveThese give vorticitiesand dimensionless pressures

D Boundary and interface conditions

The solution for the Darcy pressure is the truncated series of eqn (16):where J₀ is the zero order Bessel function of the first kind, λ_n is the nth zero of J₀, and . To identify the values of the coefficients, we apply the boundary and interface conditions (8)–(14)via the following orthogonality condition for Bessel functions:⁹⁶

For each n, we define

where we have evaluated the integral using the definition of J₀(x)and the fact that . We hence obtain the following equation for the coefficients, independently for each value of n,where

Provided with values for the permeabilities, we solve these equations numerically for each n, thereby obtaining the coefficients A^(j)_n and B^(j)_n and hence the pressure P. However, it is generally more numerically robust to write this system as

where L₀ = 0 and L₃ = 1, with corresponding matrix equationwhere the matrix M̃ is given byUnder these rescalings, both P and M̃ contain only terms that are at most [scr O, script letter O] (1) (or more accurately, [scr O, script letter O] (n)), greatly improving numerical accuracy. In practice, we find that the coefficients Ã_n⁽³⁾ and B̃_n⁽¹⁾ decay exponentially with n, whilst the remaining coefficients decay like n^−3/2. We find that, upon reaching n = 50, all coefficients have decayed to less than 2% of that coefficient's peak value across all values of n, and the decay is usually significantly greater than this, around 0.1%. From this, we conclude that a truncation of N = 50 is sufficient to produce accurate values for the pressure field, and hence the velocity fields.

E FEM computations

The full computations (finite element method) were performed dimensionally within COMSOL Multiphysics Version 5.6,⁶⁴ using a stationary study and creeping flow physics interface. We used COMSOL's option for a normal mesh size, with a maximum element size of 47.4 nm, maximum element growth rate of 1.15, curvature factor of 0.6 and resolution of narrow regions of 0.7. However we manually reduced the minimum element size to 6.5 nm, allowing filament cross-sections to be approximated as octagons by the mesh. We applied no-slip conditions on the surfaces of filaments and disks, and symmetry conditions at z = 0 when the M disk was absent. The Z disk and thin filaments had a contraction speed of 1000 nm s⁻¹, and the dynamic viscosity of the fluid was set to that of water at room-temperature, 0.001 Pa s, though both of these are largely irrelevant due to linearity. Symmetry conditions allowed us to consider only a 30° sector of the sarcomere. We approximated a stress-free far-field condition by extending the Z disk and M disk boundaries (whether present or absent), and their corresponding boundary conditions, to twice the sarcomere radius, effectively producing a sarcomere with twice the radius, but with filaments present within only the inner half-radius. We then applied a stress-free condition on the enlarged outer radial boundary.

F Advection–reaction–diffusion

Experiments have shown that diffusion of ATP inside the sarcomere is surprisingly slow, with diffusive coefficient estimated to be D ≈ 2 × 10⁻¹⁵ m² s⁻¹, several orders of magnitude lower than what would be observed in water.⁶⁶ This is likely due to the dense ecology of substrates and enzymes within the sarcoplasm.^52,66 For contraction consistent with human cardiac sarcomeres, this gives a relevant Péclet number of Pe ≈ 125, suggesting that advection has significant impact on ATP distribution.^{9,34,35,42,90–93,97} We model the reaction of ATP as occurring in the overlap region only, and during contraction only. We quantify the reaction using Michaelis–Menten kinetics, with reaction rate V_mC/(K_m + C), where C is the ATP concentration, maintained at 10 mM outside the sarcomere, and Michaelis constants V_m = 1 mM and K_m = 0.01 mM.⁶⁶ We non-dimensionalised the system, scaling all lengths with sarcomere radius R, and scaling ATP concentration with the external ATP concentration of 10 mM. We set boundary conditions C = 1 at r = 1, ∂C/∂r = 0 at r = 0 and ∂C/∂z = 0 at z = 0 and z = L, where L is the dimensionless sarcomere length, and solved for the evolution of the system using an advective grid method that stretches with the sarcomere as it contracts and relaxes sinusoidally, with continuity of both C and substrate flux applied at interfaces between regions, until a periodic state was reached. Comparisons between this periodic state, and the state achieved in the absence of fluid advection, are shown in Fig. 9 (A: no flow; B: advection by flow included) at various stages of contraction (i, ii, iii). We clearly see that the fluid flow significantly reduces the sizes of ATP “deadzones”, i.e. regions of very low ATP concentration, and overall increases the average reaction rate by approximately 26%, representing a significant improvement in sarcomere activity. These preliminary calculations warrant further investigations.

Fig. 9 ATP concentration over the course of a contraction, in a human cardiac sarcomere with no advection (A) and advection by the fluid flow (B). Snapshots of ATP concentration are shown at low (i), intermediate (ii) and high (iii) levels of human cardiac sarcomere contraction. Horizontal axis represents the dimensionless axial coordinate of the sarcomere, and the vertical axis represents the dimensionless radial coordinate. Dashed lines indicate transitions between the key regions of the sarcomere, and reaction occurs in the middle region only.

G Varying permeabilities in the thin region

In mammalian sarcomeres, the presence of the giant protein titin and the disorganisation of the thin filaments within the thin region may give pause for concern regarding the calculated permeabilities k_‖⁽³⁾ and k_⊥⁽³⁾. As discussed in the main text, a semi-analytic calculation of these permeabilities would be unfeasible, however they may be calculated numerically (provided sufficiently accurate physical models of the molecular microstructure are available) or possibly measured experimentally. To demonstrate the potential for the Darcy model to be applied to this new system, we plot in Fig. 10 the total radial efflux leaving the sarcomere, in each of the three regions (A, B, C for regions 1, 2, 3 respectively), for the same sarcomere as in Fig. 7c, consistent with mammalian sarcomeres, with the only changes being variation of k_‖⁽³⁾ and k_⊥⁽³⁾ from their standard values in Table 2 by dimensionless factors ranging between 0.3 (making them comparable to the permeabilities in the overlap region) and 2. These results demonstrate how the Darcy model can be applied to a continuous range of parameters (as well as highlighting the independence of the radial efflux in the thick and overlap regions on the lateral permeability in the thin region, k_⊥⁽³⁾). Importantly, the Darcy model can be applied to any system with specified, even arbitrary permeabilities, without the need to resolve the precise physical structure of the filaments.

Fig. 10 Total dimensionless radial efflux in regions 1 (A), 2 (B) and 3(C) due to variations of k_‖⁽³⁾ and k_⊥⁽³⁾ by factors between 0.3 and 2, compared to their values in Table 2. The porosity Φ⁽³⁾ in the thin region is fixed at its established value used in previous examples in the main text.

Acknowledgements

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement 682754 to EL), and from the European Commission's Erasmus+ programme. We thank Sage Malingen, Richard Tyser, Grae Worster and Marco Vona for useful discussions and advice.

Notes and references

1. B. Alberts, J. Wilson, T. Hunt, R. Heald, A. Johnson, D. Morgan, M. C. Raff, K. Roberts and P. Walter, Molecular biology of the cell, W.W. Norton & Company, 2022.
2. R. W. Lymn and E. W. Taylor, Biochemistry, 1971, 10, 4617–4624.
3. H. Al-Khayat, Glob. Cardiol. Sci. Pract., 2013, 36.
4. H. E. Huxley, Biochim. Biophys. Acta, 1953, 12, 387–394.
5. H. E. Huxley, J. Cell Biol., 1957, 3, 631–648.
6. K. Powers, K. Nishikawa, V. Joumaa and W. Herzog, J. Exp. Biol., 2016, 219, 1311–1316.
7. M. Hou, Biomech. Model. Mechanobiol., 2018, 17, 1797–1810.
8. N. González-Morales, Y. S. Xiao, M. A. Schilling, O. Marescal, K. A. Liao and F. Schöck, eLife, 2019, 8, e50496.
9. T. Shimomura, H. Iwamoto, T. T. Vo Doan, S. Ishiwata, H. Sato and M. Suzuki, Biophys. J., 2016, 111, 1295–1303.
10. OpenStax, CC BY 4.0 https://creativecommons.org/licenses/by/4.0, via Wikimedia Commons, edited.
11. User:Sameerb, Copyrighted free use, via Wikimedia Commons, edited .
12. R. Craig and R. Padrón, Molecular Structure of the Sarcomere, 2004, ch. 7, pp. 129–166.
13. J. S. Wray, P. J. Vibert and C. Cohen, Nature, 1975, 257, 561–564.
14. S. G. Page, J. Cell Biol., 1965, 26, 477–497.
15. L. D. Peachey and A. F. Huxley, J. Cell Biol., 1962, 13, 177–180.
16. P. K. Luther, J. Struct. Biol., 2000, 129, 1–16.
17. P. K. Luther, J. S. Barry and J. M. Squire, J. Mol. Biol., 2002, 315, 9–20.
18. J. P. Schroeter, J. P. Bretaudiere, R. L. Sass and M. A. Goldstein, J. Cell Biol., 1996, 133, 571–583.
19. K. C. Holmes, D. Popp, W. Gebhard and W. Kabsch, Nature, 1990, 347, 44–49.
20. T. D. Pollard, Curr. Opin. Cell Biol., 1990, 2, 33–40.
21. H. Huxley, J. Mol. Biol., 1963, 7, 281–308.
22. J. Trinick and A. Elliott, J. Mol. Biol., 1979, 131, 133–136.
23. P. Knight and J. Trinick, J. Mol. Biol., 1984, 177, 461–482.
24. H. Huxley and W. Brown, J. Mol. Biol., 1967, 30, 383–434.
25. R. W. Kensler and M. Stewart, J. Cell Biol., 1983, 96, 1797–1802.
26. G. G. Knappeis and F. Carlsen, J. Cell Biol., 1968, 38, 202–211.
27. P. Luther and J. Squire, J. Mol. Biol., 1978, 125, 313–324.
28. P. K. Luther and R. A. Crowther, Nature, 1984, 307, 566–568.
29. R. A. Crowther and P. K. Luther, Nature, 1984, 307, 569–570.
30. J. R. Sommer, Z. Naturforsch., C:J. Biosci., 1982, 37, 665–678.
31. A. S. Zot and J. D. Potter, Annu. Rev. Biophys. Biophys. Chem., 1987, 16, 535–559.
32. R. J. Solaro and H. M. Rarick, Circ. Res., 1998, 83, 471–480.
33. A. M. Gordon, E. Homsher and M. Regnier, Physiol. Rev., 2000, 80, 853–924.
34. H. E. D. J. ter Keurs, N. Diao and N. P. Deis, Ann. N. Y. Acad. Sci., 2010, 1188, 165–176.
35. E. K. Rodriguez, W. C. Hunter, M. J. Royce, M. K. Leppo, A. S. Douglas and H. F. Weisman, Am. J. Physiol., 1992, 263, H293–H306.
36. S. Shankar and L. Mahadevan, Nat. Phys., 2024, 20, 1501–1508.
37. A. F. Huxley and R. Niedergerke, Nature, 1954, 173, 971–973.
38. H. Huxley and J. Hanson, Nature, 1954, 173, 973–976.
39. S. G. Page and H. E. Huxley, J. Cell Biol., 1963, 19, 369–390.
40. G. Elliott, J. Lowy and B. Millman, J. Mol. Biol., 1967, 25, 31–45.
41. H. Huxley, A. Faruqi, M. Kress, J. Bordas and M. Koch, J. Mol. Biol., 1982, 158, 637–684.
42. R. J. Stenger and D. Spiro, J. Cell Biol., 1961, 9, 325–351.
43. G. Hoyle, Am. Zool., 1967, 7, 435–449.
44. T. Toumanidou, in Biomechanics of the Spine, ed. F. Galbusera and H.-J. Wilke, Academic Press, 2018, pp. 141–166.
45. J. A. Cass, C. D. Williams, T. C. Irving, E. Lauga, S. Malingen, T. L. Daniel and S. N. Sponberg, Biophys. J., 2021, 120, 4079–4090.
46. S. A. Malingen, K. Hood, E. Lauga, A. Hosoi and T. L. Daniel, Arch. Biochem. Biophys., 2021, 706, 108923.
47. P. P. de Tombe and H. E. ter Keurs, J. Physiol., 1992, 454, 619–642.
48. B. C. W. Tanner, G. P. Farman, T. C. Irving, D. W. Maughan, B. M. Palmer and M. S. Miller, Biophys. J., 2012, 103, 1275–1284.
49. D. M. Swank, W. A. Kronert, S. I. Bernstein and D. W. Maughan, Biophys. J., 2004, 87, 1805–1814.
50. N. Thekkethil, J. Köry, M. Guo, P. S. Stewart, N. A. Hill and X. Luo, Biomech. Model. Mechanobiol., 2024, 23, 1551–1569.
51. J. Najafi, S. Dmitrieff and N. Minc, Proc. Natl. Acad. Sci. U. S. A., 2023, 120(9), e2216839120.
52. M. Hargreaves and L. L. Spriet, Nat. Metab., 2020, 2, 817–828.
53. K. Luby-Phelps, Int. Rev. Cytol., 1999, 189–221.
54. J. C. van der Werff, C. G. de Kruif, C. Blom and J. Mellema, Phys. Rev. A:At., Mol., Opt. Phys., 1989, 39, 795–807.
55. J. van der Werff, C. de Kruif and J. Dhont, Phys. A, 1989, 160, 195–204.
56. R. Verberg, I. M. de Schepper and E. G. Cohen, Phys. Rev. E:Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1997, 55, 3143–3158.
57. F. Kono, S. Kawai, Y. Shimamoto and S. Ishiwata, Sci. Rep., 2020, 10, 16372.
58. T. Washio, S. A. Shintani, H. Higuchi, S. Sugiura and T. Hisada, Sci. Rep., 2019, 9, 9355.
59. D. Hayes, M. Huang and C. R. Zobel, J. Ultrastruct. Res., 1971, 37, 17–30.
60. S. A. Malingen, A. M. Asencio, J. A. Cass, W. Ma, T. C. Irving and T. L. Daniel, J. Exp. Biol., 2020, 223, jeb224188.
61. J. Bear, Dynamics of fluids in porous media, Dover Publications, Inc., 2018.
62. E. Guyon, Physical hydrodynamics, Oxford University Press, 2015.
63. P. G. Saffman, Stud. Appl. Math., 1971, 50, 93–101.
64. COMSOL Multiphysics® 5.6, COMSOL AB, Stockholm, Sweden, 2020.
65. V. A. Selivanov, S. Krause, J. Roca and M. Cascante, Biophys. J., 2007, 92, 3492–3500.
66. A. E. Alekseev, R. Guzun, S. Reyes, C. Pison, U. Schlattner, V. A. Selivanov and M. Cascante, Biochim. Biophys. Acta, 2016, 1860, 2269–2278.
67. T. Kalwarczyk, N. Ziebacz, A. Bielejewska, E. Zaboklicka, K. Koynov, J. Szymanski, A. Wilk, A. Patkowski, J. Gapinski and H. Butt, et al. , Nano Lett., 2011, 11, 2157–2163.
68. K. Kwapiszewska, K. Szczepański, T. Kalwarczyk, B. Michalska, P. Patalas-Krawczyk, J. Szymański, T. Andryszewski, M. Iwan, J. Duszyński and R. Hołyst, J. Phys. Chem. Lett., 2020, 11, 6914–6920.
69. A. Huxley, Reflections on muscle, Liverpool University Press, 1980.
70. H. Wackerhage, U. Hoffmann, D. Essfeld, D. Leyk, K. Mueller and J. Zange, J. Appl. Physiol., 1998, 85, 2140–2145.
71. A. S. Sangani and A. Acrivos, Int. J. Multiphase Flow, 1982, 8, 193–206.
72. H. Huxley and J. Hanson, Biochim. Biophys. Acta, 1957, 23, 229–249.
73. J. Hanson and H. E. Huxley, Biochim. Biophys. Acta, 1957, 23, 250–260.
74. C. C. Gregorio, H. Granzier, H. Sorimachi and S. Labeit, Curr. Opin. Cell Biol., 1999, 11, 18–25.
75. L. Tskhovrebova and J. Trinick, Nat. Rev. Mol. Cell Biol., 2003, 4, 679–689.
76. H. Higuchi, Y. Nakauchi, K. Maruyama and S. Fujime, Biophys. J., 1993, 65, 1906–1915.
77. J. Trinick, P. Knight and A. Whiting, J. Mol. Biol., 1984, 180, 331–356.
78. K. Wang, R. Ramirez-Mitchell and D. Palter, Proc. Natl. Acad. Sci. U. S. A., 1984, 81, 3685–3689.
79. A. Nagy, P. Cacciafesta, L. Grama, A. Kengyel, A. Malnasi-Csizmadia and M. S. Kellermayer, J. Cell Sci., 2004, 117, 5781–5789.
80. R. Yamasaki, M. Berri, Y. Wu, K. Trombitás, M. McNabb, M. S. Z. Kellermayer, C. Witt, D. Labeit, S. Labeit and M. Greaser, et al. , Biophys. J., 2001, 81, 2297–2313.
81. P. Bianco, A. Nagy, A. Kengyel, D. Szatmári, Z. Mártonfalvi, T. Huber and M. S. Z. Kellermayer, Biophys. J., 2007, 93, 2102–2109.
82. M. Kulke, S. Fujita-Becker, E. Rostkova, C. Neagoe, D. Labeit, D. J. Manstein, M. Gautel and W. A. Linke, Circ. Res., 2001, 89, 874–881.
83. W. A. Linke, M. Ivemeyer, S. Labeit, H. Hinssen, J. C. Rüegg and M. Gautel, Biophys. J., 1997, 73, 905–919.
84. W. A. Linke, M. Kulke, H. Li, S. Fujita-Becker, C. Neagoe, D. J. Manstein, M. Gautel and J. M. Fernandez, J. Struct. Biol., 2002, 137, 194–205.
85. K. Trombitas and H. Granzier, Am. J. Physiol., 1997, 273, C662–C670.
86. S. Dutta, C. Tsiros, S. L. Sundar, H. Athar, J. Moore, B. Nelson, M. J. Gage and K. Nishikawa, Sci. Rep., 2018, 8, 14575.
87. A. D. Liversage, D. Holmes, P. J. Knight, L. Tskhovrebova and J. Trinick, J. Mol. Biol., 2001, 305, 401–409.
88. L. A. Tskhovrebova, J. Muscle Res. Cell Motil., 1991, 12, 425–438.
89. J. M. Squire, H. A. Al-khayat, C. Knupp and P. K. Luther, Adv. Protein Chem., 2005, 17–87.
90. L. Burbaum, J. Schneider, S. Scholze, R. T. Böttcher, W. Baumeister, P. Schwille, J. M. Plitzko and M. Jasnin, Nat. Commun., 2021, 12, 4086.
91. D. H. Moore and H. Ruska, J. Cell Biol., 1957, 3, 261–268.
92. W. A. Linke, V. I. Popov and G. H. Pollack, Biophys. J., 1994, 67, 782–792.
93. T. Oda and H. Yanagisawa, Commun. Biol., 2020, 3, 585.
94. https://github.com/John-Severn/Fluid-mechanics-of-sarcomeres-as-porous-media—MATLAB-code .
95. E. M. Sparrow and A. L. Loeffler, AIChE J., 1959, 5, 325–330.
96. M. Abramowitz and I. Stegun, Handbook of mathematical functions, Dover Publishing Inc., New York, 1970.
97. K. C. Vinnakota and J. B. Bassingthwaighte, Am. J. Physiol.: Heart Circ. Physiol., 2004, 286, H1742–H1749.

---