Open Access

Research articles

Capacity and limitations of microfluidic flow to increase solute transport in three-dimensional cell cultures

Willy V. Bonneuil

Department of Engineering Mechanics, KTH Royal Institute of Technology, Stockholm, Sweden

Contribution: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Writing – original draft, Writing – review and editing

Google Scholar

Find this author on PubMed

Search for more papers by this author

Neeraj Katiyar

Department of Materials Science and Engineering, Uppsala University, Uppsala, Sweden

Science for Life Laboratory, Uppsala University, Uppsala, Sweden

Contribution: Investigation, Methodology, Visualization

Google Scholar

Find this author on PubMed

Search for more papers by this author

Maria Tenje

Department of Materials Science and Engineering, Uppsala University, Uppsala, Sweden

Science for Life Laboratory, Uppsala University, Uppsala, Sweden

Contribution: Conceptualization, Funding acquisition, Supervision, Writing – review and editing

Google Scholar

Find this author on PubMed

Search for more papers by this author

and

Shervin Bagheri

https://orcid.org/0000-0002-8209-1449

Department of Engineering Mechanics, KTH Royal Institute of Technology, Stockholm, Sweden

Contribution: Conceptualization, Formal analysis, Funding acquisition, Investigation, Methodology, Supervision, Writing – review and editing

Google Scholar

Find this author on PubMed

Search for more papers by this author

Willy V. Bonneuil^†

Willy V. Bonneuil

https://orcid.org/0000-0002-6147-5937

Department of Engineering Mechanics, KTH Royal Institute of Technology, Stockholm, Sweden

[email protected]

Contribution: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Writing – original draft, Writing – review and editing

Google Scholar

Find this author on PubMed

Search for more papers by this author

Neeraj Katiyar

Department of Materials Science and Engineering, Uppsala University, Uppsala, Sweden

Science for Life Laboratory, Uppsala University, Uppsala, Sweden

Contribution: Investigation, Methodology, Visualization

Google Scholar

Find this author on PubMed

Search for more papers by this author

Maria Tenje

Department of Materials Science and Engineering, Uppsala University, Uppsala, Sweden

Science for Life Laboratory, Uppsala University, Uppsala, Sweden

Contribution: Conceptualization, Funding acquisition, Supervision, Writing – review and editing

Google Scholar

Find this author on PubMed

Search for more papers by this author

and

Shervin Bagheri

https://orcid.org/0000-0002-8209-1449

Department of Engineering Mechanics, KTH Royal Institute of Technology, Stockholm, Sweden

Contribution: Conceptualization, Formal analysis, Funding acquisition, Investigation, Methodology, Supervision, Writing – review and editing

Google Scholar

Find this author on PubMed

Search for more papers by this author

Published:29 January 2025https://doi.org/10.1098/rsif.2024.0463

Abstract

Culturing living cells in three-dimensional environments increases the biological relevance of laboratory experiments, but requires solutes to overcome a diffusion barrier to reach the centre of cellular constructs. We present a theoretical and numerical investigation that brings a mechanistic understanding of how microfluidic culture conditions, including chamber size, inlet fluid velocity and spatial confinement, affect solute distribution within three-dimensional cellular constructs. Contact with the chamber substrate reduces the maximally achievable construct radius by 15%. In practice, finite diffusion and convection kinetics in the microfluidic chamber further lower that limit. The benefits of external convection are greater if transport rates across diffusion-dominated areas are high. Those are omnipresent and include the diffusive boundary layer growing from the fluid–construct interface and regions near corners where fluid is recirculating. Such regions multiply the required convection to achieve a given solute penetration by up to 100, so chip designs ought to minimize them. Our results define conditions where complete solute transport into an avascular three-dimensional cell construct is achievable and applies to real chambers without needing to simulate their exact geometries.

1. Introduction

The morphology and behaviour of living cells fundamentally change with the dimensionality of the environment they grow in. Three-dimensional cell cultures more accurately replicate the conditions that cancerous cells, migrating cells (e.g. immune cells) [1,2] or stem cells [3] experience in vivo [4]. The understanding of organ development and tumour growth has been furthered by cell-culture techniques that promote the self-assembly of cells into spheroids or organoids (stem-cell-derived three-dimensional tissue constructs) and where biological processes can be modulated by controlled physicochemical stimuli [5–7]. These techniques also show promise to accelerate drug testing [8]. Here, we define three-dimensional tissue constructs grown in controlled environments as multi-cellular engineered living systems (M-CELS). M-CELS acquire a broader physiological relevance as they grow larger because that lets them mimic more advanced development stages. Most M-CELS are engineered without a perfusable vasculature, which induces diffusion-limited growth. The incorporation of perfusable vasculatures into M-CELS is actively researched, but is complex and lacks a generally applicable strategy [9–11]. Nutrient transport and waste removal in avascular M-CELS rely on diffusion from their surface. Cells beyond a certain depth hence do not receive enough nutrients and become necrotic, i.e. they prematurely die. Necrosis affects the phenotype of the remaining live cells, e.g. by releasing inflammatory or chemo-attractant markers into the intercellular space [12] and altering the mechanical properties of the live rim [13]. The balance between proliferating cells on the surface of M-CELS and dying cells at their centre confers M-CELS a diffusion-limited maximum diameter [14]. Both phenotype alteration and size limitation restrict the stage of tumour growth or organ development up to which avascular M-CELS are a relevant biological model.

Microfluidic technology has enabled the culture of avascular M-CELS of larger diameters and greater compartmentalization, i.e. of more advanced development [15–18]. The fluid velocities over M-CELS vary widely, as have the M-CELS-confining structures which are necessary for their retention. The magnitudes of the improvements brought about by microfluidic culture have been mixed. In one brain-organoid study, organoid diameters increased by a fifth under a steady fluid flow of 0.1 mm s⁻¹ (our calculation, from parameters listed in the paper), without necrotic-core prevention [15]. In another brain-organoid study, diameters increased by a third under the oscillatory fluid flow of 1 mm s⁻¹ (our calculation), with necrotic-core prevention [18]. Fluid flow has been treated as an ‘on–off’ parameter in many studies comparing static with microfluidic culture. Its varying effects encourage the mechanistic elucidation of its influence on solute transport in M-CELS.

Previous modelling efforts on nutrient transport into M-CELS have yielded analytical expressions for nutrient concentration in systems with spherical [19,20] or ellipsoidal [21] symmetry, isolated from their culture environment. Numerical work has been conducted on nutrient supply to spheroids in specific culture set-ups, including culture traps [22] or confinement pillars [23]. The dimensions of defined chamber geometries have been optimized to increase nutrient supply to M-CELS of a given cell type [24,25]. These numerical studies have not established general principles of chamber design and operation by which to harness the benefits of microfluidic flow and which can be quantitatively applied at the experiment-design stage. Here, we achieve this by combining theoretical and numerical analyses of nutrient supply to M-CELS. This study accounts for the distance between M-CELS and a nutrient source and M-CELS structural confinement. It makes the influence of chamber geometry and flow rate on solute transport and necrosis occurrence more predictable. Although nutrient supply is a primary interest, this investigation applies to any solutes, including drugs, that diffuse into an M-CELS and whose freely diffusing form has a sink term: consumption or binding to an agonist receptor. We estimate not only the conditions for the prevention of necrosis, but more generally those for an increase in penetration distance. Our results highlight the penalizing role of large diffusion-dominated areas such as corners or hydrogel capsules and of fast diffusion within M-CELS, which both delay the solute-transport enhancement brought by convection.

2. Methods

2.1. Computational domains

Let a computational domain $Ω$ comprise an M-CELS as an inclusion $Ω_{2}$ in a straight fluid channel $Ω_{1}$ , such that $Ω_{1} \cup Ω_{2} = Ω$ . The interface between the domains is noted ${\bar{Ω}}_{1} \cap {\bar{Ω}}_{2} = Γ$ . The fluid channel has an inlet $I$ and an outlet $O$ . The channel walls were noted $W_{0}$ at $y = 0$ and $W_{h}$ at $y = 2 h$ and their union $W$ . $Ω_{2}$ was assumed porous and rigid with homogeneous properties (table 1).

