We model the cell culture media as a Newtonian viscous fluid, with viscosity μ and density ρ. We neglect the effect of cells and scaffold deformation, and assume that the scaffold is a homogeneous rigid porous material. We note that the analysis may be readily extended to include two or more distinct porous layers. The heights of the inlet reservoirs, as well as the valve configurations can be varied to control fluid flow through the bioreactor system. Throughout this paper we use tube system to mean the combined channel, resistor tubing, valve and valve tubing. Flow is governed by the incompressible Navier-Stokes equations in the upper and lower tube systems, and by Darcy’s law in the scaffold. We impose no-slip and no-flux conditions on the impermeable rigid boundaries of the channels, tubes and valves. To model valves opening and closing, we prescribe time-dependent wall geometry. At the channel-scaffold interface we impose continuity of normal flux, continuity of normal stress, and the Beavers-Joseph slip condition (Beavers and Joseph, 1967).

In this paper, we present a simple model that captures the flow in the numerous different regions of the bioreactor system. In particular, we assume that the aspect ratio and flow rate are everywhere small enough that we can adopt the lubrication approximation, and thus we neglect inertia in the system. We note that for junctions in which there are sharp changes in cross-sectional area, the lubrication approximation breaks down. However, as discussed at the end of Section 2.2.4, we can systematically justify imposing continuity of flux and pressure conditions at these junctions.

The (equal) centreline arclengths of the upper and lower tube systems are denoted l*. The upper and lower tube systems are labelled by i = u, ℓ respectively, throughout the following. We denote by h u * ( h ℓ * ) the height of the inlet of the upper (lower) tube system above that of its outlet, so that gravity generates a head of ρ g h u * ( ρ g h ℓ * ) across the upper (lower) tube system; see Figure 2. We adopt an orthogonal curvilinear coordinate system (x*, y*, z*) in each of the two tube systems, with the z* − axis running along the axial centreline from inlet to outlet, and x* and y* being cross-sectional coordinates.

The upper and lower tube systems have four regions differing in their cross section, shown in Figure 3: channels, with fixed, rectangular cross-sections of width 2c* and height 2 b c * ; resistor tubing, with fixed, circular cross-sections of radius b r * ; valve tubing, with fixed, circular cross-sections of radius b v * , and valves, with time-dependent cross-sections, idealised for modelling purposes to be elliptical.

We denote by z j * for j = 0, … , 9 the z-coordinate of the junctions between the components, as illustrated in Figure 3. Valve tubing occupies z * ∈ ( z 0 * , z 1 * ) ∪ ( z 8 * , z 9 * ) and each piece has length l t * , valves occupy z * ∈ ( z 1 * , z 2 * ) ∪ ( z 7 * , z 8 * ) and each valve has length l v * , resistor tubing occupies z * ∈ ( z 2 * , z 3 * ) ∪ ( z 6 * , z 7 * ) and each piece has length l r * , and channels occupy z * ∈ ( z 3 * , z 6 * ) and have length l c * . Thus, the cross-sectional areas in the valve tubing, resistor tubing, and channel regions of the bioreactor are given by A i ∗ = π b v * 2 for z * ∈ z 0 * , z 1 * ∪ z 8 * , z 9 * , π b r * 2 , for z * ∈ z 2 * , z 3 * ∪ z 6 * , z 7 * , 4 b c * c * for z * ∈ z 3 * , z 6 * . (1)

The elliptical cross-sections of the valves have semi-major axes a i j * ( z * , t * ) and semi-minor axes b i j * ( z * , t * ) , where j = 1, 7, and t* denotes time. We parameterise the valve walls as ( x * , y * , z * ) = ( a i j * ⁡ cos ⁡ ϕ , b i j * ⁡ sin ⁡ ϕ , z * ) , for ϕ ∈ 0,2 π , z * ∈ ( z j * , z j + 1 * ) , where j = 1, 7, so that the cross-sectional valve areas are A i * = π a i j * b i j * . The semi-minor and semi-major axes are determined by prescribing the area and circumference of the valve cross-section at each axial location and time. At each point in time, the cross-sectional area of each valve located in z 1 * < z * < z 2 * or z 7 * < z * < z 8 * is prescribed to take one of the forms A i * = δ π b v * 2 if the valve is closed, A i * = π b v * 2 if the valve is open, and A i ∗ = π b v ∗ 2 1 − 1 − δ 4 1 − cos 2 π z * − z j ∗ z j + 1 * − z j ∗ 1 − cos 2 π t * T , j = 1,7 (2) if the valve is moving between being closed and open, with period T for one full opening and closing cycle. The small parameter δ controls how completely the valve can close. (For numerical convenience, the valves are prevented from closing completely.) Note that as the exact functional form of the time dependency in (2) is unknown, for modelling purposes we have chosen an illustrative cosine form for simplicity. The valve is most heavily constricted at its centre, where A i * can fluctuate between π b v * 2 (fully open) and δ π b v * 2 (closed). The ellipses’ circumferences are given by 4 a i j * E ( e i j * ) , where e i j * = 1 − b i j * 2 / a i j * 2 is the eccentricity, and E is the complete elliptic integral of the second kind. The circumferences are fixed to be constant in time, and equal to the circumference of the valve tubing.

