Optimal vein density in artificial and real leaves
 *Department of Organismic and Evolutionary Biology,
 ^{§}School of Engineering and Applied Sciences, and
 ^{¶}Arnold Arboretum, Harvard University, Cambridge, MA 02138; and
 ^{†}Laboratoire de Physique de la Matière Condensée, Centre National de la Recherche Scientifique–Unité Mixte de Recherche 6622, Université de NiceSophiaAntipolis, Parc Valrose, 06108 Nice Cedex 2, France
See allHide authors and affiliations

Edited by Karl J. Niklas, Cornell University, Ithaca, NY, and accepted by the Editorial Board March 31, 2008 (received for review September 27, 2007)
Abstract
The long evolution of vascular plants has resulted in a tremendous variety of natural networks responsible for the evaporatively driven transport of water. Nevertheless, little is known about the physical principles that constrain vascular architecture. Inspired by plant leaves, we used microfluidic devices consisting of simple parallel channel networks in a polymeric material layer, permeable to water, to study the mechanisms of and the limits to evaporationdriven flow. We show that the flow rate through our biomimetic leaves increases linearly with channel density (1/d) until the distance between channels (d) is comparable with the thickness of the polymer layer (δ), above which the flow rate saturates. A comparison with the plant vascular networks shows that the same optimization criterion can be used to describe the placement of veins in leaves. These scaling relations for evaporatively driven flow through simple networks reveal basic design principles for the engineering of evaporation–permeationdriven devices, and highlight the role of physical constraints on the biological design of leaves.
Leaves are the main pumps for the upward movement of fluids in plants, which takes place in a network of conduits collectively called xylem (1). On warm, sunny days, the leaves of a large tree can transport hundreds of liters of water from the soil, corresponding to several milliliters per day per leaf. Although photosynthesis, rather than water loss, is the desideratum of living leaves, the shared diffusional pathway for CO_{2} and water vapor and the regulation of this pathway by water status in the leaf implies that these two fluxes are closely coupled. Indeed, maximum rates of photosynthesis scale linearly with the transport capacity of the xylem (2–4). A key attribute of vascular design within leaves is to ensure that no portion of the leaf is undersupplied and thus unable to support the water losses associated with CO_{2} uptake (5–8). However, studies of leaf venation have largely treated leaves as twodimensional structures (8). Here, we account for the fact that leaves have a finite thickness and address the placement of veins within leaves in terms of both their spacing and their distance from the evaporative surface.
Evaporationdriven flow in leaves is controlled by parameters such as wind, temperature, and stomatal aperture, all of which have been studied extensively (1, 9). However, the internal constraints imposed on evaporative transport by the architecture of the leaf hydraulic system have received relatively little attention (5, 6). Water utilizes two successive pathways as it flows from the veins of a leaf to the outside air, passing from the liquid to the vapor state. The first pathway is associated with capillary suction through the porous materiallike network of cells and cell walls. One can describe the flux J using Darcy's law:
In addition, a key element in the water balance of leaves is the presence of turgorcontrolled pores (stomata) that control the rate at which water vapor is lost via transpiration. Any attempt to understand the constraints on water delivery imposed by vein placement must account for active stomatal control. The dependence of stomatal aperture on Ψ means that Darcy's law describing liquid flow from veins to the sites of evaporation and Fick's law governing diffusion of water vapor through stomatal pores are coupled. The outcome of this coupling is that both the rate of evaporation controlled by stomata and the evaporation rate driven by diffusion are inversely related to the distance over which the gradients operate. This analogy allows us to mimic water transport in leaves using a simple physical model in which water moves from supply channels (veins) to the leaf surface solely by diffusion.
To examine the design constraints on evaporatively driven hydraulic systems (Fig. 1 and Experimental Methods) we use microfluidic systems that use permeationdriven flows (10, 11). The ability to manufacture and manipulate complex microfluidic geometries (12) using poly(dimethylsiloxane) (PDMS) allowed us to explore analogies to leaves in a quantitative, measurable way using our knowledge of water permeation through PDMS in microchannels. As shown recently (11), the transport mechanism is dominated by molecular diffusion of vapor driven by concentration gradients. Water molecules move in small clusters through the hydrophobic PDMS network from one pore to another over a scale of a few angstroms, just as water vapor diffuses in air, but with a much smaller diffusion coefficient because of the presence of hydrophobic barriers. This simple picture is valid in the immediate vicinity of a single channel of water in the PDMS leaf analog.
For the case of multiple sources or channels shown in Fig. 1, if the thickness of the PDMS layer is small compared with the spacing between the channels, i.e., δ ≤ d, each channel provides an evaporative surface flux that is independent of its neighbor. However, when the layer thickness is large compared with the distance between channels, δ ≥ d, the channels collectively act as a single one. With these simple ideas in mind, our goal was to determine, using both experiment and theory, the scaling relations that constrain evaporatively driven transport through simple microfluidic systems and to compare the design optima for artificial networks with the internal architecture of real leaves, by determining the role of vein (or channel) geometry, spacing, leaf thickness, and external humidity in characterizing flow.
Results
Effects of Environment and Geometry.
By varying the external humidity imposed at the PDMS surface, we find that our devices behaved as free evaporative surfaces, with flow rates proportional to the relative dryness (= 1 − RH) imposed at the surface (Fig. 2). However, the rate of water movement was sensitive to channel density; at low relative humidities, a greater density of channels resulted in higher flow rates. The minimum channel density that can fully satisfy the evaporative demand is the key design parameter for natural and artificial pervaporative pumps. An oversupply of channels represents a significant inefficiency, both in terms of construction costs and the consumption of space that might otherwise be devoted to other functions. At low channel density (δ ≪ d, where d is the distance between channels), each channel functions independently such that the total flow rate is proportional to channel number. The concentration of water within the PDMS layer decreases toward the air, with the exact concentration profile dependent on channel shape and size. Except in the region very close to the channel walls, the concentration decreases logarithmically as a function of distance from the channel, as expected from the solution of the Laplace equation (see Model: Analytical Approximations and Model: Numerical Simulation). The total flow rate (J _{T}) is the sum of the flux from each channel, which can be written as where D is the diffusivity of water vapor in PDMS, c _{0} and c _{1} are, respectively, the concentrations at the channel and evaporative surface, L is the channel length, h is the channel height, d is the distance between channels, W is the width of the microfluidic device, a is the equivalent radius of a semicircular channel (see Model: Analytical Approximations and Model: Numerical Simulation), and N ∼ W/d is the number of channels.
Conversely, for high channel density (δ ≫ d), the vapor concentration decreases linearly toward the PDMS–air interface such that the flow rate is independent of channel density and the flux is then In Fig. 3, we show that flow rate through artificial leaves of varying thickness is in good agreement with theory; for a given thickness, it increases with channel densities and eventually saturates. Furthermore, the maximum flow rate J _{m} ∝ δ^{−0.92} (Fig. 4), close to the simple prediction for diffusively limited transport in a 1D medium.
Optimization.
Our preceding results suggest that for a given PDMS layer thickness, there is a characteristic channel spacing (d _{c}) below which there is no additional gain in the flow rate because the system behaves as if it were effectively a quasi onedimensional medium. Equating the flux when d ≫ δ (Eq. 1 ) to the flux when d ≪ δ (Eq. 2 ), with n = W/d, we find
In Fig. 3, one can check that the characteristic channel spacing above which the flow saturates for each thickness is in good agreement with Eq. 3 . In the general case, corresponding to different channel shapes and dimensions, this exact result is replaced by a simple, intuitive scaling law, as expected on dimensional grounds for the transition from twodimensional behavior to onedimensional behavior that occurs when the channel spacing becomes comparable with the PDMS layer height and the channels act collectively rather than as individuals. Indeed, not only is there no gain in flow by adding more channels for pervaporative pumps of a given thickness, but when the channels are too close to each other, the surface area over which PDMS adheres with glass decreases so that the maximum pressure sustainable by the device is lowered. This suggests that for this particular artificial leaf, a compromise between adhesive strength and flow rate determines the optimal channel spacing. More generally, if there is a cost associated with making veins (e.g., ensuring that adhesion is large enough, or in a biological context, a metabolic cost), the scaling relation Eq. 4 provides a simple optimization principle for the design of pervaporative devices.
Real Leaves.
The diversity of leaf venation found in nature provides an ideal opportunity to test the generality of the scaling relations demonstrated for artificial leaves. Most leaves exhibit a hierarchical and reticular venation, with the larger veins serving as major supply lines for the smaller veins that serve as the principal sites where water moves from the xylem to the mesophyll. Because we are interested in the ability of the venation to supply the evaporative demand, rather than distribute water across the leaf as a whole, we focused on the placement of only the highestorder veins (7). Although there is a large literature on the shape (13) and sizes of venation patterns, there are few studies (7) ^{‖} that simultaneously report leaf thickness, and only a recent paper (14) provides information on the physically meaningful parameter of the distance from the xylem to the evaporative surface. We therefore measured both the average distance between the highestorder veins and the leaf evaporative surface (δ) and the average distance between adjacent veins (d), for 32 species characteristic of open, unshaded environments (see Experimental Methods). The fact that a leaf is reticulate or parallel has no bearing on either the determination of d or on the interpretation of these data using our model.
The venation of this diverse assemblage of leaves shown in Fig. 5 is well described by the simple law (Eq. 4 ), suggesting that leaves do not produce veins in excess of what are needed to allow maximum rates of water delivery to the leaf surface. Factors contributing to this optimum include both the costs of constructing the veins and the utilization of space that could otherwise be filled with photosynthetic cells. A simple fit gives d = 1.08 δ, R = 0.989, quantitatively consistent with Eq. 4 . However, one has to keep in mind that our model only provides a scaling relationship, and a calculation of the exact numerical prefactor of δ remains beyond the ability of any simple model unless we understand the multiple parameters that affect the design of real leaves. Nevertheless, our scaling law provides a theoretical framework for the recently reported inverse correlation between photosynthesis and distance for water movement between veins and stomata (14).
Discussion
In this study, we have explored only the simplest topology of channels in our artificial leaves, ignoring variations in the morphology and architecture that characterize the complex venation systems of real leaves. The general scaling principles that our biomimetic leaves follow in response to the constraints from the evaporative surface reflects the design requirements for optimal channel placement. Taken together with the strong correlation between density of veins and their distance to the leaf surface in real leaves, this suggests that design features oriented toward optimal utilization of veins in providing water for evaporation are universal and thus independent of the morphological diversity among leaves. Because of the analogy of transpiration rate variation with distance due to active stomatal control and of diffusion in PDMS leaves, our study also provides further evidence of the homeostatic regulation of leaf water status via stomatal control. For natural systems, the existence of architectural optima describing the placement of veins in leaves provides an important foundation for understanding leaf function (5) as well as the physical/physiological limits imposed by the distribution system itself on development (8). For biomimetic systems, our study points toward principles for their design based on an understanding of the underlying physical limits.
Methods and Models
Experimental Methods.
Microchannels were fabricated by using PDMS, as described in ref. 15, with 10:1 PDMS:crosslinker mixture (Sylgard 184), spincoated onto a negative photoresist SU8 mold, degassed for 15 min at a 4inch Hg vacuum pressure, and cured at 65°C for 4 h. The channels (width ω = 15 μm, height h = 30 μm) were made in a thin layer bonded to a glass slide by plasma treatment, its thickness adjusted by the spincoating speed. A thick layer was bonded by plasma treatment on top of the thin layer to make the inlet for the liquid and to act as a frame, avoiding wrinkles in the PDMS film. This layer presents a rectangular window the size of the evaporative surface. The channels were slowly filled with water, the air escaping through the PDMS, because the channels had sealed ends. The devices were stored in a wet environment to keep the channels filled with water because contact with air results in the PDMS surfaces regaining their hydrophobicity (16), which makes refilling them more difficult. The channels (length L = 35 mm), were connected over a width W = 30 mm to a large inlet channel (100 μm wide) that feeds them. The spacing between two channels d ∈ [30 μm to 1,500 μm]; and the PDMS thickness above the channel δ ∈ [45 μm to 760 μm] (Fig. 1). The channel density is then n = 1/d, and the number of channels over 30 mm is n = W/d.
The devices were connected below them (30 cm) to a beaker resting on the plate of a precision balance (Sartorius BP 211 D) to ensure that water was pumped upward. Connected to a computer, this gives a direct access to the mass flow rate J, provided that it is sufficiently high (>0.1 μg·s^{−1}). Typical value were between 0.1 and 10 μg·s^{−1} (or between 0.1 and 10 nl·s^{−1}), the latter corresponding to an average speed in the main channel of several centimeters per second. We controlled the evaporative driving force imposed at the surface by blowing air with precisely regulated humidity (20–100% RH) from a dew point generator (LICOR 1600). The air stream produced by the dew point generator was directed across the top PDMS surface at a rate of 1 liter/min. A Plexiglas plate was placed 1 mm above the PDMS surface to reduce variation in the evaporation rate due to disruptions of the boundary layer.
Freehand crosssections of fresh leaf material were examined under a microscope, and the average distance between veins and from veins to the abaxial leaf surface was measured by using AxioVision 4.5 (Zeiss), calibrated for each magnification (×40–400). Species examined were as follows: Acer rubrum L., Anacardium occidentale L., Apeiba membranacea Spruce ex Benth., Araucaria araucana (Molina) K. Koch, Avicennia germinans L., Begonia sp. L., Calophyllum brasiliense Cambess., Castilla elastica Sessé subsp. costaricana (Liebm.) C. C. Berg, Citrus sinensis (L.) Osbeck, Clusia sp. L., Crassula ovata (Miller) Druce, Dipteryx panamensis (Pittier) Record & Mell, Eichhornia crassipes (Mart.) Solms, Euphorbia milii Desmoul., Guzmania monostachia (L.) Rusby ex Mez., Hernandia didymantha Donn. Sm., Hoya carnosa (L. f.) R. Br., Kalanchoe sp. Adans., Liriodendron tulipifera L., Nothofagus alpina (Poepp. & Endl.) Oerst., Nothofagus dombeyi Mirb. (Oerst.), Passiflora sp. L., Phalaenopsis sp. Blume, Pouteria calistophylla (Standl.) Baehni, Quercus rubra L., Rhizophora mangle L., Simarouba amara Aubl., Simira maxonii (Standl.) Steyerm., Sterculia recordiana Standl., Swartzia simplex (Sw.) Spreng., Warszewiczia coccinea (Vahl) Klotzsch, and Zamia pumila L.
Model: Analytical Approximations.
Assuming that the length of the channels is much longer than both the distance between them and the height of the PDMS film, we can neglect lengthwise variations in the water concentration profile. Of course, this is no longer true near the entrance regions and near the end of the channels, but these edge effects are small, and we neglect them here. Then, we may consider a 2D problem for the evaporation of water from the channels that satisfies the Laplace equation, i.e., with c(x, z, t) being the concentration, z being the coordinate along the direction of evaporation, and x being the coordinate perpendicular to it, along the channel width. The associated boundary conditions are (i) at the air interface c = c _{1}, which we will implement eventually using the technique of images, and (ii) there is no flow on the boundaries x = d/2 and x = −d/2.
For the general case of a periodic variation of the concentration in the PDMS, we can decompose the solution of the above equation in a Fourier series of the form Σ_{n=1} ^{∞} a_{n} cos(nπx/d)exp(−nπz/d) + bz, where the first term describes variations due to the individual channels, and the second term is the simple linear profile due to a constant gradient of concentration. We see that for thicknesses larger than the distance between channel, δ > d, the periodic terms vanishes and only the 1D term remains so that the concentration decreases linearly toward the air interface, showing why the transition between a 2Dlike and 1Dlike diffusion occurs for δ ∼ d. When d ≫ δ, however, the channels are effectively decoupled, so that it is sufficient to consider a unit cell around one channel (−d/2 < x < d/2). Replacing the channel of rectangular crosssection with an effective channel with a semicylindrical one (for simplicity), of radius a = (2Wh/π)^{1/2}, yields the result that the concentration decreases logarithmically with distance from the channel. In reality, channels are rectangular, but their effective radius a can be determined by using numerical simulation (described later) and leads to very similar results for the range of values for δ/d we have used. Then, one finds that the flux at the film–air interface for a single channel is so that the total flux for N = W/d channels is For d ≫ δ, this leads to a linear relation for the flux as a function of the number of channels: Extrapolating Eq. 6 into the limit d ≪ δ (or N → ∞), we find that This expression is close to the one obtained from Eq. 2 (with a numerical coefficient ln [(δ + h)/a][(δ+h)/δ]. This is not surprising as the channels are no longer independent in this limit; nevertheless it still gives the appropriate form of a linear increase followed by a plateau as a function of channel density.
Model: Numerical Simulation.
We used a finite element method implemented in Matlab (pdetool) to solve the steady state diffusion equation for the water vapor concentration in two dimensions: Δc = 0, using the particular geometry of our system. We assume that the concentration and fluxes are periodic in the x direction. We further assume that there are no edge effects on the three water–PDMS interfaces of the channels so that c(x, z) = c _{0}. On the PDMS–glass interface of the channels, the vertical flux is assumed to be zero, J = −D(dc/dz) = 0, while on the two sides of the unit cell of the artificial leaf, symmetry dictates that the transversal flux J = −D(dc/dx) = 0. Finally, on the PDMS–air interface, we assume c = c _{1} = % RH × c _{0}. The parameter values used are (11) D = 8.5 10^{−10} m^{2}/s and c _{0} = 40 mol/m^{3}. The solid curves in Fig. 3 represent the results of our simulation for the flux J at the surface. We find good agreement with our measurements, both for the saturation curve and the absolute flow rate (it is slightly lower in our measurements at low thicknesses). For comparison with the analytical approximation, we took the radius of the cylinder to be a = 18 μm to mimic the rectangular channel of dimensions ω × h = 30 μm × 15 μm (to keep the area constant, a = (2ωh/π)^{1/2}, giving a = 17 μm), leading to results that match the numerical simulation well over the entire range of thicknesses.
Acknowledgments
We thank Jacques Dumais for his generous support and advice, Abraham Stroock for his comments on the manuscript, the Materials Research Science and Engineering Center at Harvard University and National Science Foundation (IOB0517071) for funding, and the Center for Nanoscale Systems for technical support.
Footnotes
 ^{‡}To whom correspondence should be addressed. Email: xavier.noblin{at}unice.fr

Author contributions: X.N., L.M., D.A.W., N.M.H., and M.A.Z. designed research; X.N., L.M., I.A.C., N.M.H., and M.A.Z. performed research; X.N., L.M., N.M.H., and M.A.Z. analyzed data; and X.N., L.M., N.M.H., and M.A.Z. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission. K.J.N. is a guest editor invited by the Editorial Board.

↵ ‖ In this reference, qualitative correlation is made between vein density and tissue proportions, which is not necessarily related to δ.
 © 2008 by The National Academy of Sciences of the USA
References

↵
 Taiz L ,
 Zeiger E
 ↵
 ↵
 ↵

↵
 RothNebelsick A ,
 Uhl D ,
 Mosbrugger V ,
 Kerp H
 ↵
 ↵
 ↵
 ↵
 ↵

↵
 Randall GC ,
 Doyle PS
 ↵

↵
 Bohn S ,
 Andreotti B ,
 Douady S ,
 Munzinger J ,
 Couder Y

↵
 Brodribb TJ ,
 Feild TS ,
 Jordan GJ
 ↵
 ↵