## Abstract

The placenta is an essential component of the life-support system for the developing foetus, enabling nutrients and waste to be exchanged between the foetal and maternal circulations. Maternal blood flows between the densely packed branches of villous trees, within which are foetal vessels. Here, we explore some of the challenges in modelling maternal haemodynamic transport using homogenization approaches. We first show how two measures can be used to estimate the minimum distance over which the distribution of villous branches appears statistically homogeneous. We then analyse a simplified model problem (solute transport by a unidirectional flow past a distribution of point sinks) to assess the accuracy of homogenization approximations as a function of governing parameters (Péclet and Damköhler numbers) and the statistical properties of the sink distribution. The difference between the leading-order homogenization approximation and the exact solute distribution is characterized by large spatial gradients at the scale of individual villi and substantial fluctuations that can be correlated over lengthscales comparable to the whole domain. This study highlights the importance of quantifying errors owing to spatial disorder in multi-scale approximations of physiological systems.

## 1. Introduction

The placenta is a vital organ for the foetus. It performs multiple simultaneous functions, including exchange of blood gases between the maternal and foetal circulations and supply of nutrients to, and removal of waste products from, the foetus. In humans, where the embryo implants within the interstitium of the uterus, the placenta is highly invasive and haemochorial in organization: in other words, maternal blood emerging from spiral arteries in the uterine wall flows freely around chorionic villous trees containing the foetal vessels [1,2] (figure 1*a*). This evolutionary strategy is shared between humans and only a few other species. The close maternal apposition makes human placental development and functioning vulnerable to maternal disease. For example, in pre-eclampsia (associated with intra-uterine growth restriction), inadequate invasion, raised maternal blood pressure and restricted maternal blood flow into the placenta are thought to impair development of foetal blood vessels and lead to reduced chorionic villous branching [3]. In maternal diabetes, there is increased angiogenesis and increased villous tree branching [4]. In either case, normal transport processes are compromised. Damage to the foetus during pregnancy can have immediate health implications as well as legacy effects which manifest in adulthood [5]. There are therefore obvious benefits in developing systems-level computational models of placental flow and transport that will help explain the origins and implications of structural abnormalities in disease, and help reveal the interactions between geometry, blood flow, growth and nutrient uptake that are central to the function of this organ.

Despite its relative accessibility (using a freshly delivered term-placental perfusion system) and its medical importance, theoretical models of placental function (summarized in [6]) remain relatively primitive. A particular challenge is the intricate geometrical structure of the villous trees (figure 1*b*), which is hard to characterize, highly variable and prohibitively complex to simulate in detail. A natural starting point in modelling maternal placental blood flow is to exploit the theory of flow in porous media [7,8]. We have recently used this approach to demonstrate the importance of the calibre of maternal spiral arteries and decidual veins in determining the overall flow resistance in a functional placental unit (or ‘placentone’, containing a single villous tree) and to examine the factors that optimize nutrient delivery [6].

Behind porous medium models lies the theory of homogenization [9–11], which provides a systematic tool for characterizing the bulk properties of materials in terms of their fine-scale structure. This method is particularly well developed for media with periodic microstructure. It is used for example to derive Darcy's law for viscous flow in a porous medium [12] and its extensions that account for inertia [13] and poroelasticity [14]. In physiology, Darcy's law is widely used to simulate interstitial flows; it has also recently been applied to describe homogenized flow and transport within microcirculatory networks [15] when leakiness allows communication between nearby vessels. Using volume-averaging procedures, approximations have also been developed for materials with random microstructure [16]. For materials with statistically homogeneous and ergodic microstructure, virtually all leading-order results for periodic media are directly applicable [12], including Darcy's law [17]; rigorous bounds on effective material properties for composite materials are available [12]. However, fewer results are available connecting the statistical properties of the microstructure to the properties of the homogenization residue (the difference between the leading-order approximation and the exact solution) [18,19].

In order to understand how the placenta functions as an organ of nutrient exchange, we here examine the combined effects of transport of solutes such as glucose or oxygen by maternal blood flow and uptake by the foetal circulation in villous trees. We can therefore call on prior studies [20–23] of the competition between advection, diffusion and reaction in porous media or heat exchangers, characterized by a Péclet number *Pe* (relating the strength of advection to diffusion) and by a Damköhler number *Da* (relating the rate of reaction or uptake to diffusion). In addition to physical parameters, it is equally important to characterize the geometry of villous branches. Careful placental stereology (using systematic random sampling at the microscale) has provided estimates of bulk quantities (e.g. total villous volumes, surface areas and lengths) and local measures such as star volumes (the mean volume of all parts of a space which are visible when viewed in all directions from a given point within it) [24]. However, further measures of placental anatomy, particularly of its statistical variability, are necessary in order to develop comprehensive models of placental transport.

In this paper, we draw together a number of techniques and ideas that are useful in developing multi-scale stochastic models of nutrient transport in the human placenta. Using histological images (figure 1*b*), we illustrate in §2*a* how methods of spatial statistics [25] can be used to characterize some of the important underlying lengthscales in villous trees. We then use an idealized theoretical model (transport under flow past a one-dimensional array of sinks, §2*b*) to assess how governing parameters (*Pe*, *Da*) and the statistical properties of the sink distribution together affect the accuracy of homogenization approximations (§2*c,d*). We show in particular how randomness leads to large fluctuations in solute distributions that can be correlated over surprisingly long distances.

## 2. Methods