The velocities of the walls in the bioreactor are U i ∗ = ∂ / ∂ t * ( a i j * ⁡ cos ⁡ ϕ , b i j * ⁡ sin ⁡ ϕ , z * ) for z * ∈ ( z j * , z j + 1 * ) for j = 1, 7, and U i * = 0 for z * ∈ ( z 0 * , z 1 * ) ∪ ( z 2 * , z 7 * ) ∪ ( z 8 * , z 9 * ) .

The scaffold has centreline arclength, l s * , and rectangular cross-section with height h s * and width 2 b s * . It occupies z 4 * < z * < z 5 * .

The cross-sections of the upper tube system, lower tube system, and scaffold are labelled as Ω _u , Ω _ℓ and Ω _s , respectively, with boundaries ∂Ω _u , ∂Ω _ℓ , ∂Ω _s , and interfaces Γ _i = ∂Ω _i ∩ ∂Ω _s between scaffold and channel cross-sections.

Typical values for h i * , l * , l v * , l c * , l t * , c * , b c * , b v * , b r * and b s * are given in Table 1. The values reflect those of the bioreactor in Shepherd et al. (2018), and the wider literature on the permeability of collagen scaffolds (Mohee et al., 2019). The model is valid for a wider range of parameters than those in Shepherd et al. (2018) (discussed in Section 2.2.4); thus, to explore bioreactor design considerations, we provide illustrative results both for the parameter values of Shepherd et al. and for a wider range of parameters.

We consider flow of a viscous Newtonian fluid of viscosity μ and density ρ. Neglecting the effects of curvature, as the ratio of transverse lengthscale to radius of curvature of the tube centreline is small, the flow in the upper and lower tube systems, Ω i × ( z 0 * , z 9 * ) , is governed by the incompressible Navier-Stokes equations ρ ∂ u i ∗ ∂ t * + u i ∗ ⋅ ∇ u i ∗ = − ∇ p i ∗ + μ ∇ 2 u i ∗ , ∇ ⋅ u i ∗ = 0 , (3) and the flow in the scaffold, Ω s × ( z 4 * , z 5 * ) , is governed by Darcy’s law with the continuity equation u s ∗ = − K μ ∇ p s ∗ , ∇ ⋅ u s ∗ = 0 . (4)

Here, u i ∗ are the velocities in the upper and lower tube systems, with components ( u i ∗ , v i ∗ , w i ∗ ) in the (x, y, z) directions, and the reduced pressures are p i ∗ = − p atm + ρ g Y * + P i ∗ in terms of Y*, the elevation in the lab frame, and the absolute pressure P i ∗ . The velocity in the scaffold is u s ∗ , with components ( u s ∗ , v s ∗ , w s ∗ ) in the (x, y, z) directions, the reduced pressure in the scaffold is p s ∗ = − p atm + ρ g Y * + P s ∗ , and the scaffold permeability is K. Note that transverse flow through the scaffold refers to vertical flow, v s * , through the scaffold, and axial flow in the scaffold refers to horizontal flow, w s ∗ .

At the solid walls of the tube systems we impose the usual no-slip and no-flux conditions given by u i ∗ = U i ∗ on ∂ Ω i . (5)

We prescribe normal stress and no tangential velocity at the inlet and outlet of the upper and lower tube systems, viz. u u * ⋅ t = 0 , n ⋅ σ u * ⋅ n = − ρ g h u * at z * = 0 , (6a) u u * ⋅ t = 0 , n ⋅ σ u * ⋅ n = 0 at z * = l * , (6b) u ℓ * ⋅ t = 0 , n ⋅ σ ℓ * ⋅ n = − ρ g h ℓ * at z * = 0 , (6c) u ℓ * ⋅ t = 0 , n ⋅ σ ℓ * ⋅ n = 0 at z * = l * , (6d) where t is any tangent to the cross-section, σ i ∗ = − p i ∗ I + μ ( ∇ u i ∗ + ( ∇ u i ∗ ) T ) / 2 is the stress tensor in the pipe, and g is the gravitational acceleration. On the impermeable scaffold walls we have the no-flux condition, so that u s ∗ ⋅ n = 0 on ∂ Ω s \ Γ u ∪ Γ ℓ . (7)