Table 1. Symbol table.
Collapse
symbol	meaning
$Ω_{1}$	microfluidic channel
$Ω_{2}$	M-CELS
${\hat{Ω}}_{i}$	domain $i$ in the radial-linear analytical geometry
$Γ$	fluid–M-CELS interface
$I$	channel inlet
$O$	channel outlet
$W$	wall boundary
$S$	symmetry boundary
$v$	fluid velocity
$p$	fluid pressure
$v_{0}$	peak inlet fluid velocity
$a$	M-CELS radius
$L$	distance between M-CELS and channel inlet
$k$	Darcy permeability (homogeneous to a surface area)
$ρ$	fluid density
$μ$	dynamic viscosity
$c$	solute concentration
$c_{0}$	inlet solute concentration
$D_{i}$	solute diffusivity in domain $Ω_{i}$
$R_{\max}$	maximum solute consumption rate
$Φ_{L}$	fraction of cells receiving enough solute in $Ω_{2}$
subscript 0	relative to the no-fluid-flow asymptote
subscript $\infty$	relative to the maximal-supply asymptote

To isolate the effects of M-CELS confinement, we simulated two geometries. First, an ‘unconfined’ geometry where $Ω_{2}$ is a disc surrounded by fluid (figure 1) and where $Ω$ is bisected by a symmetry line $S$ . The flow and the M-CELS are confined by the channel walls. The term ‘unconfined’ refers to the M-CELS having no contact with any solid. Second, a confined geometry where $Ω_{2}$ is a disc cut such that its height from the bottom of the channel is 90% of its diameter (figure 2). $Ω_{2}$ is in contact with $W_{0}$ through its linear boundary. This confinement is ‘geometrically minimal’ in that it includes fluid-recirculation areas upstream and downstream of obstacles (in $Ω_{1}$ close to the contact points $Ω_{1} \cap Ω_{2} \cap W$ ) and an area of contact between the M-CELS and a solid wall, but does not include wall-normal obstacles. The surface area of contact $Ω_{2} \cap W_{0}$ corresponds to the one observed in an experimental human lung-fibroblast spheroid (PB-CH-450-0811, PELOBiotech, Germany) cultured between confinement pillars in a microfluidic channel (electronic supplementary material, figure S1). In both geometries, $Ω$ is described by a Cartesian coordinate system $(x, y)$ centred on the centre $O$ of the smallest circle containing $Ω_{2}$ . A radial coordinate system $(r, θ)$ is also defined from $O$ .

Figure 1. Computational domain (unconfined). M-CELS as a disc $Ω_{2}$ of radius $a$ in a rectangular channel $Ω_{1}$ . $Ω$ bisected by a symmetry line $S$ , parallel to the direction of flow. Coordinate origin at the centre of the M-CELS.

Figure 2. Computational domain (confined). M-CELS as a partial disc $Ω_{2}$ cut such that its height above the bottom wall of the fluid channel is 90% of its diameter. Coordinate origin at the centre of the smallest circle containing $Ω_{2}$ .

2.2. Quasi-one-dimensional analytical domain

To acquire an analytical understanding of the problem, we define a quasi-one-dimensional ‘radial-linear’ domain $\hat{Ω}$ as the union of a radial section of a disc of radius $a$ and angle $d θ$ , named ${\hat{Ω}}_{2}$ and centred on $\hat{O}$ , and the constant-width extension of that section, named ${\hat{Ω}}_{1}$ (figure 3). $\hat{Γ}$ and $\hat{I}$ and the coordinate systems $(O, \hat{x}, \hat{y})$ and $(O, \hat{r}, \hat{θ})$ are defined in analogous ways to their two-dimensional counterparts. The planes ${\hat{S}}_{r} = {\hat{θ} = \pm d θ / 2}$ are symmetry planes for ${\hat{Ω}}_{2}$ and the planes ${\hat{S}}_{l}$ of normal $\hat{y}$ that bound ${\hat{Ω}}_{1}$ are symmetry planes for ${\hat{Ω}}_{1}$ . If $d θ ≪ 1$ , we assume that, for a scalar $ξ$

Quasi−one-dimensional ``radial-linear’ reduction — Figure 3. Quasi-one-dimensional ‘radial-linear’ reduction of $Ω$ into $\hat{Ω}$ . The angular section ${\hat{Ω}}_{2}$ has an angle $d θ ≪ 1$ .

{\frac{d ξ}{d \hat{r}} \cdot \hat{r} |}_{\hat{Γ}} = {\frac{d ξ}{d \hat{x}} \cdot \hat{x} |}_{\hat{Γ}} .

(2.1)

The assumptions of symmetry and narrowness of ${\hat{Ω}}_{2}$ allow equations in $\hat{Ω}$ to be reduced to a quasi-one-dimensional ‘radial-linear’ form where variables depend only on $\hat{x}$ in ${\hat{Ω}}_{1}$ and only on $\hat{r}$ in ${\hat{Ω}}_{2}$ .

2.3. Equations for fluid flow

Fluid flow in $Ω_{1}$ is modelled by the steady Navier–Stokes equations with a parabolic inlet velocity and normal outlet flow (electronic supplementary material, Section M1). The normalized momentum conservation equation is

R e (v_{1}^{*} \cdot \nabla^{*}) v_{1}^{*} = - \nabla^{*} p^{*} + \nabla^{*} \cdot (\nabla^{*} v_{1}^{*} + \nabla v_{1}^{* T}),

(2.2)

where the Reynolds number has been defined as $R e = ρ v_{0} a / μ$ , local acceleration has been neglected, and non-dimensional variables are noted with asterisks. Microfluidic conditions involve peak inlet velocities of the order of 1 mm s⁻¹, M-CELS radii are at most of the order of 1 mm, and cell culture medium has similar density and viscosity to those of water [26], so $R e ≲ 1$ . We therefore did not neglect convective acceleration.

Fluid flow in $Ω_{2}$ is modelled by the Darcy equations,

\begin{aligned} 0 = \frac{k}{μ} \nabla p_{2} + v_{2} \end{aligned},

(2.3a)

\begin{aligned} \nabla \cdot v_{2} & = 0 \end{aligned} .

(2.3b)

We assumed that the M-CELS had homogeneous properties, which is a valid assumption at least at the early stages of M-CELS growth when little cellular differentiation has occurred. Continuity of pressure and velocity are prescribed at the fluid–porous interface,

\begin{aligned} p_{1} & = {p_{2} |}_{Γ}, \end{aligned}

(2.4a)

\begin{aligned} v_{1} & = {v_{2} |}_{Γ} . \end{aligned}

(2.4b)

Pressure continuity may be applied at a fluid–porous interface when there is a separation of scales in the porous domain (here, the inter-cell distance of approximately 10 nm is much smaller than the M-CELS radius of at least 100 µm [27]). The continuity of velocity between its Navier–Stokes and Darcy solutions results from mass conservation and the assumption that the slip length in $Ω_{2}$ is zero. The latter scales with $\sqrt{k}$ so is at most 100 nm and therefore, much smaller than a cell (M-CELS are expected to be less permeable than the most permeable tissues, where $k \sim 10^{- 14} m^{2}$ [28]).

2.4. Equations for nutrient transport

Solute transport is simulated by steady advection–diffusion–reaction equations. For $i \in {1,2}$ ,

0 = \nabla_{i} \cdot (D_{i} \nabla c_{i}) - \nabla_{i} \cdot (v_{i} c_{i}) - R (c_{i}) .

(2.5)

The velocities $v_{i}$ are obtained from the fluid-flow equations. Diffusivities are uniform in each domain, at $D_{i}$ . Solute consumption is 0 in $Ω_{1}$ and follows Michaelis–Menten kinetics in $Ω_{2}$ , as in models of tumour spheroids [23,29],

R (c_{2}) = \frac{R_{m a x} c_{2}}{c_{2} + K_{1 / 2}} .

(2.6)

The consumption rate increases continuously from 0 in the absence of solute to $R_{\max}$ at high solute concentration ( $c_{2} ≫ K_{1 / 2}$ ). The half-rate constant $K_{1 / 2}$ verifies $R (K_{1 / 2}) = R_{\max} / 2$ . The consumption parameters are uniform in the M-CELS, which is reasonable at the early stages of development. An unlimited solute supply, $c = c_{0}$ , was assumed at $I$ and $O$ . This represents the common placement of a culture-medium reservoir on either side of a microfluidic chamber.

Let non-dimensionalized variables, noted with asterisks, be defined as $x_{1} = x_{1}^{*} L$ , $x_{2} = x_{2}^{*} a$ , $v_{1} = v_{1}^{*} v_{0}$ and $c_{i} = c_{i}^{*} c_{0}$ for $i \in {1,2}$ . To non-dimensionalize $v_{2}$ , we combine the pressure drop across a sphere in Stokes’ flow ( $Δ p_{S} \sim μ v_{0} / a$ ) and Darcy’s law ( $v_{2} \sim (k / μ) (Δ p_{S} / a)$ ) to define $v_{2} = v_{2}^{*} k v_{0} / a^{2}$ . We also define $r = r^{*} a$ and $K_{1 / 2} = K_{1 / 2}^{*} c_{0}$ . In $Ω_{1}$ , the non-dimensionalized transport equation is