### (a) Imaging and spatial statistics

A peripheral lobule of the human placenta was used for the histological study, as it tends to contain a single chorionic villous tree matched with a spiral artery, and is therefore more representative of a single placental circulatory unit [2,6]. A single maternal lobule, obtained from a normal full-term placenta (delivered by elective caesarian section), was frozen in nitrogen-cooled isopentane before cutting according to a systematic random sampling protocol. Sections (8 μm thick) were stained with tolonium chloride (toluidine blue), revealing the chorionic villi and the medium-to-large blood vessels within them. An example is shown in figure 1*b*, using a QImaging MicroPublisher 5 Mpx camera and a Zeiss Axioplan microscope. We used ImageJ to identify villous branch outlines (figure 2*a*) and compute their area fraction *ϕ* as a function of window size *W* [26]. Centres of mass of each branch (figure 2*b*) enable us to interpret villous branch locations as a spatial point process to which we can apply appropriate statistical methodology [25].

Let *N*(*A*) represent the number of points lying within an area *A*. The first- and second-order intensities of this process are
2.1
in the limit , where d*A*(**x**) is an infinitesimal area around a given point **x**. The density *λ* is the expected number of points per unit area. The two-point correlation function *λ*_{2} is the normalized joint number of points expected around two given points in space. For a stationary isotropic point process, for which *λ*=constant, *λ*_{2}(**x**_{1},**x**_{2})=*λ*_{2}(|**x**_{1}−**x**_{2}|) [12]. The *K*-function [25] is the average number of points within a given distance *r* from an arbitrary point **x**_{0}, scaled with *λ*, where
2.2
For a completely random (Poisson) spatial point process, *K*=*πr*^{2} and *λ*_{2}=*λ*^{2}. For *n* points in a given rectangular domain of area |*A*| and perimeter *P* (as in figure 2*b*), we employ Ripley's estimate for the *K*-function defined as
2.3
where I is an indicator function and *w*_{ij} is Ripley's trigonometric weighting function to correct for edge effects [25].

We use the *K*-function to estimate the cross-correlation of intervillous distances and to test the pattern for regularity or clustering, when compared with the complete spatial randomness represented by a Poisson point process [25]. For the latter, we use
2.4
where *a*_{1}(*r*)=|*A*|^{−2}(0.21*Pr*^{3}+1.3*r*^{4}), *a*_{2}(*r*)=|*A*|^{−3}(0.24*Pr*^{5}+2.62*r*^{6}) and *b*(*r*)=*πr*^{2}|*A*|^{−1}(1−*πr*^{2}|*A*|^{−1})+|*A*|^{−2}(1.0716*Pr*^{3}+2.2375*r*^{4}) are Lotwick and Silverman's polynomials [27]. We estimate in equations (2.3) and (2.4) using the Kest function of the Spatstat package for `R` [28]; equation (2.4*b*) is exact for *r* not exceeding 1/4 of the smallest side of a rectangular domain [27]. Measurements of area fraction *ϕ* and the *K*-function for placental tissue are discussed in §3*a*.

### (b) A model problem for solute uptake

As a simple representation of solute uptake by a villous tree, such as within lobule L_{1} in figure 1*a*, we consider a one-dimensional array of *N* identical point sinks of strength *q*_{0}, distributed across a domain 0≤*x**≤*L*. The distance *L* represents a typical pathlength between a spiral artery and decidual vein, passing multiple villous branches; see [6]. The typical distance between two adjacent sinks is *l*=*L*/(*N*+1) and we assume *ε*≡*l*/*L*≪1. A solute with steady concentration distribution *C**(*x**) is swept past the sinks by a uniform flow field *u*_{0} and diffuses between them with diffusion coefficient *D*; the concentration *C*_{0} at the inlet (representing a spiral artery) is prescribed and the concentration at the outlet is set to be zero. The concentration field is assumed non-negative, so that, if *q*_{0} is sufficiently large, an internal free boundary at *x**=*x*^{*}_{0} arises such that *C**>0 for 0<*x**<*x*^{*}_{0}, and *C**=0 for *x*^{*}_{0}≤*x**≤*L*. Writing *C**(*x**)=*C*_{0}*C*(*x*), *x**=*lx*, , the dimensionless solute distribution satisfies
2.5a
2.5b
2.5c
where *Pe*=*u*_{0}*l*/*D* is a Péclet number and *Da*=*q*_{0}*l*/(*DC*_{0}) is a Damköhler number, each defined using the inter-sink distance. Solute uptake is assumed to follow zeroth-order kinetics, the sinks being represented by
2.6
for the ordered sink locations 0<*ξ*_{1}<*ξ*_{2}<⋯<*ξ*_{N}<*ε*^{1}. For convenience, we define *ξ*_{0}=0 and *ξ*_{N+1}=*ε*^{1} and define Δ_{i}≡*ξ*_{i}−*ξ*_{i−1} (*i*=1,…,*N*+1).

Motivated by the spatial patterns in figure 2*b*, we consider four sink distributions, allowing us to mimic structural features that are potentially of relevance to villous trees:

*f*=*f*_{p}; say, a periodic distribution with*ξ*_{i}=*i*, Δ_{i}=1 for 1≤*i*≤*N*;*f*=*f*_{u}; say, a uniformly random distribution, where*ξ*_{i}are ordered values drawn from independent distributions, implying*εξ*_{i}∼Beta(*i*,*N*+1−*i*) and*ε*(Δ_{1},…,Δ_{N})∼Dirichlet(|*ξ*), where*α*=(*ξ**ξ*_{1},…,*ξ*_{N}) and [29];*f*=*f*_{h}(*d*); say, a Matérn hard-core type II distribution, whereby sinks are drawn sequentially from a uniformly random distribution and are accepted provided they do not fall within a distance*d*of an existing sink or boundary, the process being continued until*N*sinks are reached [25]; and*f*=*f*_{n}(*σ*); say, a normally perturbed periodic distribution satisfying , for some variance*σ*^{2}, with periodic conditions imposed on sinks falling outside the domain. Sinks are reordered if they swap positions (as occurs for sufficiently large*σ*).

In the limit , *f*_{h} tends to *f*_{u}; in the limit , *f*_{h} exhibits increasing regularity, where *d*_{cr} is an upper bound estimated empirically to be close to [30]. All distributions have the first-order intensity *λ*=1. The population autocorrelation of the inter-sink distance, defined by
2.7
satisfies for *n*>0 and small *σ* for *f*=*f*_{n} (reflecting independent sink locations) and (using standard properties of Dirichlet distributions [31]) for *n*>0 for *f*=*f*_{u} (i.e. a large gap between sinks is compensated by reduced gaps elsewhere); simulations (taking the mean from an ensemble of realizations of each process) confirm that for all three random distributions (when *d*=0.65, *σ*=10; data not shown), implying that there is no substantial long-scale correlation in the sink distributions.

We simulate realizations of each distribution and compute the transport problem (2.5) over each distribution using a semi-analytical method, matching concentration fluxes at each sink. Mesh convergence was checked by validation against the results obtained by the finite-element solver COMSOL Multiphysics and exact solutions in simple cases. Results are presented in §3*b,c*, tested against predictions from homogenization approximations that we now develop for periodic and random sink distributions.

### (c) Homogenization approximation: periodic sinks

We develop asymptotic approximations of equation (2.5) in the limit using homogenization. When *f*=*f*_{p}, a standard approach may be adopted, described in appendix A, the outcomes of which we summarize briefly here. The concentration field is represented by a two-scale expansion of the form *C*(*x*)=*C*^{(0)}(*X*)+*εC*^{(1)}(*x*,*X*)+*ε*^{2}*C*^{(2)}(*x*,*X*)+⋯ , where *X*=*εx* characterizes the slow variation of the solution over lengthscales comparable to the whole domain, while *x* characterizes rapid variations over a ‘unit cell’ surrounding an individual sink.

Key features of the concentration field over a periodic distribution of sinks are summarized in figure 3. Diffusive, advective and uptake fluxes across the whole domain balance when *DC*_{0}/*L*∼*u*_{0}*C*_{0}∼*q*_{0}*N*, which corresponds to *Da*=*O*(*ε*^{2}), *Pe*=*O*(*ε*). For given *ε*, this defines an organizing centre in (*Pe*, *Da*)-parameter space (figure 3*a*). A second organizing centre lies at *Pe*=*O*(1), *Da*=*O*(*ε*). By deriving distinct asymptotic limits for *C* around each organizing centre, we can obtain a comprehensive overview of the homogenization approximation across parameter space. This provides a useful perspective before focusing on parameters appropriate for specific solutes.

We first set *Pe*=*εp*, *Da*=*ε*^{2}*q* and assume *p*,*q*=*O*(1) as *ε*→0. *C*^{(0)} then satisfies (see equation (A5)). The two-parameter solution for *C*^{(0)} (given by equation (A6)) encompasses four distinct limits, denoted by U_{A}, U_{D}, D and A in figure 3*a*, in which different physical effects dominate. Representative solutions, illustrating the form of *C*^{(0)} in each region of (*Pe*, *Da*)-parameter space, are shown in figure 3*b* at points (1), (2), (3), (4), respectively. Briefly:

— in region D (

*p*≪1,*q*≪1, illustrated by point (3)),*C*^{(0)}≈1−*X*, representing diffusion between the boundary source and sink, the distributed sinks and advection being too weak to have any effect at the leading order;— in region U

_{D}(, illustrated by point (2)),*C*^{(0)}≈(1−(*X*/*X*_{0}))^{2}for and*C*^{(0)}=0 otherwise, representing uptake by the distributed sinks balancing diffusion;— in region U