At the interfaces Γ _u and Γ _ℓ between scaffold and channels, we prescribe continuity of normal flux, continuity of normal stress, and the Beavers-Joseph-Saffman condition: u i ∗ ⋅ n = u s ∗ ⋅ n on Γ i , (8a) n ⋅ σ i ∗ ⋅ n = − p s ∗ on Γ i , (8b) n ⋅ ∇ u i ∗ + ∇ u i ∗ T ⋅ t = α K u i ∗ ⋅ t on Γ i , (8c) where α is the Beaver-Joseph slip coefficient, and t is any tangent to the interface. The system is closed by imposing the initial fluid velocities u u * and u ℓ * at t* = 0 everywhere except in the scaffold region, but we shall consider a quasi-steady reduction that does not require the imposition of initial conditions (the original initial conditions being satisfied in a short temporal boundary layer).

We exploit the slender geometry of the bioreactor to derive systematically a reduced model, where we consider the lubrication regime in which the aspect ratio and reduced Reynolds number are small in each tubing section. For the purposes of non-dimensionalisation, we introduce the aspect ratio ϵ = b r * / l * ≪ 1 , which is the ratio of the smallest transverse lengthscale the largest axial lengthscale. We non-dimensionalise the governing Eq. 3, and boundary conditions (5), (Section 2.2.3) in the upper and lower tube systems using the following scalings x * , y * , z * = l * ϵ x , ϵ y , z , p i ∗ = ρ g h u * p i , u i ∗ , v i ∗ , w i ∗ = U ϵ u i , ϵ v i , w i , t * = l * U t . (9)

The reduced pressures are non-dimensionalised with respect to the pressure head across the upper tube system. We choose the velocity scaling U to reflect a balance between viscous dissipation and the axial pressure gradient, so that U = ϵ 2 ρ g h u * l * μ . (10)

The governing Eq. 4 and boundary conditions (7) and (Section 2.2.3) in the scaffold are non-dimensionalised using the following scalings x * , y * , z * = l * ϵ x , ϵ y , z , p s ∗ = ρ g h u * p s , u s ∗ , v s ∗ , w s ∗ = U s u s , v s , ϵ w s . (11)

The scale for the reduced pressure is the same as in the tube systems, motivated by the continuity of stress condition (8b) at the interfaces Γ _u,ℓ . Darcy’s law (4) sets the cross-sectional velocity scaling as U s = K ρ g h u * μ h s ∗ . (12)

In addition to ϵ and the various dimensionless lengths in the bioreactor (see Tables 2, 3 for lists of dimensionless lengths and coordinates, respectively), the resulting dimensionless governing equations are characterised by the dimensionless parameters R = ϵ 2 ρ l * U μ , A = K α l * , B = K ϵ 4 l * 2 S = l * U T . (13)

The reduced Reynolds number, R , represents the ratio of inertia to viscous effects based on the transverse velocity and lengthscale in the tube systems. The dimensionless Beavers-Joseph slip coefficient, A , is the ratio of the slip length K / α to the system length l*. The dimensionless permeability, B , measures the ratio of the transverse scaffold velocity, U _s , to the size of transverse tubing velocities, ϵU. The Strouhal number, S , is the ratio of the timescale of flow driven by the upper reservoir, and the timescale of valve motion.

We work in the physically relevant regime in which the reduced Reynolds number R = O ( 1 ) , A = O ( 1 ) ε , and S = O ( 1 ) as ϵ → 0, and retain only leading order terms, so that errors are O (ϵ ²). The material parameters and valve timescale must be chosen to ensure that we are in the correct regime.

We reduce the system to a network of pressure—cross-sectional area—flux relations at leading order as follows. Integrating the continuity Eq. 3 over the cross-sections Ω _i , imposing the kinematic boundary conditions (5) and the interface flux condition (8a), and using Reynold’s transport theorem, we find that the relationship between the axial flux and cross-sectional areas in the upper and lower tube systems is given, for 0 < z < 1, by ∂ A i ∂ t + ∂ Q i ∂ z = J i . (14)

Here, the cross-sectional area, axial flux and net flux out of the adjoining scaffold are given by A i z , t = ∫ ∫ Ω i d x d y , Q i z , t = ∫ ∫ Ω i w i d x d y , (15a) J i z , t = − B ∫ Γ i ∂ p s ∂ n d s for z ∈ z 4 , z 5 , 0 , otherwise, (15b) with ∂/∂n being the derivative in the direction of the outward unit normal to the scaffold cross-section, so that positive J _i indicates flux from the scaffold into the tube system.