0 = \nabla_{1}^{* 2} c_{1}^{*} - {P e}_{1} \nabla^{*} \cdot (v_{1}^{*} c_{1}^{*}),

(2.7)

where the Péclet number is ${P e}_{1} = v_{0} L / D_{1}$ . In $Ω_{2}$ , the non-dimensionalized transport equation is

0 = \nabla_{2}^{* 2} c_{2}^{*} - {P e}_{2} \nabla \cdot (v_{2}^{*} c_{2}^{*}) - D a \frac{c_{2}^{*}}{c_{2}^{*} + K_{1 / 2}^{*}},

(2.8)

where the Péclet number is ${P e}_{2} = k v_{0} / (a D_{2})$ and the Damköhler number is $D a = R_{\max} a^{2} / (c_{0} D_{2})$ .

The most favourable situation for convection in $Ω_{2}$ shows the latter is negligible. Microfluidic flows rarely have a peak velocity higher than 1 mm s⁻¹, the permeability of M-CELS is assumed under 10⁻¹⁴ m² (that of the most permeable tissues [28]), and the question of solute transport sufficiency in M-CELS usually arises when those are larger than 100 μm. The Stokes–Einstein diffusivity in water of solutes of interest in M-CELS development is above $10^{- 11} m^{2} s^{- 1}$ (which is the value for proteins of 100 kDa). We assume that $D_{2} / D_{1} > 0.01$ for non-lipophilic solutes (those that do not cross cellular membranes), so $D_{2} > 10^{- 13} m^{2} s^{- 1}$ . The value of 0.01 is the order of the lowest effective glucose diffusivity over a range of tumour spheroids [30]. This yields ${P e}_{2} < 0.1$ , so solute transport in $Ω_{2}$ may be described by a diffusion–reaction equation,

0 = \nabla_{2}^{* 2} c_{2}^{*} - D a \frac{c_{2}^{*}}{c_{2}^{*} + K_{1 / 2}^{*}} .

(2.9)

A no-flux condition is prescribed on the symmetry plane and walls,

{(\nabla c_{i} \cdot n) |}_{W \cup S} = 0, i \in {1, 2} .

(2.10)

At the interface $Γ$ , conservation of nutrient quantity implies continuity of concentration. Let the interface concentration be noted $γ$ ,

{c_{1} |}_{Γ_{1}} = {c_{2} |}_{Γ_{2}} = γ .

(2.11)

The evaluation of a quantity on $Γ_{i}$ is understood as the limit of that quantity when approaching $Γ$ from $Ω_{i}$ . Conservation of nutrient quantity also implies continuity of normal diffusive fluxes,

D_{1} {\frac{\partial c_{1}}{\partial n} \cdot n |}_{Γ_{1}} = D_{2} {\frac{\partial c_{2}}{\partial n} \cdot n |}_{Γ_{2}},

(2.12)

with $n$ the normal to $Γ$ (of arbitrary orientation).

2.5. Nutrient balance at the fluid–multi-cellular engineered living systems interface

Three additional non-dimensional quantities are defined to describe physical balances of importance on $Γ$ .

The ratio of diffusive solute fluxes along the normal $n$ (oriented into $Ω_{2}$ ) is

\begin{aligned} \frac{J_{1}}{J_{2}} & = \frac{- D_{1} {\frac{\partial c_{1}}{\partial n} \cdot n |}_{Γ_{1}}}{- D_{2} {\frac{\partial c_{2}}{\partial n} \cdot n |}_{Γ_{2}}} \\ = \frac{a D_{1}}{L D_{2}} \frac{\frac{\partial c_{1}^{*}}{\partial n^{*}} \cdot n^{*}}{\frac{\partial c_{2}^{*}}{\partial n^{*}} \cdot n^{*}} . \end{aligned}

(2.13)

This ratio is equal to 1, so the ratio of interface concentration gradients is characterized by

\frac{\frac{\partial c_{2}^{*}}{\partial n^{*}} \cdot n^{*}}{\frac{\partial c_{1}^{*}}{\partial n^{*}} \cdot n^{*}} = \frac{a D_{1}}{L D_{2}} .

(2.14)

Let us define

R_{d} = \frac{a D_{1}}{L D_{2}},

(2.15)

which characterizes whether diffusion into $Ω_{2}$ is supply-limited (low $R_{d}$ ) or rate-limited (high $R_{d}$ ).

Let us approximate $γ$ as uniform to consider the solute quantity balance in a radial section of $Ω_{2}$ of angle $d θ$ centred on $O$ . The corresponding section of $Γ$ receives

n_{d} = - D_{1} a d θ \frac{d c_{1}}{d n} \cdot n

(2.16)

by diffusion and

n_{c} = a d θ γ v_{1} \cdot n

(2.17)

by convection. Let $r_{n}$ be the radius where $c_{2}$ falls to zero. In a radially symmetric geometry, $Γ_{n}$ is the circle $r = r_{n}$ . We allow $r_{n} = 0$ if $c_{2}$ is positive everywhere. The radial section loses the following solute quantity due to consumption:

\begin{aligned} n_{R} & = R_{m a x} \int_{r_{n}}^{a} r d θ d r \\ = \frac{R_{m a x}}{2} (a^{2} - r_{n}^{2}) d θ . \end{aligned}

(2.18)

Let us now normalize the solute quantities and consider the ratio between diffusive supply to $Γ$ and consumption in $Ω_{2}$ :

\frac{n_{d}}{n_{R}} = \frac{2 D_{1} c_{0}}{R_{m a x} a L} \frac{1 - r_{n}^{* 2}}{\frac{d c_{1}^{*}}{d n} \cdot n},

(2.19)

from which we define the diffusive supply number as

S_{d} = \frac{2 c_{0} D_{1}}{R_{m a x} a L} .

(2.20)

$S_{d}$ describes how the diffusive supply of solutes to $Γ$ affects the quantity of solutes consumed in $Ω_{2}$ . If $S_{d} < 1$ , solute consumption is supply-limited, i.e. the quantity consumed by the M-CELS is limited by the quantity that it receives from its environment. If $S_{d} > 1$ , solute consumption is rate-limited, i.e. the M-CELS consumes as much as its biology demands. $S_{d}$ may be expressed as

S_{d} = \frac{2 R_{d}}{D a} .

(2.21)

Let us consider the ratio between convective solute supply to $Γ$ and consumption in $Ω_{2}$ ,

\frac{n_{c}}{n_{R}} = \frac{2 c_{0} v_{0}}{R_{m a x} a} \frac{1 - r_{n}^{* 2}}{γ^{*} v_{1}^{*} \cdot n},

(2.21)

from which we define the convective supply number as

S_{c} = \frac{2 c_{0} v_{0}}{R_{m a x} a} .

(2.23)

$S_{c}$ may be described in an analogous way to $S_{d}$ and expressed as

S_{c} = \frac{2 R_{d} {P e}_{1}}{D a} .

(2.24)

The ratios $R_{d}$ , $S_{d}$ and $S_{c}$ are measures of the excess or deficit of solute supply (diffusive or convective) relative to solute kinetics in the M-CELS (diffusive or consumptive).

2.6. Numerical methods

The simulations of coupled fluid flow and nutrient transport were implemented in the commercial software COMSOL Multiphysics^TM (COMSOL AB, Sweden, 6.1). The electronic supplementary material, Section M2 presents the details of that implementation and the dimensional parameter values.

2.7. Measures of solute transport

The goal of the simulations is to quantify the sufficiency of solute transport to cells in the M-CELS. That is measured by considering the shape and the surface area of the subset of $Ω_{2}$ where $c_{2}$ is sufficiently high. The boundary of sufficient transport (or necrotic boundary if the solute is a nutrient) is defined as

Γ_{n} = {x_{2} | c_{2} = c_{n}},

(2.25)

where $c_{n}$ is a threshold concentration that depends on the metabolic needs of the cells. We use $c_{n} = 10^{- 3} c_{0}$ for generality. The centre of mass of the zone of insufficient transport (or necrotic centre if the solute is a nutrient) has the coordinates

{{\bar{x}}_{n}, {\bar{y}}_{n}} = {\frac{\int_{c_{2} < c_{n}} x_{2} d A}{\int_{Ω_{2}} d A}, \frac{\int_{c_{2} < c_{n}} y_{2} d A}{\int_{Ω_{2}} d A}} .

(2.26)