_{A}(1≪*p*≪*ε*^{−1},*p*≪*q*≪*p*^{2}, illustrated by point (1)),*C*^{(0)}≈1−(*X*/*X*_{0})−(1/*pX*_{0})*e*^{−p(X0−X)}for 0≤*X*≤*X*_{0}, with*ε*≪*X*_{0}≈*p*/*q*≪1, representing uptake by the distributed sinks balancing advection (outside a diffusive boundary layer near*X*=*X*_{0}); and— in region A (1≪

*p*≪*ε*^{−1},*q*≪*p*, illustrated by point (4)),*C*^{(0)}≈1−e^{−p(1−X)}, representing dominant advection, supplemented by diffusion close to the outlet.

In the neighbourhood of *Pe*=*O*(*ε*), *Da*=*O*(*ε*^{2}), the first correction *εC*^{(1)} vanishes, whereas the second correction *ε*^{2}*C*^{(2)} (see equation (A8)) is proportional to *Da*. This rises to *O*(*ε*) as *Da* approaches *O*(*ε*), violating the assumed structure of the original asymptotic expansion, and indicating breakdown of the homogenization approximation along *Pe*≪1, *Da*=*O*(1). Furthermore, the correction within the unit cell develops an internal boundary layer (in *x*) of width 1/*Pe* adjacent to the sink as *Pe* increases through unity (compare points (1) and (5), figure 3*c*), signalling the development of ‘staircases’ in the concentration profile, occurring in regions A^{s} and U^{s}_{A} in (*Pe*,*Da*)-space (figure 3*a,b*). To capture these, a new expansion must be constructed near the second organizing centre at *Pe*=*O*(1),*Da*=*O*(*ε*). Writing *Da*=*εq*_{1} with *q*_{1}=*O*(1) as , we find that here *C*^{(0)} decreases linearly from the inlet (see equation (A10)), before terminating in a diffusive boundary layer either upstream of, or at, the downstream boundary (illustrated, respectively, by points (1) and (4) in figure 3*b*); the correction *C*^{(1)} satisfies equation (A12).

To assess the difference between the leading-order homogenized and exact solutions to equation (2.5), we define the homogenization residue (or corrector) as
2.8
For periodic sink distributions, we estimate the residue using *r*^{ε}≈*εC*^{(1)} for *Pe*=*O*(1) (and *r*^{ε}≈*ε*^{2}*C*^{(2)} for *Pe*=*O*(*ε*)) and assess its magnitude under the (weak) mean-squared and (strong) Sobolev norms,
2.9
The convergence of *C*^{(0)} to *C* across parameter space is discussed in §3*b* .

### (d) Homogenization approximation: random sink distributions

For the random sink distributions (*f*=*f*_{u}, *f*_{h} or *f*_{n} in equation (2.5)), the same leading-order problem (A5) emerges for *C*^{(0)}. This is readily shown via averaging over long spatial scales, following Holmes [32]. Thus, we can again make direct use of the parameter-space map in figure 3*a*.

We use two approaches to determine the size and the statistical properties of the residue (2.8). When *Pe*=*O*(*ε*) and *Da*=*O*(*ε*^{2}) (and diffusion dominates advection at the microscale), we can determine the statistics of *r*^{ε} directly (see appendix B). The fluctuations *r*^{ε} are related to linear combinations of the sink locations *ξ*_{i}, weighted by the slowly varying coordinate *X*. Evaluation of (see equation (B9)) yields the homogenized ordinary differential equation (A5*a*), while we find that the point-wise variance satisfies
2.10
(The equivalent expression for *f*=*f*_{h} is harder to obtain and is left for further work.) The variance varies smoothly with *X*, vanishing at either end of the domain where the concentration is prescribed. When the sinks are almost periodic (*f*=*f*_{n} in equation (2.10)), fluctuations about *C*^{(0)} are of magnitude *O*(*ε*^{3/2}) (in the *L*_{2} norm), whereas these rise substantially (to *O*(*ε*^{1/2})) when there is greater disorder in the sink locations (*f*=*f*_{u}).

For larger *Pe*, we compute *r*^{ε} for *f*=*f*_{u}, *f*_{h} and *f*_{n} using Monte Carlo simulations of equation (2.5). We estimate its convergence in mean by evaluating integrals with trapezium quadrature. We also estimate Var(*r*^{ε}) and the transverse covariance
2.11
which characterizes the degree to which fluctuations are correlated across the domain. Results are presented in §3*c*.

## 3. Results

### (a) Imaging

Figure 2*c* shows how the villous area fraction *ϕ* depends on the size *W* of the window used to compute it. When *W* is comparable with the diameter of terminal villi (below 100 μm), we see vigorous oscillations of *ϕ*, depending on whether the window falls in villous tissue or in intervillous space. As *W* increases, variations in *ϕ* fall and are significantly reduced for . Above this threshold, we can reasonably treat the intervillous space in this sample as a continuous medium of uniform (or at least slowly varying) area fraction, as was assumed for example, by Chernyavsky *et al.* [6].

In figure 2*d*, the estimated *K*-function (2.3) for the distribution of centres of villous branches is compared with that expected of a homogeneous Poisson process (with confidence intervals). The two are indistinguishable, for the given sample, at inter-point distances mm. Figure 2*d* also shows that no points fall within approximately 25 μm of each other and that, for *r*<0.1 mm, points are more regular than they would be if distributed uniformly randomly (resembling a hard-core process). This in part reflects the finite size of terminal villous branches, of typical diameter 50 μm (figure 1*b*).

### (b) Accuracy of the homogenization approximation for periodic sink distributions

We now explore how the spatial distribution of villous branches, represented by point sinks in equation (2.5), the strength of flow past these sinks (represented by *Pe*) and the sink strength (represented by *Da*) together determine patterns of uptake. The placenta acts as an exchange organ for numerous different solutes, each of which is characterized by its own values of *Pe* and *Da*. We wish to assess the accuracy of the leading-order homogenization approximation *C*^{(0)} across parameter space.

We first review concentration profiles over a periodic array (figure 3). For a given *Pe*, there is a critical *Da* (see equation (A6*b*)), equivalent to
3.1
such that the solute is fully absorbed within the domain (i.e. for *X*≤*X*_{0}<1) for *Da*>*Da*_{cr}(*Pe*). This threshold in (*Pe*,*Da*)-space asymptotes to the boundary between asymptotic regions U and U_{D} for *Pe*≪*ε* (when *Da*_{cr}≈2*ε*^{2}) and the boundary between regions A and U_{A} (and A^{S} and U^{S}_{A}) for *Pe*≫*ε* (when *Da*_{cr}≈*εPe*). Thus, equation (3.1) demarcates a region where uptake by sinks can be considered optimal: for *Da*>*Da*_{cr}, all the solute is absorbed upstream of the outlet, making some sinks redundant; for *Da*<*Da*_{cr}, substantial solute escapes past the sinks to the outlet. Of particular relevance for common nutrients such as glucose and oxygen in the placenta is the regime in which *Pe* is of order unity or larger [6], which we illustrate using point (5) in figure 3. In this case, the correction to *C*^{(0)} develops an internal boundary layer within each unit cell. For sufficiently large *Da* and *Pe*, this leads to a ‘staircase’ structure in the exact solution (figure 3*b,c*). Leading-order solute distributions, such as those computed by Chernyavsky *et al.* [6], fail to capture these significant local variations in the solute field.

We now survey (*Pe*,*Da*)-space to assess regimes of weak or strong convergence of residue (2.8) (i.e. or as respectively). For *Pe*≪1, we find that ∥*r*^{ε}∥_{L2}=*O*(*Da*) and ∥*r*^{ε}∥_{H1}=*O*(*Da*/*ε*) (using equation (A8)), implying strong convergence for *Da*≪*ε*. This reveals a region of parameter space (*Pe*≪1, *ε*≪*Da*≪1) where the correction to *C*^{(0)}, while bounded, has large gradient. For *Pe*≫1, we find that ∥*r*^{ε}∥_{L2}=*O*(*Da*/*Pe*) and (using equation (A13)), implying strong convergence for . We can, therefore, split (*Pe*,*Da*)-space into three regions (figure 3*a*): in the unshaded region below the dashed line there is strong convergence of the homogenization approximation; the weakly shaded region above the dashed line is characterized by weak convergence in which *C*^{(0)} fails to capture large gradients at the microscale; and the homogenization approximation fails in the dark-shaded region at large *Da*. In the regime of particular physiological significance (point (5)), convergence is weak.

### (c) Accuracy of the homogenization approximation for random sink distributions

The homogenized leading-order solution (A6a) is applicable not only to a periodic array but also for a statistically homogeneous random distribution of sinks of the same average density. We now test it against simulations using the three stochastic sink distributions described in §2*b*: a uniformly random distribution (*f*_{u}), a Matérn hard-core distribution (*f*_{h}) and a perturbed periodic distribution (*f*_{n}).

We consider point (5) in figure 3*a*, for which *C*^{(0)}=1−*X*. As advection dominates diffusion at the microscale, the exact solution for any realization of the sink distribution becomes an irregular staircase, exhibiting distinct plateaux between clusters of sinks. Even though the sink distributions are stationary processes that are not correlated over multiple sinks (§2*d*), the residues *r*^{ε} for *f*=*f*_{u} and *f*_{h} (figure 4*a,b,d,e*) show long-scale correlation: the variance is proportional to *X*(1−*X*) and the transverse covariance (proportional to *X*^{2} in 0<*X*<1/2 and (1−*X*)^{2} in 1/2<*X*<1) suggests that fluctuations are correlated over lengthscales comparable to the full domain, the correlation being greatest away from boundaries. *r*^{ε} for *f*=*f*_{n} with *σ*=10 also shows evidence of long-scale correlation (figure 4*c,f*); in this case, the variance and transverse covariance are approximately uniform across the domain (outside boundary layers at *X*=0 and *X*=1).

Figure 4*g,h* shows how the magnitude of the residue depends on *d* (for *f*_{h}) and *σ* (for *f*_{n}), for different values of *ε*. (Recall that, for *d*=0, *f*_{h} is equivalent to *f*_{u}, and that *f*_{n} resembles *f*_{u} for sufficiently large *σ*.) Collapse of the data for different *ε* indicates that, in both cases, the residue is *O*(*ε*^{1/2}) for sufficiently small *d* and sufficiently large *σ* (as expected from [18]), although the error falls in magnitude as the distributions become more regular (either by increasing *d* towards *d*_{cr} or by reducing *σ* close to zero). Indeed, for *σ*=0, the residue has exactly the scaling predicted by asymptotics for the periodic sinks, namely ∥*r*^{ε}∥_{L2}≈0.451*ε* (for *Pe*=10, *Da*=*εPe*). The magnitude of the residue approaches the value of the uniformly random distribution for .

When diffusion dominates at the microscale (*Pe*=*O*(*ε*), *Da*=*O*(*ε*^{2})), the point-wise variance again varies smoothly over the whole domain, as we show analytically in equation (2.10) (see also appendix B). However, its magnitude depends strongly on the degree of periodicity in the underlying structure, with fluctuations rising from *O*(*ε*^{3/2}) for almost periodic sink distributions to *O*(*ε*^{1/2}) for uniformly random sink distributions. Correspondingly, the range of validity of the homogenization approximation when *f*=*f*_{u} is significantly smaller than in the periodic case: we estimate that this requires *Da*≪*ε*^{3/2} for *Pe*≪*ε* and *Da*≪*ε*^{1/2}*Pe* for *Pe*≫*ε*.

## 4. Discussion

Heterogeneous and disordered biological media such as the human placenta require careful consideration of statistical variability when simulating transport processes, both in characterizing the underlying geometry and in understanding the impact of randomness on physical processes. Regarding the former, we have illustrated how sampling villous area fraction and estimation of the *K*-function (figure 2*c,d*) reveal important intrinsic lengthscales in the distribution of villous branches. For the sample considered, our data show no evidence against uniformly random distribution patterns of villous branches over sufficiently large distances and no evidence of underlying periodicity that would give rise to clear steps in the *K*-function [33]. The *K*-function instead resembles a hard-core distribution at shorter lengthscales (figure 2*d*), consistent with the requirement that branches cannot overlap. Further development of semi-automated image analysis, e.g. in watershedding segmentation algorithms [26], should allow bulk processing of histological data and reduce systematic errors. Future studies can be used to fit parameters of suitable spatial models to histological data and to assess how these features of tissue architecture may vary during development, in disease and between individuals.

In the absence of highly resolved measurements of maternal blood flow and solute distributions within the intervillous space, we chose here to examine nutrient transport using a simple theoretical model that incorporates spatial disorder in sink distributions. The model ignores many important features of placental haemodynamics, such as the finite size of villous branches, complex uptake kinetics for specific nutrients (glucose, amino acids, lipids, etc.), which may be mediated by active transport mechanisms in the syncytiotrophoblast on the outer surface of villous branches [34], the non-Newtonian rheology of flow in intervillous space (and the associated haematocrit distributions that would determine the oxygen-carrying capacity of maternal blood), high-dimensional transport processes such as Taylor dispersion [35], flow and tissue inhomogeneities (e.g. near spiral artery outlets [6]) and distortion of villous branches by flow. It is important that future transport models including each of these refinements are evaluated against data as they become available. Instead, we have used our simplified model to illustrate some generic features of homogenization approximations for disordered media, which we hope will provide a useful foundation for more sophisticated simulations of placental function.

While periodic or random sink distributions share similar leading-order approximations (through the slowly varying concentration field *C*^{(0)}(*X*)), the accuracy of this approximation depends on parameter values (*Pe*, *Da* and *ε*) and sink statistics. We classified asymptotic parameter regimes in (*Pe*,*Da*)-space (figure 3*a*), in which different balances between diffusion, advection and uptake are locally dominant at the macroscale. For the regime of greatest interest physiologically (point (5)), convergence of *C*^{(0)} to the exact solution is weak (applying in the *L*_{2} but not the *H*^{1} norm), because corrections have large spatial gradients on lengthscales below the inter-sink distance. Such fine-scale features, not captured by leading-order homogenization, are likely to be of importance in models that resolve the detailed arrangement of foetal vessels within villous branches.

The magnitude of the difference between the homogenization approximation and the exact solution depends on how one chooses to measure it. In a weak (*L*_{2}) norm, the residue with a periodic sink distribution is typically *O*(*ε*) (e.g. for *Pe*=*O*(1), *Da*=*O*(*ε*)), falling to *O*(*ε*^{2}) at sufficiently low *Pe* and *Da*. However, when sinks have a uniformly random distribution, the residue (in the appropriate norm) rises to *O*(*ε*^{1/2}) in both cases (equation (2.10) and figure 4*a,b*). The magnitude of the residue falls for distributions with a greater degree of periodicity (equation (2.10) and figure 4*g,h*) but grows with increasing sink strength. Significantly, even when sink distributions are correlated only over short distances, the residues appear to be correlated over distances comparable to the domain size when advection dominates at the macroscale (figure 4*d–f*). This is also the case when diffusion dominates at the microscale, as revealed by estimates of the transverse covariance (data not shown). While the distribution of villous branches in a placental unit (figure 2) is significantly more complex than the simple distributions employed in figure 2, one can estimate *ε* crudely to be between 0.001 and 0.01, suggesting errors in homogenization approximations owing to stochasticity of up to 10 per cent that fluctuate across distances comparable to an individual lobule.

Homogenization approaches provide a powerful tool in multi-scale modelling and are likely to figure prominently in future integrative models of other tissues with fine-grained periodic or random microstructure (such as the liver [36]). However, this study shows the merits of stepping beyond the leading-order approximation in order to resolve fine-scale structures at the microscale and, perhaps more importantly, to assess carefully the magnitude and nature of cumulative (and parameter-dependent) errors that arise from stochastic variation. These errors must be interpreted using the language of non-smooth functions and distributions. Such steps will be particularly important when building complex models that integrate numerous competing processes, in order to avoid errors arising at each level of approximation from accumulating and disrupting the overall predictive capacity of the model.

## Acknowledgements

The authors would like to thank Prof. Vincenzo Capasso for helpful discussions and Dr Ruta Deshpande and Ms Marie Smith for assistance with imaging. I.L.C. was supported by the Marie Curie Network MMBNOTT (project no. MEST-CT-2005-020723) and O.E.J. by the Leverhulme Trust.

## Appendix A. Homogenization for periodic sink distributions

We seek solutions of equation (2.5) in the form , introducing a slowly varying spatial variable *X*=*εx* that becomes independent of *x* as (at least in a sense of weak two-scale convergence [11]), allowing derivatives to be expressed as:
A 1
We further assume that is an *x*-periodic function with a period 1. Substituting equation (A 1) into equation (2.5) gives
A 2a
A 2b
A 2c
A 2d
for *n*=1,2,…,*N*, setting *X*_{0}=*εx*_{0} and assuming continuity of across sinks. Here, [*f*]_{x=n}≡*f*|_{x=n+}*f*|_{x=n−}.

We first address the case *Pe*=*O*(*ε*) and *Da*=*O*(*ε*^{2}), setting *Pe*=*εp*, *Da*=*ε*^{2}*q* and assuming *p*,*q*=*O*(1) as *ε*→0. We substitute the expansion
A 3
into equation (A 2) and collect terms in powers of *ε*, assuming the *C*^{(i)} are *O*(1) as . At *O*(1), we find that *C*^{(0)}=*C*^{(0)}(*X*). The homogeneous problem at the following order has the solution *C*^{(1)}=0. Competition between advection, diffusion and uptake emerges at *O*(*ε*^{2}), where
A 4a
A 4b
A 4c
It is convenient here to consider a representative unit cell occupying −1/2<*x*<1/2 with the sink at *x*=0. Averaging equation (A 4*a*) over a unit cell, and assuming periodicity of , we recover
A 5a
A 5b
A 5c
for some 0≤*X*_{0}≤1. The solution to equation (A 5) is
A 6a
where
A 6b
The internal free boundary at *X*=*X*_{0}<1 must be considered when the sink strength is large enough for all the solute to be absorbed upstream of *X*=1. The two-parameter solution (A 6) encompasses four distinct limits (illustrated in figure 3), denoted D, A, U_{D} and U_{A}, in which different effects dominate, respectively. On the boundaries between these regions of the parameter space, we recover one-parameter solutions that capture a balance between two effects. Specifically, diffusion between the boundary source and sink balances advection along *p*=*O*(1), *q*≪1, where equation (A 6) reduces to *C*^{(0)}≈1−(e^{pX}−1)/(e^{p}−1), independent of the distributed sinks; diffusion between the boundary source and sinks balances uptake by the distributed sinks along *q*=*O*(1), *p*≪1, where equation (A 6) reduces to *C*^{(0)}≈1/2*qX*^{2}−(1/2*q*+1)*X*+1 for *q*≤2, independent of advection; and advection balances uptake by the distributed sinks along *q*=*O*(*p*), *p*≫1, where *C*^{(0)}≈1−(*q*/*p*)*X*+e^{−p}((*q*/*p*)−1)(e^{pX}−1) (for *q*≤*p*). This limit includes a boundary layer of width 1/*p* near the boundary sink at *X*=1, across which advection balances diffusion locally. A further balance arises along *q*=*O*(*p*^{2}), *q*≫1, where *q*>*Q*(*p*) and diffusion and advection locally balance uptake (independent of the downstream boundary sink); here *C*^{(0)}≈1−(*q*/*p*)*X*+(*q*/*p*^{2}) e^{−pX0}(e^{pX}−1) and *X*_{0}≈2/*p*.

Having determined *C*^{(0)}, the second correction *C*^{(2)} satisfies
A 7
subject to periodicity over the unit cell and a calibration condition 〈*C*^{(2)}〉=constant, where . Thus
A 8
We impose *C*^{(2)}|_{x=0+}=*C*^{(2)}|_{x=0−}=0, so that *C*^{(2)} satisfies the boundary conditions (A 4*b,c*), implying that 〈*C*^{(2)}〉=*q*/12. This signals the development of pronounced humps in the concentration profile as *q* increases. We investigate these by inspecting the second organizing centre in the parameter space, at *Pe*=*O*(1), *Da*=*O*(*ε*).

We set *Pe*=*O*(1) and *Da*=*εq*_{1} with *q*_{1}=*O*(1) as . Expanding as before using equation (A 2) and collecting terms in powers of *ε*, we again find at leading order that *C*^{(0)}=*C*^{(0)}(*X*). At *O*(*ε*^{1}),
A 9a
A 9b
A 9c
Averaging equation (A 9*a*) over the unit cell, assuming periodicity of and incorporating the jump in at the sink, leads to , *C*^{(0)}|_{X=0}=1, which has the linear solution
A 10
where *q*_{1}/*Pe*≤1 if the concentration profile extends to the outlet (*X*=1), and *q*_{1}/*Pe*=1/*X*_{0} when the solute concentration drops to zero at *X*=*X*_{0}<1. In order to satisfy the downstream boundary condition on *C*^{(0)} at *X*=1 or *X*=*X*_{0}, it is necessary to include a diffusive boundary which is not preserved in this scaling, but which is as described in regions U_{A} and A above.

To determine *C*^{(1)}, we note that the source terms and *q*_{1} in the linear system (A 9*a,b*) are independent of *x*, so that superposition may be used to find a solution in the form
A 11
where *a*(*x*) and *b*(*x*) are continuous with unit period and satisfy the cell problems *a*_{xx}−*Pe* *a*_{x}=0, , *b*_{xx}−*Pe* *b*_{x}=−1, . After some algebra, it follows that
A 12a
for 0<∓*x*≤1/2. As before, the global boundary conditions (A 9*c*) are used to derive the local condition *C*^{(1)}(0)=0, implying
A 12b
therefore 〈*C*^{(1)}〉≈(*q*_{1}/2*Pe*)(1−2/*Pe*) at large *Pe*≫1 and 〈*C*^{(1)}〉≈*q*_{1}/12 at small *Pe*≪1, so that *C*^{(1)} in equation (A 12a) tends to *εC*^{(2)} in equation (A 8) with error *O*(*q*_{1}*Pe*). For large *Pe*, the correction is approximately
A 13
the discontinuity across the sink being smoothed over a distance (in *x*) of 1/*Pe*. This gives rise to staircase-like solutions in regions U and A^{s} in figure 3.

We also observe from equation (A 12a) that the asymptotic approximation for *Pe*≫1 breaks down when *q*_{1}/*Pe*=*O*(*ε*^{−1}), i.e. *Da*∼*Pe*, making *εC*^{(1)} of the same order as *C*^{(0)} in equation (A 3). Overall, the homogenization approximation applies for , as shown in figure 3*a*.

## Appendix B. Homogenization for random sink distributions

When sinks are distributed non-periodically, we can derive the homogenized approximation of equation (2.5) as follows. We focus on the case *Pe*=*O*(*ε*), *Da*=*O*(*ε*^{2}), again writing *Pe*=*εp* and *Da*=*ε*^{2}*q*. We initially use equation (A 1) to rewrite equations (2.5) and (2.6) as
B 1
(For brevity, we assume here that does not fall to zero upstream of *X*=1.) Expanding using equation (A 3), we allow *C*^{(1)} and *C*^{(2)} to have fluctuations, assuming that these are not large enough to disrupt the proposed expansion. At leading order, , *C*^{(0)}|_{X=0}=1 and *C*^{(0)}|_{X=1}=0. Thus, for some and . The first term must be suppressed to avoid secular growth, so that *C*^{(0)}=*C*^{(0)}(*X*). Likewise, at the following order, we find that *C*^{(1)}=*C*^{(1)}(*X*). Collecting the terms in equation (B 1) at *O*(*ε*^{2}), we obtain
B 2
with *f* given by equation (2.6). This is to be solved subject to *C*^{(1)}=*C*^{(2)}=0 at *x*=0 and *x*=*ε*^{−1}. Thus in *ξ*_{i}<*x*<*ξ*_{i+1}, for *i*=0,1,2,…,*N*, treating *x* and *X* as independent,
B 3
for some *α*_{i}, *β*_{i}, taking *ξ*_{0}=0 and *ξ*_{N+1}=*ε*^{−1}. We define Δ_{i}≡*ξ*_{i}−*ξ*_{i−1} (for *i*=1,2,…,*N*+1) and so that and
B 4
Integrating equation (B 2) across *x*=*ξ*_{i} gives, for *i*=1,2,…,*N*,
B 5a
and
B 5b
We take *β*_{0}=0 to satisfy *C*^{(2)}=0 at *x*=0. Here, we make the assumption, validated *a posteriori*, that *F* is uniform. From equation (B5), for *i*=1,2,…,*N*,
B 6a
and
B 6b
Substituting equation (B 6) into equation (B 3), expanding and using equation (B 4) gives, after some algebra,
B 7
Imposing *C*^{(2)}=0 at *x*=*ε*^{−1}≡*N*+1 gives , and so
B 8
This expression relates solute fluctuations directly to sink distributions.

When *f*=*f*_{n}, then , and . To evaluate the distribution of *εxR*_{N}−*R*_{i}, we must sum independent distributions. Setting *x*=*X*/*ε*,
writing *i*=(*X*/*ε*)+(*i*−*x*), retaining (*i*−*x*) as an *O*(1) quantity. We therefore find from equation (B 8) that (after some algebra), for *ξ*_{i}≤*x*<*ξ*_{i+1},
B 9

To ensure that the original expansion is asymptotic, we must therefore take *F*=1 at *O*(*ε*^{−2}), yielding equation (A 5*a*) for *C*^{(0)}. Simulations indicate that the contribution at *O*(*ε*^{−1}) (and hence *C*^{(1)}) vanishes. Similarly,
B 10
we find that equation (B 10) is in agreement with simulations (not shown). Cov_{T}(*C*^{(2)}) can be determined in a similar manner. Thus, while *C*^{(2)} has *O*(1) mean, *r*^{ε} is dominated by fluctuations of relative magnitude *O*(*ε*^{3/2}). This approximation holds as long as sinks do not exchange places, which can be expected once *σ* becomes sufficiently large. Equation (B 10) suggests that the fluctuations in the case of stronger mixing of sink locations will be larger than *O*(*ε*^{3/2}).

When *f*=*f*_{u}, we referred to Matsunawa [29], who determined the distribution of linear combinations of order statistics drawn from (i.e. combinations of , where ) as a mixture of scaled beta distributions. Writing for some *a*_{i}, and defining , *L*_{N} has mean and variance
B 11
respectively [29]. Writing and setting , we have *a*_{j}=*X*−1 for 1≤*j*≤*i* and *a*_{j}=*X* for *i*+1≤*j*≤*N*, and hence
B 12
Evaluating appropriate sums of *b*_{j} and again writing *x*=*X*/*ε*, *i*=(*X*/*ε*)+(*i*−*x*) (taking |*i*−*x*|=*O*(1)), we find to leading order in *ε* using equation (B 11) that *XR*_{N}−*R*_{i} has mean and variance
B 13
respectively. We thereby recover equation (B 9) (after some algebra), requiring that *F*=1+*O*(*ε*), so once more *C*^{(0)} satisfies equation (A 5*a*). Furthermore, assuming Var(*C*^{(1)})=0, we obtain to leading order in *ε*, consistent with simulations. Thus, fluctuations about *C*^{(0)} are *O*(*ε*^{1/2}). Furthermore, simulations show that , suggesting a contribution from *C*^{(1)} which presumably must be determined by a closure condition at higher order.

## Footnotes

One contribution of 11 to a Theme Issue ‘Towards the virtual physiological human: mathematical and computational case studies’.

- This journal is © 2011 The Royal Society