The pressure-flux relations in the upper and lower tube systems are determined as follows. The transverse components of the momentum Eq. 3 imply that the leading-order pressure in both tube systems is independent of the transverse coordinates x and y.

To determine Q _i we must solve the leading-order axial momentum Eq. 3, which is given by ∂ 2 w i ∂ x 2 + ∂ 2 w i ∂ y 2 = ∂ p i ∂ z in Ω i . (16)

We solve (16) subject to no slip on the impermeable walls of the tube systems. We impose the leading-order version of the Beavers-Joseph slip condition (8c) at the scaffold-channel interface when A is small, which again reduces to the no slip condition. Hence we impose w i = 0 on ∂ Ω i . (17)

Substitution of Darcy’s law into the continuity equation in (4) and retaining leading order terms gives at leading order Laplace’s equation in the scaffold cross-section, namely ∂ 2 p s ∂ x 2 + ∂ 2 p s ∂ y 2 = 0 in Ω s . (18)

We solve (18) subject to the leading-order versions of the boundary conditions of 1) continuity of normal stress at the scaffold-channel interfaces (8b), and 2) no flux at the impermeable walls of the scaffold, which are given by p s = p u on Γ u , p s = p ℓ on Γ u , ∂ p s ∂ n = 0 on ∂ Ω s \ Γ u ∪ Γ ℓ . (19)

We solve Poisson’s Eq. 16 subject to the no-slip condition (17) to determine w _i in terms of the pressure gradient, and use (15a) to obtain the following flux-pressure relation Q i = − C i ∂ p i ∂ z , (20) where C _i are the tube system conductivities, determined by the various cross-sectional geometries, in the valve tubes, valves, resistor tubes, and channels, respectively: C i z , t = b v 4 π 8 for z ∈ z 0 , z 1 ∪ z 8 , z 9 , π a i j 4 1 − e i j 2 3 / 2 4 e i j 4 2 − e i j 2 , for z ∈ z j , z j + 1 , j = 1,7 , b r 4 π 8 for z ∈ z 2 , z 3 ∪ z 6 , z 7 , ∑ n , m = 0 ∞ 16 b c c λ n 2 ν m 2 λ n 2 + ν m 2 for z ∈ z 3 , z 6 (21) where λ _n = (n + 1/2)π/b _c , ν _m = (m + 1/2)π/c, z i = z i * / l * for i = 0, … , 9, b v = b v * / b r * , b r = b r * / b r * , b c = b c * / b r * , c = c * / b r * and a i j = a i j * / b r * . Figure 4 shows a plot of conductivities when the valves are closed.

To determine the scaffold pressure, and hence the net flux out of the scaffold into the adjoining channel J _i , we solve (18) subject to (19). Since we have assumed that the scaffold cross section is a rectangular domain with Dirichlet and Neumann boundary conditions holding on horizontal and vertical sides respectively, the solution is simply p s y , z , t = p u z , t − p ℓ z , t y h s + p ℓ z , t , (22) and hence J u = − J ℓ = − 2 c B h s p u − p ℓ for z 4 < z < z 5 0 otherwise. (23)

The dimensionless cross-sectional areas (1) in the channels, valve tubing, and resistor tubing are the same with stars dropped. The dimensionless cross-sectional area (2) of the moving valves is A i = π b v 2 1 − 1 − δ 4 1 − cos 2 π z − z j z j + 1 − z j 1 − cos 2 π S t , j = 1,7 . (24)

The cross-sectional area along the centreline, when the valves are closed, is plotted in Figure 4.

To summarise, the pressure—cross-sectional area—flux relations are given by Eqs 14, 20, with tube conductivities (21), flux (23) from scaffold to channels, valve cross-sectional areas (24) and the dimensionless version of tube cross-sectional areas (1). This system of second-order linear differential equations is solved subject to the specification of inlet and outlet tube system pressures (Section 2.2.3), which in dimensionless form are p u 0 , t = 1 , p u 1 , t = 0 , p ℓ 0 , t = h ℓ , p ℓ 1 , t = 0 , (25) and continuity of pressures p _i and fluxes Q _i (i = u, ℓ) at the interfaces between the valves, resistor tubing and channels, as illustrated in Figure 4. When there is a jump discontinuity in the cross-sectional area, these conditions can be derived systematically using a matched asymptotic analysis, as follows. We rescale the governing Eq. 3 such that the axial and radial lengthscales are both O ( ϵ ) , and the axial and radial velocity scales are both O ( 1 ) , near these interfaces. The leading order momentum equations then show that there are no leading order changes in pressure across the interfaces, and the mass equation implies no leading order changes in flux across the interfaces.