The fraction of sufficient transport is defined as the fraction of cells receiving solutes at a concentration above $c_{n}$ . It is termed the live fraction if the solute is a nutrient and is defined as

Φ_{L} = \frac{\int_{c_{2} > c_{n}} d A}{\int_{Ω_{2}} d A} .

(2.27)

For simplicity and without loss of generality, the terms associated with nutrients are used in the following sections.

2.8. Definition of asymptotes

The analysis is constructed as follows and depicted in figure 4. We first consider the effect of M-CELS confinement independently of solute supply, i.e. at $γ^{*} = 1$ . This is the maximal-supply asymptote and is indicated by an $\infty$ subscript. It is the theoretical maximum of solute transport that is achievable in a cell-culture device. We then consider the effect of finite diffusive transport in static culture, termed the static asymptote and indicated by a $0$ subscript. We finally consider the role of convective supply, i.e. the trajectory that $Γ_{n}$ and $Φ_{L}$ follow between their static and maximal-supply asymptotes as convective supply increases. The maximal-supply asymptote determines whether sufficient transport is achievable in the whole M-CELS. If that is the case, any insufficiency is termed supply-induced, and a condition is derived for its prevention. If not, insufficiency is termed rate-induced.

Figure 4. Flow chart of solute-transport analysis in M-CELS set in a culture chamber. Superscripts $u$ and $c$ respectively refer to ‘unconfined’ and ‘confined’. Capital letters refer to effects captured by the simulations but not by the analytical expressions. A: mathematical approximation and Michaelis–Menten consumption kinetics. B: two-dimensional effects, mostly the constriction between $Ω_{2}$ and $W_{h}$ . C: asymmetry of fluid flow in $Ω_{1}$ and fluid-recirculation zones near the confinement corners.

3. Results

3.1. Maximal-supply asymptote

At the radially symmetric maximal-supply asymptote, the $C^{1}$ -regularity of $c_{2}$ implies that $c_{2}$ is zero with a zero gradient at the necrotic radius $r_{n}$ . This introduction of $r_{n}$ has been done by others [19] and was shown to satisfactorily approximate Michaelis–Menten kinetics if $K_{1 / 2}^{*} ≪ 1$ [29]. This is generally verified for cellular oxygen consumption. The solution of the transport equation (2.9) approximated with a constant consumption rate yields the following equation for $r_{n}^{*} = r_{n} / a$ (details in electronic supplementary material, Section R1):

1 = \frac{D a}{4} (1 - r_{n}^{* 2}) + \frac{D a}{2} r_{n}^{* 2} \ln (r_{n}^{*}) .

(3.1)

If $r_{n}^{*} ≪ 1$ , then $- r_{n}^{* 2} \ln (r_{n}^{*}) ≫ r_{n}^{* 2}$ . Equation (3.1) simplifies into a second-order polynomial of $r_{n}^{*} \sqrt{- \ln (r_{n}^{*})}$ . It admits a real solution if $D a \geq 4$ , which is positive and tends to 0 when $D a ⟶ 4$ . Let ${D a}_{\infty} = 4$ be the onset of rate-induced necrosis. The numerical simulations verified this in the unconfined geometry ( $Φ_{L, \infty} < 1$ for $D a > 4$ ) and showed ${D a}_{\infty} = 3.4$ in the confined geometry (figure 5a). The penalization of cell survival by confinement translates into the maximum non-necrotic radius being 92% of the equivalent unconfined M-CELS.

Figure 5. Live fraction at maximal supply (a) around the onset of necrosis, which occurs at lower $D a$ in confinement, and (b) closely matching its asymptotic expression for $D a > 10$ (error under 5%).

If the necrotic area is large, $r_{n}^{*} \approx 1$ and $D a ≫ 1$ . Let $δ^{*} = 1 - r_{n}^{*}$ . The logarithm term may be expanded around $δ^{*} = 0$ as

\begin{aligned} \ln (r_{n}^{*}) & = \ln (1 - δ^{*}) \\ = - δ^{*} - \frac{δ^{* 2}}{2}, \end{aligned}

which yields the following expansion of equation (3.1) at order 2 in $δ^{*}$ :

\frac{D a}{2} δ^{* 2} - 1 = 0.

(3.2)

The asymptotic expressions for $r_{n}^{*}$ and $Φ_{L}^{\infty}$ at high $D a$ follow:

\begin{aligned} r_{n}^{*} & = 1 - {(\frac{2}{D a})}^{1 / 2}, \end{aligned}

(3.3)

\begin{aligned} Φ_{L}^{\infty} & = 2 (- \frac{1}{D a} + {(\frac{2}{D a})}^{1 / 2}) . \end{aligned}

(3.4)

The latter is an excellent predictor of the unconfined live fraction for $D a \geq 10$ (error under 5%, figure 5b). The confined live fraction in the asymptotic high- $D a$ regime was lower by around 15% for $D a \geq 10$ . This translates into necrotic radii being around 9% larger than in the equivalent unconfined M-CELS.

Figure 6 shows the shape of necrotic cores in confined M-CELS for varying $D a$ . The contact surface skews the cores towards it. Those appear along the surface’s normal bisector when $D a$ becomes larger than ${D a}_{\infty}$ . Necrosis encompasses the M-CELS barycentre from its onset (see $Γ_{n}$ at $D a = 4$ ). The normal to the contact surface passing by the barycentre thus is the region most vulnerable to necrosis.

Figure 6. Necrotic boundaries $Γ_{n}$ at maximal supply in confined M-CELS, with necrotic tissue in the area enclosed by $Γ_{n} \cup W_{0}$ . Necrosis grows perpendicularly from the confining surface for $D a \sim 1$ and increases before expanding in all directions at higher $D a$ . Positions of necrotic-core barycentres marked by $\times$ .

3.2. Static asymptote

Let us obtain an analytical approximation of $Φ_{L}^{0}$ on the quasi-one-dimensional domain $\hat{Ω}$ . The non-dimensional necrotic radius ${\hat{r}}_{n}^{*} = {\hat{r}}_{n} / a$ satisfies this equation,

- \frac{1}{S_{d}} (1 - {\hat{r}}_{n}^{* 2}) + 1 = \frac{D a}{4} (1 - {\hat{r}}_{n}^{* 2}) + \frac{D a}{2} {\hat{r}}_{n}^{* 2} \ln ({\hat{r}}_{n}^{*}),

(3.5)

whose derivation is detailed in electronic supplementary material, Section R2.

3.2.1. Prevention of necrosis

Let us investigate the prevention of supply-induced necrosis at $D a < {D a}_{\infty}$ . When the necrotic core is small, $| \ln ({\hat{r}}_{n}^{*}) | ≫ 1$ , so (3.5) is approximated as

- \frac{D a}{2} {\hat{r}}_{n}^{* 2} \ln ({\hat{r}}_{n}^{*}) = \frac{1}{S_{d}} + \frac{D a}{4} - 1 .

(3.6)

This admits a real solution if the right-hand term is positive. The necrotic core vanishes when that term becomes negative, i.e. with the unconfined ${D a}_{\infty} = 4$ ,

S_{d} (1 - \frac{D a}{{D a}_{\infty}}) \geq 1 .

(3.7)

This satisfies physical intuition about extreme values of $D a$ . When $D a$ is low, it suffices that as many nutrients arrive at the interface as the M-CELS consumes ( $S_{d} > 1$ ). As $D a$ increases, $γ^{*}$ needs to become ever closer to 1 to prevent supply-induced necrosis, which requires an excess of nutrient supply and $S_{d} > 1$ . Finally, the value of $S_{d}$ required for necrosis prevention diverges as $D a$ approaches the onset of rate-induced necrosis. The condition (3.7) is expressed using the velocity scale of diffusion, which is $v_{dif} = 2 D / d$ at a diffusivity $D$ across a domain of length $d$ ,

v_{d i f, 1} \geq \frac{v_{d i f, 2}}{2 (\frac{1}{D a} - \frac{1}{{D a}_{\infty}})} .

(3.8)

Necrosis prevention by high diffusive solute supply is limited by fast diffusion in the M-CELS. Preventing necrosis indeed requires raising the concentration on the fluid–M-CELS interface and that concentration is lower when the diffusive-flux balance favours transport away from the interface, into the M-CELS. The diffusive-supply requirement is more stringent for small M-CELS with a fast metabolism than for large M-CELS with a slow metabolism. This consideration applies to oxygen, whose diffusivity is of the order of $10^{- 9} m^{2} s^{- 1}$ in tissues but whose consumption rate varies over two orders of magnitude between cell types ([20,31]).