Using the parameters given in Table 1, the lower inlet pressure is h _ℓ = 0.5. Having computed the pressures in the tube systems, the flux through the scaffold can be computed from (23).

We analytically solve the system of second order differential equations for the velocities and pressures in each of the 18 regions, with each pair of velocities and pressures having two degrees of freedom. By imposing the conditions at the inlets, outlets, and junctions between regions, as described above, we obtain a system of 36 linear algebraic equations relating the 36 degrees of freedom. These are readily solved numerically in MATLAB using the backslash command.

Within the bioreactor, the aim is to ensure the MKs are subject to sufficient levels of shear stress to ensure effective platelet production, while the shear stress stays below the level at which platelet activation and aggregation occurs. Although the precise values of these bounds are unknown, in vivo the MKs and proplatelets in the bone marrow experience shear stress in the range 0.10−0.70 Pa (Pries et al., 1992; Mazo and von Andrian, 1999). These values provide motivation for desired levels of shear stress in vitro [although it has been shown that successful platelet production in vitro can occur at substantially higher shear rates (Dunois-Lardé et al., 2009)].

We estimate the shear stress distributions from predictions of the Darcy velocity through the porous scaffold following a model that has previously been used in the literature to relate flux to shear stress in porous scaffolds that are modelled using Darcy’s law Whittaker et al. (2009). In this model the pores are assumed to be cylindrical tubes and the geometry of the scaffold enters only via the average pore size and the scaffold tortuosity. The model does not account for scaffold deformability and elasticity, similarly to Darcy’s law.

We idealise the pores as cylinders of diameter d through which there is Poiseuille flow. Then the dimensional axial pore velocity is u pore ( r ) = 2 U pore ( 1 − 2 r / d 2 ) where r is the local radial coordinate in a pore and U _pore is the average pore fluid velocity. The dimensional shear stress is then Σ : = μ ∂ u pore / ∂ r r = d / 2 = 8 μ U pore / d . The dimensional Darcy velocity is related to U _pore via | u s * | = ϕ U pore / τ , where ϕ is scaffold porosity and τ is scaffold tortuosity (see Table 1). At leading order, | u s * | = | v s * | because from Eqs 4, 22, u s * = 0 , while w s * = O ( ϵ ) is lower order than v s * = O ( 1 ) . We therefore estimate the dimensional scaffold shear stress to be Σ = 8 μ τ U s ϕ d v s . (26)

We consider various norms of the stress as part of our analysis, and define the following shear stress metrics: σ = 1 l s ∫ z 4 z 5 Σ d z , σ ̄ = 1 T 1 ∫ 0 T 1 σ d t , and max σ = max t ∈ 0 , T 1 σ , (27) which are the average shear stress, taken over the length of the scaffold; the time-averaged shear stress over the interval t ∈ (0, T ₁) during which at least one valve is moving; and the maximum instantaneous shear stress, over the interval t ∈ (0, T ₁).

In addition to shear stress metrics, we also report the following fluxes: Q s t = − 2 b s ∫ z 4 z 5 v s z , t d z , min Q s = min t ∈ 0 , T 1 Q s t , total Q s = ∫ 0 T 1 Q s t d t , Q u , ℓ m = Q u , ℓ 1 2 , t , and Q out = Q ℓ z 9 , t , (28) which are the vertical flux through the scaffold, defined as positive when flow is from the upper to lower channels; the minimum instantaneous flux through the scaffold, over a period t ∈ (0, T ₁) during which valves are moving; the total flux through the scaffold, over a period t ∈ (0, T ₁); the axial fluxes halfway along each channel, and the flux out of the lower outlet. Dimensional versions of the above quantities are indicated with stars.

We compute the average scaffold shear stress, σ, and the scaffold flux, Q s * , for all the possible static configurations involving two, three, or four open valves, with the remaining valves closed. In Figures 5Ai,ii, we plot the average shear stress and the scaffold flux. We order the configurations by decreasing scaffold flux. As expected, configurations with the upper inlet valve open and a lower valve open give a non-negative scaffold flux, i.e., flow from the upper to the lower channel, because the upper inlet pressure is double the lower inlet pressure. Back-flow, i.e., flow from the lower to upper channel, occurs when the upper inlet valve is closed, and the upper outlet and lower inlet valves are both open. Figure 5Ai shows that shear stress is maximised by closing the upper outlet and lower inlet valves, while leaving the upper inlet and lower outlet valves open. As seen in Figure 5Aii, this same configuration also maximises scaffold flux Q s * , which promotes advective nutrient transport into the scaffold.

In Figure 5B we consider the shear-maximising configuration and explore how opening either one, or both, of the shut valves modulates the scaffold shear stress. Note that opening either one, or both, of the closed valves results in scenarios that have downward scaffold flux and open lower outlet valves, which we anticipate will promote platelet collection (see below). Recall that the degree of valve opening is controlled by the dimensionless parameter δ in (2), which varies between 10^–4 (almost fully shut) and 1 (fully open). We see that either opening the upper outlet valve, or both the upper outlet and lower inlet valves, by as little as 3%, at δ = 0.03, reduces shear stress by 44% and 31%, respectively.

We compute the average shear stress σ as we vary the scaffold permeability, K, for all the valve configurations shown in Figures 5Aii that have positive scaffold flux. For all of these configurations, Figure 5D shows that increasing K from 10^–15 m² to 10^–11 m² increases the shear stress σ by approximately three orders of magnitude, but increasing K further has minimal effect on the shear stress, as the limiting factor is no longer permeability but boundary pressures and tube system geometry. Additionally, within the range of permeabilities K = 10^–12 m² to 10^–9 m² for the scaffold used by Shepherd et al. (2018) (see Mohee et al. ,2019), the shear stress for the two shear-maximising configurations fall within the physiological target range of 0.1–0.7 Pa, as shown by the dark and light blue lines in Figure 5D. In the bioreactor operation of Shepherd et al., the configurations with all valves open (grey line), and with all but the lower inlet valve open (dark red line) were used. Figure 5D shows that the latter of these reaches the target shear stress in the permeability range K = 10^–11 m² to K = 10^–9 m², but the former does not.

For scaffold permeabilities between 10^–15 m² and 10^–12 m², opening the lower inlet valve, in addition to the upper inlet and lower outlet valves, significantly decreases scaffold shear stress, as shown by the blue dashed line in Figure 5C being visibly lower than the dark blue line. This is because, for small enough scaffold permeabilities, flow in the lower tube system dominates flow through the scaffold, in that it substantially increases pressure in the lower tube, relative to its value when the lower inlet is closed. In contrast, for large scaffold permeabilities, opening the lower inlet valve barely increases the pressure in the lower tube system, thus the shear stress is only slightly less than when the lower inlet valve is closed.

The average scaffold shear stress is also affected by varying the tube system length l*. In Figure 5D, we plot the change in scaffold shear stress as we increase the overall tube length by lengthening each of the regions: valve tubing, channels, resistor tubing, and the scaffold with its adjoining channel regions. We only present results for the shear-maximising valve configuration depicted in the bottom left inset, as we obtain the same qualitative results as all the valve configurations shown in Figures 5Aii that have non-negligible scaffold flux (i.e., all but those in the sixth column). Varying l* by lengthening one particular region affects the model in three places. First, increasing l* acts to decrease velocities in the tube systems, reflected in the velocity scaling (10). Second, increasing l* increases the dimensionless permeability B in Eq. 13, representing an increase in scaffold velocity relative to tube system velocities. Lastly, increasing l* by lengthening valve tubing or channels increases the proportion of the system with relatively wide tubing, decreasing resistance to flow and acting to increase velocities in the tube systems. Conversely, adding resistor tubing will increase the proportion of the system with relatively narrow tubing, increasing resistance to flow in the tube systems. To determine the net effect of these competing factors, we must solve the system.

Figure 5Di shows that doubling l* from 1.4 to 2.8 m by increasing either the valve tubing or channels reduces the scaffold shear stress by less than 1%. This is because the increase in axial length is approximately balanced out by the increase in proportion of relatively wide tubing in the tube systems, so that flow rates in the upper and lower channels change by less than 1%. Thus there is little change in the flux through the scaffold, so that the scaffold velocity, and therefore shear stress, are virtually unchanged. In contrast, increasing the resistor tubing from 1.4 to 2.8 m reduces the shear stress by 81%. This is because the axial lengthscale increases and the proportion of the system with relatively narrow tubing increases, so that the velocity in the upper channel decreases, as does the velocity in the scaffold.

Figure 5Dii shows that doubling the scaffold length (and also the section of channel adjoining the scaffold) reduces shear stress by 40%. The flow rates in the upper and lower channels change by less than 3%, as the increase in proportion of relatively wide tubing in the tube systems again approximately balances the increase in axial length. Thus the volume of media that passes from upper to lower channels at each instant in time stays roughly constant. As the length of the scaffold increases, the velocity in the scaffold therefore decreases, as does the shear stress.