The design constraints of experimental set-ups in practice lower the onset of rate-induced necrosis. They include space constraints on microscope stages and, often, room for the culture of multiple M-CELS. These both impose a minimal value on $L$ of at least a few millimetres. Consequently, considering that $D_{2} / D_{1}$ is typically between 0.01 and 1, $R_{d}$ is constrained under 10, and under 1 for nutrients like oxygen which readily cross cell membranes and for which $D_{2} / D_{1} \approx 1$ . For oxygen, the maximal $D a$ at which adequate nutrient supply can be achieved in static culture is $D a \approx 1$ . There, the condition (3.7) implies $R_{d} > 2 / 3$ , which is close to its limit.

The condition (3.7) is verified by numerical simulations of the unconfined geometry (figure 7a). Necrosis prevention is similarly achieved in the confined geometry, with $S_{d} (1 - D a / {D a}_{\infty}) > 1$ and ${D a}_{\infty} = 3.4$ (figure 7b). An approximately twice higher diffusive flux to $Γ$ is required to prevent necrosis relative to the unconfined geometry. This corresponds to a reduction of the maximal permissible distance between a solute reservoir and an M-CELS by half in static culture, if all other parameters are equal.

Figure 7. Live fraction at the static asymptote, for $D a \leq 3$ and $R_{d}$ varying between 0.1 and 100, in the (a) unconfined and (b) confined geometry. It suffices that the excess of diffusive solute supply satisfies $S_{d} (1 - D a / {D a}_{\infty}) > M$ for supply-induced necrosis to be prevented, with $M = 1$ in the unconfined geometry and $M = 2$ in the confined one.

3.2.2. Approach of maximal-supply asymptote

Let us investigate the live fraction as a function of diffusive supply when there is rate-induced insufficient transport and diffusive supply is ample, i.e. $S_{d} ≫ 1$ , and $R_{d} ≫ 1$ . The live fraction normalized by its maximal-supply value is analytically derived in electronic supplementary material, Section R3 as

\frac{Φ_{L}^{0}}{Φ_{L}^{\infty}} = - (R_{d} S_{d})^{- 1 / 2} + {(1 + (R_{d} S_{d})^{- 1})}^{1 / 2} .

(3.9)

This describes a sigmoid function of $S_{d}$ with an inflection point—the point where $Φ_{L}^{0} = Φ_{L, \infty} / 2$ —shifted by $R_{d}$ . As for necrosis prevention, the excess of diffusive supply relative to consumption needs to be larger when diffusion into the M-CELS is supply-limited. The variable $R_{d} S_{d}$ may be written

R_{d} S_{d} = \frac{2 c_{0} D_{1}^{2}}{R_{m a x} L^{2} D_{2}} .

(3.10)