After platelets break off from their parent cell, they should be washed out of the system quickly, so as to minimise risk of damage or activation. This requires sufficient flux down the scaffold and out of the lower outlet. In Figure 5E, we measure the flux Q out * out of the lower outlet, for different valve configurations. The greatest flux is obtained when the upper inlet valve and lower valves are open, as shown in the left column of Figure 5E. The same figure shows that the configuration with the upper inlet and lower outlet valves open has only slightly lower Q out * , and as discussed above additionally promotes scaffold shear stress and advective nutrient transport. If, as in the experiments of Shepherd et al. (2018), the valves are left open for one in 20 s, over 20 h, then with the upper inlet and lower outlet valves open approximately 35 ml of product solution will be collected at the lower outlet. Should we need more flux for platelet collection, we can use flow generated with only lower valves open to wash platelets out, for the following reasons. As seen in Figure 5E, a reasonably high flux of 5.2 μl/s is achieved when only lower valves are open, and from Figure 5Aii we additionally see that back-flow in the scaffold is avoided. By keeping the upper inlet closed, we also avoid unnecessarily using nutrient-rich media from the upper reservoir to wash platelets out. Note that the duration of opening lower valves should still be minimised as far as possible, to avoid excessive dilution of the product.

To allow diffusive nutrient transport, where there is a relatively high flux along the upper channel, while avoiding excessive collection volume at the lower outlet (i.e., dilute product), the configuration with the upper valves open and lower valves closed can be used. As seen in Figure 5E, there is no flow out of the lower outlet when the lower valves are closed. In this configuration advection in the upper tube system maintains the concentration at the inlet value throughout the upper tube system (at leading order in the high Peclet number limit), allowing diffusion across the scaffold (where there is no flow) to deliver the oxygen or nutrient to the cells, on a timescale much longer than the advective timescale.

Broadly, Figures 6Ai,Bii show that while a valve is closing (opening), the resistance of the valve increases (decreases), and flux away from (towards) the valve is enhanced. Flux down the scaffold is therefore enhanced (diminished) during most of the time that an upper valve is closing (opening) or a lower valve is opening (closing), as illustrated by Figures 6Aiii,Biii. Figures a i), a iii), b ii) and b iii) all contain negative fluxes. This indicates that back-flow, where channel flow is from outlet to inlet, or scaffold flow is from bottom to top, can be induced when outlets are closing or inlets are opening. As seen in Figures 6Aii,Bi, the movement of upper (lower) valves has minimal effect on the lower (upper) channel flux, because contribution from the change in scaffold flux is small relative to the lower (upper) channel flux induced by the prescribed pressure drop.

When the valves move on a slow timescale with S = 1 , Figures 6C,D show that all of the fluxes vary monotonically either side of a local extreme that approximately coincides with the valve being closed. When a valve closes, the flux in its channel drops, and then increases when the valve reopens. The size of the drop depends on the pressure boundary conditions. For instance, when the lower inlet valve closes, there is just a small drop in lower channel flux Q ℓ m , as shown by the dark blue line in Figure d ii). This is because the prescribed pressure drop from upper inlet to lower outlet is one, i.e., high enough to compete with the prescribed pressure drop of one from upper inlet to upper outlet. On the other hand, the light blue line in Figure d ii) shows that closing the lower outlet valve, with lower inlet valve open, reduces the lower channel flux Q ℓ m nearly to zero, as the prescribed pressure drop from upper inlet to lower inlet is only a half, i.e., less than the drop of one from upper inlet to lower outlet.

The above principles seen from moving individual valves can be applied to determine what valve combinations should be used in transitioning between static configurations. There is freedom to choose the ordering of valve movements, and the extent to which different valves move in sync with each other, by choosing when each valve starts opening or closing. We first discuss valve ordering for fast valves, and then examine slow valves.

To illustrate valve ordering for fast valves, suppose we wish to maximise shear stress in the scaffold when transitioning from the configuration with only upper valves open (promoting diffusive nutrient transport) to the configuration with only lower valves open (promoting platelet collection). Then, to take advantage of the enhanced scaffold fluxes shown in Figure 6Aiii while upper valves are closing and b iii) while lower valves are opening, the upper valves should be closed simultaneously with the lower valves opening. Closing both of the upper valves simultaneously additionally reduces back-flow in the upper channel, as the increase in upper channel flux caused by closing the inlet valve (red line Figure 6Ai) counteracts the back-flow in the upper channel caused by closing the upper outlet valve (orange line Figure 6Ai).

If the aim is relatively simple, then the principles explained above may be sufficient to pinpoint an appropriate set of valve configurations. More generally, when moving between any two static configurations, we may wish to optimise some function of Q u m , Q ℓ m , Q _s and σ. As our model is computationally inexpensive, a large subset of the space of dynamic configurations that transition between any two static configurations may be simulated in seconds, enabling optimal configurations to be explored. We will illustrate with how this might be done with two examples.

Suppose we move from the static configuration promoting shear stress, advective nutrient transport, and platelet collection, to the static configuration promotive diffusive nutrient transport, as shown in the upper right panel of Figure 7A. The upper outlet valve must open and the lower outlet valve must close. The two valves move with a relative delay time of θ ∈ [ − 0.5, 0.5], where θ = 0 corresponds to valves in-sync, θ = 0.5 is the upper outlet valve opening completely, then the lower outlet valve closing, and θ = −0.5 is the lower outlet valve closing completely, then the upper outlet valve opening. In Figure 7A we predict the proportion of time above an illustrative threshold shear stress of 6.5 × 10^–3 Pa, the total scaffold flux total (Q _s ), and the minimum instantaneous scaffold flux min (Q _s ), while varying the synchronisation of valves. The total time is considered to be the time during which at least one valve is moving. From Figure 7, we see we can for example choose θ to maximise the proportion of time that the shear stress is above the chosen threshold, under the constraint that min (Q _s ) is above some value, so that instantaneous back-flow is sufficiently weak, and the constraint that total (Q _s ) is below some chosen value, so that collection volume is sufficiently low.

In moving from the lower outlet valve closing first (θ = −0.5), to the upper outlet valve opening first (θ = 0.5), we see from Figure 7Ai that the proportion of time spent above the threshold shear stress slightly increases in θ ∈ (−0.5, −0.44), then decreases in θ ∈ (−0.44, 0), before increasing in θ ∈ (0, 0.5). This is because the earlier (later) the lower outlet valve is closed, the less (more) time there is during which flow can pass from an upper valve, to a lower valve, exerting shear stress on the scaffold. Figure 7Aii shows that the total scaffold flux is smallest when the upper and lower outlet valves are in-sync, at θ = 0. Opening the upper outlet valve and closing the lower outlet valve promotes back-flow through the scaffold, in competition with the externally imposed pressure drop from upper inlet to lower outlet. Figure 7Aiii shows that the highest degree of instantaneous back-flow, captured via min (Q _s ), occurs for θ = −0.5, where the upper outlet valve only starts opening after the lower outlet has fully closed. Away from θ = −0.5, there is little variation in min (Q _s ).

In Figure 7B, we show the proportion of time above the threshold shear, the total scaffold flux, and minimum instantaneous scaffold flux when moving from the shear-promoting static configuration, to the collection-promoting static configuration depicted on the right panel. Now θ = −0.5 is opening the lower inlet valve first, and θ = 0.5 is closing the upper inlet valve first. As predicted earlier, closing the upper inlet valve and opening the lower inlet valve pushes fluid down the scaffold, and from Figure 7Biii we see that back-flow never occurs. Closing the upper inlet valve simultaneously with opening the lower inlet valve (θ = 0) gives a high proportion of time above the threshold shear, but a lower total scaffold flux than first opening the lower inlet (θ = −0.5), as seen in Figures 7Bi,ii, respectively. The total flux is lower because the time during which at least one valve is moving is only half a period (i.e., the time for either closing or opening) for θ = 0, but is one full period for θ = −0.5.

We have so far studied advantages and disadvantages from varying the synchronisation of valves when the valves are moving quickly ( S = 80 ) ; we check whether valve synchronisation matters in slow valves ( S = 1 ). Specifically, when moving between any two of the three configurations depicted in the right panel of Figure 7, we vary the synchronisation θ between valves and measure σ , Q u m , Q ℓ m and Q _s . We find that by varying θ, we cannot increase σ , Q u m , Q ℓ m and Q _s by more than 1% above the maximum values they attain on the static configuration in Figure 5A, nor is back-flow induced. Thus with slow valves there are virtually no advantages/disadvantages to be had from varying the synchronisation of valves. Valve speeds are too slow to contribute to fluid velocity, and only valve position effects fluid velocity.

Finally, we examine how valve speed affects shear stress, illustrated in Figure 7C using the transition from the valve configuration with only upper valves open, to only lower valves open, with valves moving simultaneously. Figure 7Ci shows the instantaneous spatially-averaged shear stress during transition for three valve speeds, and Figure 7Cii shows the maximum and mean shear stress for S ∈ [ 0.1 , 80 ] . Increasing the valve speed, i.e., increasing S , increases the maximum instantaneous shear stress max(σ), at first gradually and then steeply. The time-averaged shear stress, σ ̄ , also increases, but more gently. For S = 80 , a maximum instantaneous shear stress of 0.027 Pa is attained, which is more than twice as large as the maximum shear stress of 0.011 Pa attained in the shear-promoting static configuration, with all other parameters kept constant.