Counter-intuitively, in static culture, the proximity of an M-CELS with high $D a$ to its maximal-supply asymptote does not depend on the M-CELS radius. It therefore does not depend on the M-CELS age if transport properties remain constant. The expression (3.9) also informs on how to adapt experimental parameters ( $c_{0}$ and $L$ ) based on those that may evolve with M-CELS development and health ( $R_{\max}$ and $D_{2}$ ). For example, the tight alignment of neurons of the outer cell layer of cortical-brain organoids [32] increases $D_{2}$ because the tortuosity becomes close to 1. If M-CELS diffusivity becomes $D_{2}^{'} = λ D_{2}$ , solute concentration needs to become $c_{0}^{'} = λ c_{0}$ or the distance between the solute source and $Ω_{2}$ needs to become $L^{'} = λ^{- 2} L$ to maintain solute penetration.

Equation (3.9) is a good predictor of the simulated live fraction for $D a \geq 4$ in both the unconfined and confined geometries (figure 8). It overestimates the live fraction at low $S_{d}$ because it does not account for the constriction between $Γ$ and $W_{h}$ . Solutes at low $S_{d}$ only reach the parts of $Ω_{2}$ closest to $I$ and $O$ , so the uniform- $γ$ assumption is less accurate. The live fraction progressed towards its maximal-supply asymptote along a similar trajectory in the confined geometry but with a delay. $S_{d}$ needed to be approximately 50% higher for the relative live fraction to progress to the same stage as if unconfined.

At the static asymptote, the distance between the solute reservoir and the M-CELS implies that supply to the M-CELS surface is reduced. Preventing necrosis for $D a$ close to the onset of rate-induced necrosis requires an excess of diffusive nutrient supply relative to nutrient consumption that is experimentally hard to achieve. The practical maximum $D a$ in static culture can be lower than 1, which means a division by 2 of the maximum non-necrotic radius compared with maximally supplied M-CELS. For M-CELS with rate-induced necrosis, the distance to the maximal-supply asymptote does not depend on their radius, so the efficacy of solute supply in a no-flow set-up is constant for M-CELS with age-independent transport properties.

3.3. Microfluidic culture

3.3.1. Function of increased convective supply

Convective solute supply only affects transport in the M-CELS if there is a deficiency at the static asymptote, i.e. $S_{d} ≲ 1$ . Then, convection in $Ω_{1}$ fulfils three distinct functions depending on $D a$ (figure 9). If $D a ≪ 1$ , concentration in $Ω_{2}$ is approximately uniform. Any necrosis is supply-induced and can be prevented. It suffices that $S_{c} > 1$ . If $D a ≫ 1$ , most of the M-CELS is affected by rate-induced necrosis. External convection would increase solute supply only to cells that are closest to $Γ$ . The penetration distance of solutes would increase at most until its value at the unconfined maximal-supply asymptote $\sqrt{2 / D a}$ , which still is poor transport. If $D a \sim 1$ , necrosis may be supply-induced or rate-induced, according to a threshold that depends on M-CELS confinement. Analogously to the static asymptote, supply-induced necrosis can be prevented by $S_{c} (1 - D a / 4) ≳ 1$ . Rate-induced necrosis cannot be prevented but may be reduced by high $S_{c}$ . External convection reduces the distance across which solutes need to diffuse to reach $Γ$ , so it needs to increase when diffusive transport across areas of low fluid velocity is slow. Those areas are the diffusive boundary layer over $Γ$ and the vicinity of the confinement corner ( $Γ \cap Ω_{2}$ ). Since $R_{d} S_{d}$ characterizes the efficacy of diffusive supply, we postulate that $S_{c} \sim (R_{d} S_{d})^{- 1}$ to have a microfluidic system approach the maximal-supply asymptote.

Figure 9. Qualitative effect of increased convective nutrient supply to $Ω_{2}$ . For deficient diffusive supply (low $S_{d}$ ), an ample convective supply may prevent supply-induced necrosis at low $D a$ ( $D a < {D a}_{\infty}$ ), but cannot prevent rate-induced necrosis and would only marginally increase nutrient supply to the outer cell layers at $D a ≫ {D a}_{\infty}$ .

3.3.2. Evolution of solute supply with increasing convection

Let us describe how progressively higher convection in $Ω_{1}$ affects solute supply to $Γ$ and its concentration in $Ω_{2}$ . At the static asymptote, solutes diffusing from $I$ and $O$ form a symmetric concentration field around $Ω_{2}$ (figure 10a). Low convection from $I$ to $O$ ( ${P e}_{1} > 1$ but $S_{c} < 1$ ) enhances solute transport from $I$ but inhibits it from $O$ (figure 10b). The interface concentration $γ$ increases on the upstream side of $Γ$ but decreases on its downstream side. The area of insufficient supply is shifted towards $O$ . The supply increase on the upstream side of $Γ$ mostly concerns its portion that is not affected by the confinement corner. Fluid flow near the latter indeed occurs as an infinite sequence of recirculation eddies of decreasing sizes and intensities [33]. Those were observed in the confined geometry because of the very large resistance to flow of $Ω_{2}$ , even though fluid transpires across $Γ$ . Let $Γ_{ρ}$ be the subset of $Γ$ exposed to the recirculation eddies (electronic supplementary material, figure S4). Solute transport in that region is diffusion-dominated unless the convection in $Ω_{1}$ is very large (the maximum velocity on the recirculation eddies decreases geometrically, with a ratio of approximately $10^{- 3}$ for the 60° angle of the confined geometry [33]). The efficacy of static solute supply thus remains determinant for the regions of $Ω_{2}$ near the confinement corner even when transport is convection-dominated in most of $Ω_{1}$ . As convection further increases ( $S_{c} > 1$ ), the solute flux from $I$ becomes sufficient to fill the supply deficiency on the downstream side of $Γ$ exposed to the primary flow and any remaining necrosis again becomes symmetric about $x = 0$ (figure 10c). Solute supply increases to the confinement corner as convection further increases and transport in the largest recirculation eddy eventually becomes convection-dominated (figure 10d).

Figure 10. Schematic evolution of concentration in $Ω_{1}$ under increasing convection from $I$ to $O$ . (a) Diffusive equilibrium with a symmetric concentration around ${x = 0}$ . (b) Low convection ( $S_{c} \sim 1$ ) increases concentration on the upstream side of $Ω_{2}$ but decreases it on its downstream side. The upstream increase is limited near the confinement corner. (c) Convection of higher magnitude ( $S_{c} ≫ 1$ ) fills the primary flow region in $Ω_{1}$ with solute ( $c \approx 1$ there), but not the confinement corners. (d) At ever higher convection, concentration increases in the confinement corner due to increased corner-eddy intensity.

3.3.3. Prevention of necrosis

In the quasi-one-dimensional domain $\hat{Ω}$ , the non-dimensional necrotic radius ${\hat{r}}_{n}^{*} = {\hat{r}}_{n} / a$ satisfies the following equation:

\frac{e^{- {P e}_{1}} - 1}{S_{c}} (1 - {\hat{r}}_{n}^{* 2}) + 1 = \frac{D a}{4} (1 - {\hat{r}}_{n}^{* 2}) + \frac{D a}{2} {\hat{r}}_{n}^{* 2} \ln ({\hat{r}}_{n}^{*}),

(3.11)

whose derivation is detailed in the electronic supplementary material, Section R5. Let us assume that convection is sufficiently high that $c_{2}$ is approximately symmetric and that equation (3.11) may be used on $Ω$ .

The necrosis-prevention condition on $S_{c}$ is obtained from the same reasoning as for equation (3.7) at the static asymptote. Noting that $- 1 < e^{- {P e}_{1}} - 1 < 0$ , it suffices that

S_{c} (1 - D a / 4) > 1.

(3.12)

In terms of peak inlet fluid velocity, which is the most easily adjustable experimental parameter, this condition can be written

v_{0} > \frac{v_{d i f, 2}}{4 (\frac{1}{D a} - \frac{1}{4})} .

(3.13)

This depends on M-CELS properties only and not on diffusive transport in the microfluidic chamber. Like at the static asymptote, this condition expresses that necrosis prevention is penalized by fast diffusion in the M-CELS.

For an unconfined M-CELS of radius 0.5 mm where $D a = 3$ characterizes oxygen kinetics, the inlet-velocity condition is $v_{0} \sim 10^{- 5}$ m s⁻¹. The corresponding inlet flow rate for an M-CELS culture chamber of 1 cm $\times$ 1 mm cross-section (e.g. [17,18]) is 1 ml min⁻¹. Flow rates must be limited by the physiological shear stress of the cultured cells, which is under 10⁻³ Pa at rest outside of blood vessels [34]. In order not to affect cell function in an M-CELS, peak microfluidic velocity thus should remain under 10⁻⁴ m s⁻¹ (for a fluid viscosity of 1 mPa s and a gap of 100 between M-CELS apex and channel ceiling). That corresponds to 10 ml min⁻¹ with a 1 cm $\times$ 1 mm cross-section. Necrosis prevention via external nutrient convection over an unconfined M-CELS is therefore experimentally achievable.

Numerical simulations show that necrosis prevention is achieved if $S_{c} (1 - D a / 4) > 5$ in the unconfined geometry (figure 11a). The difference with the analytical condition (3.12) results from the diffusive boundary layer that grows from $Γ$ into $Ω_{1}$ , which nutrients need to cross to reach $Ω_{2}$ . In the confined geometry, nutrients also need to diffuse across the corner recirculation areas, so a much higher excess of convection is required to prevent necrosis. That excess increases with lower $R_{d}$ and higher $D a$ , such that $S_{c} \sim 10^{4}$ is insufficient to prevent necrosis at $D a = 3$ (figure 11b). This value is not advisable in experiments because it implies unphysiological shear stress on cells. For oxygen in a confined M-CELS of radius 0.5 mm, the convection limit due to shear stress is at $S_{c} \sim 100$ . This means that necrosis prevention by microfluidic flow is possible if $R_{d} \geq 0.5$ at $D a = 2$ , but impossible if insufficient diffusive supply means that it exists at $D a = 3$ .

Figure 11. Simulated live fraction in microfluidic culture for $D a$ such that any necrosis is supply-induced ( $D a < {D a}_{\infty}$ ). (a) *Unconfined geometry*: necrosis prevention achieved if $S_{c} (1 - D a / 4) > 5$ for $D a = 2$ and $D a = 3$ . (b) *Confined geometry*: no excess of convective supply suffices for general necrosis prevention. The required excess increases with a more deficient diffusive supply, i.e. with lower $R_{d}$ and higher $D a$ .

3.3.4. Transport asymmetry under moderately high convection

The asymmetry of solute supply to $Γ$ when ${P e}_{1} > 1$ and $S_{c} < 1$ translates into asymmetric necrotic cores, as the example of $D a = 5$ and $R_{d} = 0.5$ shows (figure 12). The shape of the cores follows the progressive supply increase to $Γ$ as convection in $Ω_{1}$ increases. The cores first shrink on the upstream half of $Ω_{2}$ ( $Γ \cap {x < 0}$ ), before shrinking along $y$ as the supply deficiency is filled in the constriction between $Γ$ and $W_{h}$ , and finally shrinking on the downstream half of $Ω_{2}$ . In the confined geometry, the third phase is followed by a reduction of the cores near the confinement corners. In that phase, ever-increasing convection in the primary flow region increases the intensity of the first recirculation eddy and therefore convective solute fluxes near the parts of $Γ$ exposed to the recirculation region (figure 12b,d).

To estimate the configuration of maximum asymmetry of solute concentration, we consider the quasi-one-dimensional domain $\hat{Ω}$ and derive the conditions where ${\hat{γ}}^{*} = 0.5$ (electronic supplementary material, Section R6). Those are

{D a}^{1 / 2} S_{c} = 4 \sqrt{2} .

(3.14)

The quantity ${D a}^{1 / 2} S_{c}$ is a good predictor of the maximum asymmetry of simulated necrotic cores (figure 13). Using dimensional parameters,

Figure 13. Trajectory of necrotic barycentres in microfluidic culture ( $x$ -coordinate). The value of ${D a}^{1 / 2} S_{c}$ determines an asymmetry peak. (a) Unconfined geometry and (b) confined.

{D a}^{1 / 2} S_{c} = 2 v_{0} \sqrt{\frac{c_{0}}{R_{m a x} D_{2}}} .

(3.15)

The symmetric nature of solute concentration thus only depends on solute kinetics inside M-CELS and on boundary conditions on $I$ . Since it does not depend on the M-CELS diameter or the microfluidic channel dimensions, the peak inlet velocity could be used as a single parameter to induce persistent asymmetric concentrations inside M-CELS and investigate directional responses. The maximum downstream shift of the simulated necrotic barycentres occurs for ${D a}^{1 / 2} S_{c} \in [1,4]$ in the unconfined geometry (figure 13a), showing that ${D a}^{1 / 2} S_{c}$ is a good predictor of supply asymmetry. The shift was delayed and less marked in the confined geometry because solute supply increases more slowly near the recirculation areas: its maximum occurs for ${D a}^{1 / 2} S_{c} \in [3,10]$ (figure 13b).

3.3.5. Approach of maximal-supply asymptote

Let us consider cases of rate-induced necrosis ( $D a > {D a}_{\infty}$ ) and assume that convection is sufficiently high that the asymmetry peak has passed. We consider this to be ${D a}^{1 / 2} S_{c} > 8$ in the unconfined geometry and ${D a}^{1 / 2} S_{c} > 20$ in the confined geometry, i.e. twice the asymmetry peak. The approach of the maximal-supply asymptote depends on the rate at which solute supply increases on $Γ$ via the diffusion-dominated areas. Let us consider the recirculation areas and define a corner Péclet number ${P e}_{ρ}$ based on the recirculation properties and the quasi-one-dimensional no-flow solution,

{P e}_{ρ} = \frac{γ v_{ρ, 1}}{D_{1} {\frac{d c}{d \hat{x}} |}_{\hat{Γ}}},

(3.16)

where $v_{ρ, 1}$ is the maximum velocity of the largest eddy, with $v_{ρ, 1} = q v_{0}$ and $q \sim 10^{- 3}$ (fig. 5 in [33]). From the no-flow expression of concentration in $\hat{Ω}$ , obtained by substituting electronic supplementary material equations (22) into (20), we obtain

{P e}_{ρ} = q {P e}_{1} (\frac{S_{d}}{1 - {\hat{r}}_{n, 0}^{* 2}} - 1) .

(3.17)

Substituting ${\hat{r}}_{n, 0}^{*}$ with $δ_{0}^{*}$ using electronic supplementary material, equation (29) and expanding at order 1 in $δ_{0}^{*}$ gives

{P e}_{ρ} = q {P e}_{1} (\frac{R_{d} S_{d}}{2 (\sqrt{1 + R_{d} S_{d}}) - 1} - 1) .

(3.18)

Under deficient diffusive supply ( $R_{d} S_{d} ≪ 1$ ), successive Taylor expansions lead to

{P e}_{ρ} = q R_{d} S_{c} .

(3.19)

The balance between convective supply and consumption on $Γ_{ρ}$ is defined, analogously to $S_{c}$ , as

\begin{aligned} S_{c, ρ} & = \frac{2 R_{d} {P e}_{ρ}}{D a} \\ = q R_{d} S_{d} S_{c} . \end{aligned}

(3.20)

A similar reasoning may be followed for the diffusive boundary layer covering all of $Γ$ . The velocity of the first eddy $q v_{0}$ would be replaced by the velocity at the entrance of the boundary layer, which also scales with $v_{0}$ . Solute supply to $Γ$ thus is governed by $R_{d} S_{d} S_{c}$ when $S_{c} ≫ 1$ . The excess of convective supply required to maintain a given proportion of live cells ( $Φ_{L}$ ) scales with $(R_{d} S_{d})^{- 1}$ . The corresponding peak inlet fluid velocity scales like

v_{0} \sim \frac{v_{d i f, 2}}{S_{d}^{2}} .

(3.21)

The approach of maximal supply via external convection thus is delayed by deficient diffusive supply to the M-CELS and by fast diffusion within it. Using dimensional parameters,

v_{0} \sim a D_{2} {(\frac{R_{m a x}}{c_{0}})}^{2} {(\frac{L}{D_{1}})}^{2} .

(3.22)

This velocity evolves linearly with the M-CELS diameter, so it would need to increase with M-CELS age if solute penetration distance were to be maintained for an experiment. In the example of neurons in a cortical organoid, the increase in radial diffusivity resulting from neuronal alignment would also contribute to $v_{0}$ needing to increase with M-CELS age. The dependency of $v_{0}$ on the square of the inverse of the inlet concentration and the square of the distance to the solute source provides a means to reduce experimentally required flow rates. These come with the caveats that inlet concentrations are often physically or biologically constrained (e.g. dissolved gases like oxygen have a saturation limit) and that the distance to the solute source is constrained.

The numerical simulations verify that the rate of approach of maximal supply is governed by $R_{d} S_{d} S_{c}$ . In the unconfined geometry, overcoming the diffusive boundary layer on $Γ$ to achieve $Φ_{L} / Φ_{L, \infty} = 0.9$ requires $R_{d} S_{d} S_{c} > 10$ (figure 14a). In the confined geometry, the additional diffusive barrier constituted by the corner recirculation areas make that condition $R_{d} S_{d} S_{c} > 10^{3}$ . This limits the biological configurations where supply deficiency can be filled within the constraints of shear stress on cells.

3.4. Comparison with previous models and experiments

The results presented above shed new light on some previously published simulations and experiments (table 2). Let us note that the predicted onset of necrosis in the unconfined geometry ( ${D a}_{\infty} = 4$ ) is in agreement with calculations in cylinders in [20] and [29] and that spheres are characterized by ${D a}_{\infty} = 6$ , as these studies demonstrate. Similarly, the solute quantity balances on $Γ$ are defined as $S_{d} = 3 R_{d} / D a$ and $S_{c} = 3 R_{d} {P e}_{1} / D a$ for spheres.

Table 2. Application of results to previous studies. Subscripts 1 and 2 are applied to the fluid or gel outside the M-CELS of interest and to the M-CELS, respectively (as in the models). Details of the values used for the comparisons are in the electronic supplementary material, table S3.
Collapse
M-CELS type	confining geometry	study results^a
human midbrain organoid [15]	adhesion to bottom of reservoir and encapsulation in GelTrex	10% reduction in necrotic area under continuous fluid flow vs. a shaking platform (from 16% to 14%).
Application: $D a = 0.5$ for oxygen in the 400 µm diameter spherical organoids. The 1 mm thick gel capsule imposes a diffusion-limited zone which can be represented by $S_{d} = 0.75$ if fluid flow is large enough that $c_{1}^{*} = 1$ outside of it. $Φ_{L}$ is very sensitive to $S_{d}$ in that region, all the more so as $D a$ is low, so it is coherent that necrosis would be large enough to stop growth for such a value of $S_{d}$ (figure 7b). Experimentally, $Φ_{L} = 0.85$ . The model (not simulated for $D a = 0.5$ ) would slightly underestimate that as there are quiescent cells near the necrotic area that consume less oxygen and thus induce deeper oxygen penetration than the model predicts.
human retinal organoid [35]	ellipsoidal well, separated from fluid flow by a membrane	fluid flow increases the simulated penetration of oxygen and divides the experimental number of dead cells (Caspase3+) by 2.5.
Application: $D a = 4$ for oxygen in the spherical 400 µm diameter organoids. In static conditions, $S_{d} = 0.3$ and the model predicts^b $Φ_{L, 0} = 0.4$ . The distance between the organoid and the membrane implies that $S_{d} = 1.5$ at best if fluid flows fast enough in the upper channel that $c_{1}^{*} = 1$ there^c. This would leave $Φ_{L} = 0.9$ under fluidic conditions and there would indeed be some Caspase3+ cells left. The given Caspase3+ ratio rather suggests $Φ_{L} = 0.75$ , which corresponds to figure 11b with the experimental^d $D a = 4$ , $R_{d} = 0.4$ , and $P e = 20$ .
HepG2 tumour spheroid [36]	cylindrical well, separated from fluid flow by membrane	spheroids grow to a maximum diameter of 200 µm in static culture and to at least 300 µm in fluidic culture
Application: the maximum diameter corresponds to $D a = 1.6$ and $S_{d} = 0.13$ for oxygen, for which the model predicts $Φ_{L, 0} = 0.2$ (figure 7b)^e. Quiescent cells make this an underestimation, but it does indicate a large necrotic area which could be sufficient to stop growth [14]. In fluidic culture, $S_{c} \sim 10^{2}$ which is ample enough to avoid necrosis until spheroid diameters of at least 300 µm, as the authors observe.

^aOnly results that are traceable to the supply of a solute to M-CELS. Cell activity markers are therefore not considered here.

^bWe assume that a sphere and a cylinder behave identically if they have the same Da/Da∞ and we follow the Da = 3 curve on figure 7.

^cS_d is calculated from the chip inlet in the static case and from the membrane in that hypothetical case where the upper channel is maximally supplied.

^dDa = 3 on the figure, see note b.

^eDa = 1 on the figure, see note b.

4. Discussion

Multi-cellular engineered living systems (M-CELS) are situated at a certain distance from the source of soluble nutrients, growth factors or drugs that need to be supplied to them. Since solute transport is diffusion-dominated within M-CELS, their internal solute distribution is linked to their external solute transport via the solute concentration at the fluid–M-CELS interface only. It thus is key to rapidly cross the source–M-CELS distance. Both slow external transport relative to M-CELS metabolic needs and a large surface area of contact with the culture-chamber walls limit solute penetration into M-CELS.

We quantified the maximum consumption-to-diffusion kinetics ratio (noted ${D a}_{\infty}$ ) that an M-CELS can reach without necrosis as a function of its structural confinement. Finite transport kinetics in microfluidic chambers imply the need for an excess of solute supply to the M-CELS surface. That excess is a hyperbolic function of the ratio $D a / {D a}_{\infty}$ . The cost of preventing necrosis thus diverges when $D a$ approaches ${D a}_{\infty}$ , e.g. requiring infinite flow rates. This renders the realization of a non-necrotic M-CELS at ${D a}_{\infty}$ unfeasible and places a lower limit on the maximum consumption-to-diffusion ratio that can be achieved in a given experimental set-up. That limit depends on the design constraints of the set-up and cellular sensitivity to shear stress.

External convection brings about a higher solute supply to the fluid–M-CELS interface only after the solutes have crossed zones where their transport is diffusion dominated. These cover the whole fluid–M-CELS interface and include the diffusive boundary layer that grows from the interface into the fluidic chamber and the recirculation area near each confinement corner. That area exists at the intersection of the M-CELS and the wall in the cylindrical geometry we simulated. It would be present near any corner of angle inferior to 146° in the confinement structure of a spherical M-CELS [33]. Slow diffusion across those areas relative to solute consumption in the M-CELS limits the benefits of external convection. This arises when the diffusion-dominated areas are thick, e.g. when there are narrow corners, and when the diffusive supply in an equivalent static configuration is deficient, e.g. when the distance between the solute source and M-CELS is large. In both cases, slow diffusion kinetics in the culture chamber lower concentration at the fluid–M-CELS interface, and external convection needs to be all the higher to overcome the diffusive barrier. The same reasoning applies to M-CELS encapsulated in hydrogels. Solute diffusivity across hydrogels is close to their diffusivity in water because hydrogel porosities are high, but capsule thickness is detrimental to increased solute penetration via external convection. One should thus aim for capsules to be as thin as possible so as not to impede solute supply while maintaining the cell–matrix interactions and the protection against shear stress that these provide.

The easiest way to increase solute supply to fluid–M-CELS interfaces is to modulate transport velocities, either by design for diffusion (short distance between solute source and M-CELS) or by operation for convection (peak fluid velocity). The effect of both is penalized by fast diffusion within the M-CELS, independently of how the internal-diffusion kinetics compare with consumption and external-diffusion kinetics. Fast internal diffusion means that any increase in solute supply to the fluid–M-CELS interface is more rapidly distributed throughout the M-CELS, slowing down the increase in interface concentration that is necessary to increase solute penetration depth. Achieving non-necrotic M-CELS (or a given penetration depth) at a given $D a$ thus requires less external convection in large M-CELS with low consumption than in small M-CELS with high consumption.

The fundamental understanding of how microfluidic chamber design and operation affect solute supply to M-CELS that we have presented allows us to state how much larger rigid avascular M-CELS can grow by adding fluid flow to the culture set-up. In the oxygen example that illustrated our results, the maximum no-flow $D a$ was a quarter of ${D a}_{\infty}$ , i.e. the maximum M-CELS radius was divided by 2. Fluid flow around an unconfined M-CELS allows to recover a maximal $D a$ close to ${D a}_{\infty}$ (figure 11a). Peak fluid velocity is constrained in experimental set-ups because of the shear-stress sensitivity of cells, which capped the maximum achievable $D a$ at ${D a}_{\infty} / 2$ in the presence of recirculation areas, i.e. an increase in maximal radius of around 40% compared with static culture (figure 11b). Similarly, approaching maximal supply in an M-CELS with incomplete solute penetration required a 100 times faster fluid flow because of the need for solutes to cross the recirculation areas (figure 14b). This highlights the importance of eliminating thick diffusion-dominated areas in confinement structures. These exist in widely used U-shaped barriers and wells (e.g. [37]), but are absent from designs such as slit U-shaped wells [38] or porous confinement scaffolds [39]. These results also point to solute-transport benefits being inhibited by the encapsulation of M-CELS in hydrogels.

Real M-CELS culture systems contain three additional complexities that our models do not explicitly address. First, spheroids and many organoids (e.g. brain, liver) are spherical or close to it. That implies a larger ${D a}_{\infty}$ (6 for a sphere [20,29]), but the reduction in ${D a}_{\infty}$ by confinement and the shape of necrotic cores would be similar. In culture, spherical M-CELS receive solute from the third dimension in the same ways as from the one orthogonal to the inlet–outlet axis ( $y$ ) in our two-dimensional models. Those are diffusion, whose fluxes are reduced by constriction, and convection from shear flow. The necrosis prevention conditions would therefore be the same within constants. Maximal-supply approach rates would be the same. Second, spherical M-CELS may be confined without recirculation areas, such as behind U-shaped traps with a slit [38]. This is an intermediary situation between the two studied here. The static and maximal-supply asymptotes would follow the same laws as in the confined geometry, but necrosis prevention and the approach of maximal supply by convection follow the unconfined geometry because they are not delayed by diffusion-dominated areas. Third, M-CELS may be cultured in arrays to increase experimental throughput, which is particularly beneficial on drug-testing platforms [8]. This format is detrimental to solute supply in static culture because of the increased distance to the solute source with each row of the M-CELS array. In fluidic culture, it suffices that the M-CELS are sufficiently far apart that solute mixing in the wake of the upstream M-CELS restores solute concentration in the fluid to its source value. We analyse concentration near an M-CELS in electronic supplementary material, Section R7 and show that the necrosis-prevention condition is transferable to an array of spheroids if the distance between spheroids is a tenth of the entry length of the culture channel, so typically of the order of 1 mm. This is more than the few hundred micrometres that often separate spheroids cultured in arrays [36,40,41]. Developing a predictive model of solute supply to arrays of M-CELS that incorporates the array dimensions would contribute to increasing the repeatability of assays performed on M-CELS. Biologically, the assumptions of uniform properties hold for tumour spheroids and early-stage spherical organoids. Later-stage ‘rough’ organoids, such as wrinkled brain organoids [32,42], do not because the troughs in their surface are diffusion highways that enable deeper penetration of solutes, particularly of fast-diffusing ones like oxygen and glucose. Organoids with differentiated cell layers, e.g. an outer layer of high-consumption neurons, would present a different maximal-supply asymptote [20], but the analysis of external supply would still apply to them. Current experimental efforts in organoid research focus on incorporating a perfusable vasculature to better mimic in vivo mechanical forces and solute-transport patterns. This study of the link between the geometry and operation of the culture set-up and solute concentration inside M-CELS can inform the boundary conditions at the entrance of the vessels. This would, in turn, contribute to quantifying the inhomogeneity of solute transport in perfusable M-CELS as well.

Ethics

This work did not require ethical approval from a human subject or animal welfare committee.

Data accessibility

The simulation results and processing codes are on the Zenodo repository [43].

Supplementary material is available online [44].

Declaration of AI use

We have not used AI-assisted technologies in creating this article.

Authors’ contributions

W.V.B.: conceptualization, data curation, formal analysis, investigation, methodology, writing—original draft, writing—review and editing; N.K.: investigation, methodology, visualization; M.T.: conceptualization, funding acquisition, supervision, writing—review and editing; S.B.: conceptualization, formal analysis, funding acquisition, investigation, methodology, supervision, writing—review and editing.

All authors gave final approval for publication and agreed to be held accountable for the work performed therein.

Conflict of interest declaration

We declare we have no competing interests.

Funding

This project has received funding from the Olle Engkvist Foundation (OES 213-0231), the Knut and Alice Wallenberg Foundation (KAW 2021.0172) and the European Research Council under the European Union’s Horizon Europe research and innovation programme (PHOENIX grant agreement 101043985) via grants led by MT.

Acknowledgements

All authors gratefully acknowledge support from T. Christian Gasser. WB and SB gratefully acknowledge support from the Department of Engineering Mechanics, School of Engineering Sciences, KTH.

Footnotes

Electronic supplementary material is available online at https://doi.org/10.6084/m9.figshare.c.7634294.

^†Present address: Institut de Physique de Rennes, Université de Rennes, Rennes, France

Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.

VIEW FULL TEXT DOWNLOAD PDF

This Issue

Journal of The Royal Society Interface cover image

January 2025

Volume 22Issue 222

Article Information

DOI:https://doi.org/10.1098/rsif.2024.0463
PubMed:39875093
Published by:Royal Society
Online ISSN:1742-5662

History:

Manuscript received08/07/2024
Manuscript accepted05/11/2024
Published online29/01/2025

License:

Citations and impact

Keywords

Subjects

Your password must have 8 characters or more and contain 3 of the following:

Capacity and limitations of microfluidic flow to increase solute transport in three-dimensional cell cultures

Abstract

1. Introduction

2. Methods

2.1. Computational domains

2.2. Quasi-one-dimensional analytical domain

2.3. Equations for fluid flow

2.4. Equations for nutrient transport

2.5. Nutrient balance at the fluid–multi-cellular engineered living systems interface

2.6. Numerical methods

2.7. Measures of solute transport

2.8. Definition of asymptotes

3. Results

3.1. Maximal-supply asymptote

3.2. Static asymptote

3.2.1. Prevention of necrosis

3.2.2. Approach of maximal-supply asymptote

3.3. Microfluidic culture

3.3.1. Function of increased convective supply

3.3.2. Evolution of solute supply with increasing convection

3.3.3. Prevention of necrosis

3.3.4. Transport asymmetry under moderately high convection

3.3.5. Approach of maximal-supply asymptote

3.4. Comparison with previous models and experiments

4. Discussion

Ethics

Data accessibility

Declaration of AI use

Authors’ contributions

Conflict of interest declaration

Funding

Acknowledgements

Footnotes

This Issue