M Leocmach

C P Royall

Hajime Tanaka

Novel zone formation due to interplay between sedimentation and phase ordering

Keywords: 82, 70, Dd -Colloids 47, 57, ef -Sedimentation and migration 64, 60, Cn -Order-disorder transformations

Usually, we expect large particles to sediment faster than small of the same material. Contrary to this intuition, we report a dynamical competition between sedimentation and phase ordering which leads to smaller particles settling faster than larger ones. We access this phenomenon using suspensions of polymers and two colloidal species which we image with confocal microscopy. Polymers mediate attractions between colloids, leading to phase separation and crystallisation. We find that the dynamical interplay between sedimentation, phase separation and crystal nucleation underlie this phenomenon. Furthermore, under certain conditions we find a kinetic pathway leading to an apparent coexistence between a one component crystal and binary fluid of equal buoyancy. These findings may be relevant to the basic understanding of sedimentation-induced zone formation in nature and industrial applications.

Introduction. -Under a gravitational field, solid particles of hard (metals, semiconductors and oxides) or soft matter (colloids, macromolecules and biological materials) are inhomogeneously distributed along the vertical direction: This process is called sedimentation. Key early advances include the discovery of sedimentation-diffusion equilibrium by Jean Perrin one century ago [START_REF] Perrin | Les Atomes[END_REF], where the number density ρ as a function of height z is given by [START_REF] Russel | Colloidal Dispersions[END_REF] P (z)/k B T = 1 lg h z ρ(z)dz where P (z) is the (osmotic) pressure, k B is Boltzmann's constant and T is temperature. Sedimentation may be characterised by the gravitational length l g = k B T /mg, where m is the buoyant mass of the colloidal particle and g is the gravitational acceleration. l g may be recast as a Peclet number (P e = τ D /τ S ), the ratio between the time τ D it takes a particle to diffuse its own radius and the time τ S to sediment its own radius. The Peclet number characterises the colloidal (P e 1) and granular limit (P e ≫ 1). An isolated spherical particle of diameter σ settles under gravity, at an average (Stokes') velocity u 0 = mg 3πησ . Here m = δρπσ 3 /6 with δρ the mass density difference between colloids and solvent and η the viscosity of the solvent. This leads to u 0 ∼ σ 2 : Larger particles sediment faster.

From this one-particle picture, we expect that a system with two or more species with different Stokes' velocity may separate during the process of sedimentation. However, reality is not necessarily that simple. In a suspension of a finite concentration, the kinetics of sedimentation are strongly influenced by hydrodynamic interactions between particles [START_REF] Russel | Colloidal Dispersions[END_REF][START_REF] Segre | [END_REF], and direct attractive interactions between the sedimenting particles can have a drastic effect and even qualitatively alter the simple picture above [4]. For example, in a colloid-polymer mixture (C+P), complex couplings between sedimentation and the polymer concentration come into play [5]. The polymer chemical potential drives the depletion attraction between the colloids, yet the polymer free volume is itself coupled to state point of the system [6]. Reflecting this feature, predictions from free volume theory suggest that highly nonintuitive behaviour occurs even in equilibrium in the form of a colloidal liquid floating above a colloidal gas [7]. Furthermore, it was shown that a system of C+P undergoes flocculation and then each floc settles much faster than individual particles, and can couple to phase ordering such as crystallisation during sedimentation [8]. In the sticky sphere limit (where the size of the polymer is very small compared to the size of the colloid) flocculation can lead to dense colloidal crystals or kinetically arrested gels [9].

Here we consider a system of binary colloid+polymer mixtures (C+C+P) with a very small overall concentrap-1 tion of one of the colloidal species. While the sedimentation behaviour of colloid-polymer mixtures has received some attention, its effect on binary hard-spheres is little understood. In bulk equilibrium binary hard spheres have been studied via experiment, [10][11][12], computer simulation [13] and theory [14][15][16][17], and exhibit a plethora of plastic crystals [18], binary crystal superlattices [13,17], and Laves phases [19]. Here the interplay between sedimentation and the equilibrium phase diagram (see fig. 1(a) [5])and phase ordering leads to fascinating unexpected kinetic pathways. We emphasise that C+C+P is not merely a complication of C+P, but a minimal system exhibiting both phase condensation, driven by depletion, and frustration against crystallisation due to polydispersity. The competition between these two effects is crucial to understand the formation of sediment layers in nature. Moreover, we will show that it is in fact possible to gain much insight into the physics of C+C+P mixtures where small colloids dominate the large by considering the unperturbed C+P case for the initial kinetics, and the binary hard sphere case to describe the dense sediment.

Experimental. -We used poly(methyl methacrylate) (PMMA) colloids, sterically stabilised with polyhydroxyl steric acid. The large (L) and small colloids (S) have diameters of σ L = 1.20 µm and σ S = 0.680 µm (ξ SL ≡ σ S /σ L = 0.57) and polydispersity 5% and are labeled with 4-chloro-7-nitrobenzo-2-oxa-1,3-diazol (NBD) and rhodamine isothiocyanate, respectively. The solvent was cis-decalin whose relative dielectric constant is 2.4. Results on the electrostatic charging in solvents of comparable dielectric constants [20] suggest that we expect only a handful (< 10) charges per particle in this system, implying electrostatic interactions between colloids which are inferior to k B T .

The polymer used was polystyrene (PS) with a molecular weight of 3.1 × 10 7 (M w /M n = 1.3). We estimate the polymer radius of gyration in an ideal solvent as R G = 160 nm: the theta temperature of a PS/cis-decalin mixture is 12˚C; following [21] we assume that our system is reasonably well-described assuming polymer ideality. Some deviations are observed, particularly at higher colloid volume fraction [22,23], however, we argue that a scaled particle approach is appropriate for the level of this work. This assumption was discussed in more detail in [23]. The polymer-large colloid size ratio is estimated as ξ P L ≡ σ P /σ L = 0.27 implying a polymer-small colloid size ratio of ξ P S = 0.47. According to [6], a S+P system would undergo gas-liquid phase separation, but not a L+P system, however experiments showed gas-liquid phase separation for a size ratio around 0.25 [24]. Neither can be considered as the sticky sphere limit.

Gravitational lengths for large and small colloids are l L g = 1.75 µm and l S g = 9.64 µm respectively. Peclet numbers are 0.684 and 0.0705 respectively. We set the respective quantities of large and small particles in order to have φ 0 S = 50φ 0 L . Here the superscript 0 denotes the volume fraction at time=0, i.e., the mean colloid volume fraction. Similarly, we set the polymer concentration to keep φ 0 P = 4R 3 G ρ P /3 = 1.8φ 0 S , where ρ P is the polymer number density. The solvent volume fraction is our experimental variable. Here we specify this solvent volume fraction by the polymer concentration φ 0 P . We used a Leica SP5 confocal microscope, using 488 nm and 532 nm laser excitation for the NBD and rhodamine labelled colloids respectively. Local volume fractions of large colloids were estimated by particle tracking. For the small colloids we used intensity measurements which were calibrated calibrated against homogeneous samples of known concentrations. In both pure fluid and pure crystal systems of small particles, we found linear dependencies of measured intensity against colloid volume fraction respectively. The volume occupied by the crystals were given by image segmentation. Rectangular sample cells (height: 40 mm, width: 5.0 mm and depth 2.0 mm) were used. Sedimentation is a slow process and we use the term "final state" to express an apparent steady state after the time needed for a single particle to sediment from the middle to the bottom of the cell (a few weeks). There is practically no observable change after two weeks. This final state is not the lowest free energy state, the time needed to reach the free energy global minimum runs to years, due to the slow diffusion of the particles in the dense sediment.

Results. -Since the vast majority of the colloids are small, it is instructive to consider the S+P system [fig. 1(a)], whose kinetic (and equilibrium) behaviour has been studied [5]. We identify three regimes of initial composition, I leading to colloidal fluid or phase separation into a colloid-rich and colloid-poor fluid phases, II leading to triple coexistence [5] and III leading to arrested phase separation (gel). We distinguish regimes I and II by whether the system undergoes homogeneous nucleation or not. In fact, regime I has the potential for colloidal gas-liquid type phase separation to dominate the early kinetics, which is similar to regime II. However, as we discuss below, it is in the longer-time kinetics that a clear difference between regimes I and II is seen.

In regime I, the final state resembles a binary hard sphere fluid (BHS): fluid at the top with a small φ L and crystals of S particles (XS) at the bottom with very few grain boundaries and some large particles included individually as defects in the XS lattice. These defects were more concentrated at the bottom of the crystals, as expected. Due to the small φ 0 L , no XL is observed. In other words, without phase separation of the majority (small) species, the behaviour is consistent with the naive expectation above [fig. 1(b,c,d Regime II. -More intriguing is the final state of C+C+P between the two previous cases, regime II, to which we devote the remainder of the discussion. A simple equilibrium expectation would be a final state very similar to regime I, but with a gas-liquid interface. However, we observed the following multi-zone structure from bottom to top (see fig. 2). At the bottom of the sediment, we find a binary amorphous solid (zone 1) [fig. 2(f)], followed by a coexistence between XS and a binary fluid (zone 2) [fig. 2(e)]. Crystals of small particles are limited in extent (20∼30 µm) and the binary fluid between them percolates from the upper to the lower bound of this zone. Neither melting nor growth of the crystals were detected, suggesting that crystal and fluid are really coexisting 'thermodynamically'. This is further supported by Supplementary movie 1 which reveals diffusion in the binary fluid, showing that the system has some degree of thermal motion and is not completely jammed. (zone 3) Crystals of small particles almost fill space but with a lot of grain boundaries [fig. 2(d)]. The size of each single crystal is 20-100 µm. φ L is very small. We refer to this zone as the "ice pack" because it is a crystalline state situated over a phase with a lot of fluid. (zone 4) The liquid phase is dominated by small colloids. At high altitude, we recover the gas phase of small colloids (zone 5) [fig. 2(c)].

The composition profile of each phase as a function of height is shown in fig. 2(a), and the volume occupied by the crystals in fig. 2(b). The present data were taken before the complete sedimentation of zone 4, but samples observed 3 month later showed no quantitative changes in their lower part (zones 1-3). Due to further crystallisation, p-3 the upper bound of the "ice pack" was higher and the liquid extended only ∼ 10 gravitational length before the gas-liquid interface. This is consistent with simple C+P mixtures [5].

The most striking feature of the multi-zone structure is that the compositions of the coexisting phases (φ L = 0.22 ± 0.02, φ S = 0.34 ± 0.02, φ X = 0.56 ± 0.02) are almost constant over zone 2. We note that the macroscopic vertical extent of zone 2, (∼ 2 mm), is thousands times of the colloidal gravitational length. This indicates that the two coexisting phases, XS and the binary fluid have the same density, otherwise they should be separated under the gravitational field. This is confirmed graphically on fig. 2(a) by the correspondence between φ L + φ S and φ X .

Phase behaviour. -Now we consider how such floating crystals may be rationalised, taking a binary hard sphere mixture (BHS) as a minimal system exhibiting this phenomenon. The more involved C+C+P, whose essential features are similar for these purposes to BHS, due to the very low polymer concentration in the colloid-rich phase [6], will be presented elsewhere. However, we note that the inclusion of polymers tends to destabilise binary crystals, which we neglect here. Furthermore, no binary crystals were observed in the experiments. Following [17], in fig. 3 we focus on the binary fluid (index A), the two possible one component crystals (index X) and their respective equilibria, including the three-phase eutectic region.

Let us consider a monodisperse fluid-crystal coexistence of hard spheres of diameter σ S = 1. The buoyancy difference between the two coexisting phases is proportional to their volume fraction difference ∆φ = φ X -φ A > 0. If we add one sphere of diameter σ L = 2σ S in the fluid and keep the osmotic pressure balance against the same crystal, we have to remove one small sphere to keep the same number density. But the large particle occupies more volume than the small particle that was removed so the volume fraction in the fluid increases. So ∆φ decreases with the introduction of large particles in the fluid. This effect should lead ultimately to ∆φ = 0 for a given φ * L . If both colloidal species are the same material, both phases have the same buoyancy, just as a zeppelin floats in the air. This is a single point on the phase diagram because φ * S and φ * X are fixed due to coexistence relations: the "floating point". For a crystal of large particles coexisting with a binary fluid, A L+S ↔ X L , adding smaller particles to the fluid enhances ∆φ, consistent with [25]. More quantitatively, it is possible to integrate a binary hard sphere fluid equation of state (here Hansen-Goos-Roth [16]) and a hard sphere crystal equation of state (here Alder-Wainwright [26]) to obtain the reduced excess chemical potentials in both phases [17]. For the fluid we get

µ ex A,L (φ L , φ S ) = φ 3 1 -φ 0 + 5 -φ 3 1 -φ 0 (φ 1 + φ 2 ) - 2(1 -3φ 0 ) φ 0 (1 -φ 0 ) 2 φ 1 φ 2 + 3 1 -3φ 0 + φ 2 0 φ 0 (1 -φ 0 ) 2 φ 2 1 + 2 -5φ 0 + 6φ 2 0 -φ 3 0 φ 2 0 (1 -φ 0 ) 3 φ 3 1 - (φ 0 -φ 1 )(φ 2 0 + 2φ 2 1 -φ 0 (φ 1 -2φ 2 )) (1 -φ 0 ) 3 ln(1 -φ 0 ), ( 1
)
µ ex A,S (φ L , φ S ) = ξ 3 LS φ 3 1 -φ 0 + 5 -φ 3 1 -φ 0 ξ LS (φ 1 + ξ LS φ 2 ) - 2(1 -3φ 0 ) φ 0 (1 -φ 0 ) 2 ξ 3 LS φ 1 φ 2 + 3 1 -3φ 0 + φ 2 0 φ 0 (1 -φ 0 ) 2 ξ 2 LS φ 2 1 + 2 -5φ 0 + 6φ 2 0 -φ 3 0 φ 2 0 (1 -φ 0 ) 3 ξ 3 LS φ 3 1 - (φ 0 -ξ LS φ 1 )(φ 2 0 + 2ξ 2 LS φ 2 1 -φ 0 (φ 1 -2ξ LS φ 2 )) (1 -φ 0 ) 3 × ln(1 -φ 0 ). (2)
with

φ i = φ L + φS ξ i

LS

, i = 0...3, the weighted volume fractions. At A L+S ↔ X L coexistence, we have equality of the pressure and of the chemical potential of L in each phase:

µ ex A,L (φ L , φ S ) = µ ex XL (φ XL ) + ln φ XL φ L , (3)
p A (φ L , φ S ) = p XL (φ XL ). (4)
Starting from the known monodisperse coexistence, we integrate eqs.( 3) and (4) numerically using eq.((1)). This gives us the parametric representation of the XL crystallisation line in the (φ L , φ S ) plane. The XS crystallisation line is similarly obtained. Being interested in the respective buoyancy of the phases, we plot the phase diagram (fig. 3) in the (φ 0 , θ) plane, with φ 0 = φ L + φ S the total volume fraction and θ = φ L /φ 0 the volume proportion of L particles. Of course, θ = 1 in the crystal XL and 0 in XS. In order to display the coexisting phases in fig. 3, the melting lines of the two crystals were drawn as follows: the crystal compositions φ XL and φ XS are drawn with respect to θ in the coexisting fluid. Therefore, on the abscissa we see the respective buoyancies of the fluid and coexisting crystal. The crossing of the two curves indicates the floating point.

As expected, ∆φ is broadened by adding S in the A L+S ↔ X L coexistence. In the A L+S ↔ X S coexistence, adding L shrinks the buoyancy gap, leading to a floating point whose position depends on the size ratio ξ LS (see fig. 4). Moreover, for ξ LS > ξ max LS ≃ 0.464, the floating point is not located in the stable A L+S ↔ X S coexistence region, but unstable against further phase separation (X L crystallisation). This is the case (ξ SL = 0.57) our experimental system belongs to [see fig. Kinetics. -We now describe the dynamic sequence by which a sample in regime II reaches the final multizone structure. The initially homogeneous sample first undergoes a gas-liquid phase separation, induced by the depletion attraction. Then, droplets of colloidal liquid sediment. After 2 days, the sample is separated into a colloidal gas in the upper part of the sample and binary liquid in the lower part. We did not see any heterogeneous nucleation of crystals at the bottom, unlike regime I. We believe this is because the bi-dispersity dramatically suppresses nucleation of crystals in the bottom binary fluid. Without flocculation, particles in this fluid should sediment individually. After another 2 days, the higher sedimentation speed of the large particles leads to a zone at the top of the liquid, almost devoid of large particles. The sample is thus divided in three zones: colloidal gas at the top, liquid with few large particles in the middle, and a binary fluid rich in large particles at the bottom.

Above the sedimentation front of large particles, φ L is low so homogeneous nucleation of S crystals proceeds. Here we emphasize that this removal of large particles from a liquid phase of small particles due to sedimentation is crucial for the initiation of homogeneous nucleation of S crystals and the resulting inverted distribution of large and small particles. Reduction of frustration against crystallization is a necessary condition for frequent crystal nucleation to be realised [27,28]. Considering the buoyancy difference between a crystal nucleus of diameter ασ S and of volume fraction φ X and the surrounding S fluid (φ S ), the Peclet number of the nucleus can be expressed function of the Peclet number of a single particle S: P e X = α 4 (φ X -φ S )P e S . For φ X -φ S = 0.1, a crystal nucleus of a diameter larger than 3.4 σ S has P e X > 1. Thus, such crystals experience little Brownian motion and simply sediment quickly. However, when these crystals reach the sedimentation front of the large particles [fig. 5(a reduces, as the binary fluid has a higher total volume fraction and is therefore denser than the almost monodisperse fluid. The relevant Peclet number for the crystals is then P e X = 10 4 (φ X -φ 0 )P e S . If the fluid and crystal compositions are close to the floating point, φ X ≈ φ 0 , which means P e X ≪ 1, the crystals float. After another 3 days, the crystals piled up from the bottom to the top of the large sphere sediment [fig. 5(b)]. The crystals did not grow after entering the binary fluid zone and we conclude that XS and the binary fluid appear to "coexist".

When the pile of the floating crystals reaches the top of the large-sphere-rich sediment, crystal nucleation is still active in the upper monodisperse S fluid. So crystals continue to pile up above the large sphere sediment. In the monodisperse fluid, however, there is no mechanism which prevents further growth of the crystals, and thus they fill all the available space to form the "ice pack" [see fig. 2(c)]. This fluid-crystal phase separation finally stops when the monodisperse liquid layer reaches a stable depth of a few gravitational lengths, which leads to the final state.

Conclusions. -To sum, we found that the introduction of attractive interactions between binary colloid mixtures leads to complex dynamic interplay between sedimentation and phase ordering. This kinetic pathway is fundamentally different from that of sedimenting binary hard sphere mixtures, where sedimentation is followed by heterogeneous nucleation: crystallisation is followed by sedimentation. Despite the complex nature of this multicomponent C+C+P system, we are able to rationalise many of our findings by appeal to previous C+P work, due to the low overall concentration of the larger colloids, and to binary hard spheres, due to the low polymer concentration in the colloid-rich phase after phase separation. The dynamic sequence we have observed leads to the formation of a multi-zone structure, including an exotic floating crystal zone and an ice pack zone. This final multi-zone state is a metastable state, likely connected to a hidden floating point. Our study indicates that the introduction of phase ordering in sedimentation leads to a rich variety of the final zone formation, which may contribute to the deeper understanding of zone formation industrial applications and the natural world.

Fig. 1 :

 1 Fig. 1: (color online) (a) Schematic phase diagram of C+C+P showing the three regimes considered here, I, II and III. Long dashed line represents the values accessed by changing volume only. Systems initially in the shaded regions II eventually reach triple coexistence [5]. (b-f) Final states without floating crystals, corresponding to regimes I (φ 0 P = 0.38) (b,c,d) and III (φ 0 P = 0.73) (e,f) in (a). Here (b) is fluid, and (c,d) are crystalline. Large particles are embedded in the crystal, a higher density towards the bottom (d) and reducing rapidly to a much lower density (c), reminiscent of a 'sedimentation front'. (e,f) are sedimented gel with limited local crystallisation and trapped polymer-rich regions (black). Approximate heights of images are indicated on the vertical axes. Gravity points downwards.

  )]. Experimentally, by changing the polymer concentration from φ 0 P ∼ 0.40 to ∼ 0.43, we cross the boundary between regime I and II. Note that in the underlying S+P phase diagram [5], the polymer volume fraction of φ 0 P ∼ 0.40 is required for triple coexistence. In regime III, with a high polymer volume fraction (φ 0 P ≥ 0.73), we observed a (XS) gel with some evidence of local crystallinity coexisting with a colloidal gas, [fig. 1(e,f)].

Fig. 2 :

 2 Fig. 2: (color online) (a) Sedimentation profile of a C+C+P sample (regime II) during sedimentation, with its bottom part already in the final state. Symbols represents the measured local volume fractions of each component in each phase. Lines are the theoretical floating volume fractions. (b) The fraction of the volume occupied by XS. Overall composition:φ 0 L = 0.004, φ 0 S = 0.186 and φ 0 P = 0.45. Total size the sample: 3 cm (only the bottom is shown here). (c) gas (zone 5)-liquid (zone4) interface. (d) ice pack (zone 3). (e) floating crystals (zone 2). (f) glass with a few crystallites (zone 1). Bright (red) regions are XS. Very bright (green) dots are large particles.

Fig. 3 :Fig. 4 :

 34 Fig. 3: (color online) Theoretical phase diagram of binary hard spheres for ξSL = 0.45 and ξSL = 0.57. Melting curves are represented by continuous lines as a function of θ = φL/φ0 in the coexisting fluid. Dashed curves are prolonging phase boundaries and melting curves unstable to further phase separation.

Fig. 5 :

 5 Fig. 5: (color online) The kinetics of phase ordering and sedimentation in a phase separating C+C+P system, regime II: (φ 0 P = 0.36) (a) Sedimentation front of large particles in the fluid. Time-lapse (every 10 min) picture showing crystals of small particles falling into the binary fluid. Total time = 3 days, total image width 0.1 mm. (b) The same place in the final state (one month), covered by the ice pack.

We thank Jeroen van Duineveldt for helpful discussions. ML thanks MEXT for a scholarship. CPR acknowledges the Royal Society for financial support. HT acknowledges a grant-in-aid from MEXT.


Dr Phillip Kent
email: p.kent@ic.ac.uk

The Visibility of Models: Using technology as a bridge between mathematics and engineering

Phillip Kent* (Imperial College) and Richard Noss** (Institute of Education) University of London, U. K.

Submitted to the proceedings of the ICMI Study Conference, "On the Teaching and Learning of Mathematics at University Level", Singapore, December 1998.

ABSTRACT: Engineering mathematics is traditionally conceived as a set of unambiguous mathematical tools applied to solving engineering problems, and it would seem that modern mathematical software is making the toolbox metaphor ever more appropriate. We question the validity of this metaphor, and make the case that engineers do in fact use mathematics as more than a set of passive tools-that mathematical models for phenomena depend critically on the settings in which they are used, and the tools with which they are expressed. The perennial debate over whether mathematics should be taught by mathematicians or by engineers looks increasingly anachronistic in the light of technological change, and we think it is more instructive to examine the potential of technology for changing the relationships between mathematicians and engineers, and for connecting their respective knowledge domains in new ways.

Introduction

The statement that mathematics plays a central role in engineering and science is certainly true, and-in its bare bones formulation-certainly a truism. "Mathematical models" are everywhere, "modelling" is a central activity. However, in this paper we would like to ask in what ways is a model "mathematical", and, at the same time, in what ways is it scientific, or part of engineering? How might these different aspects be connected in the minds of learners and experts?

Traditionally, these questions have unproblematic answers. If "service" mathematics is essentially a set of tools whose workings need not be visible to the user, then the difficulty is simply one of teaching "the mathematics" and learning to apply it later.

The metaphor of application is ubiquitous. But what is it that is applied? And what, exactly, is it applied to? Further, if the toolbox metaphor is to be helpful, we need to have some idea of what different people will see when they look inside the box. Will they see the same thing, will it have the same structure? Will it have the same function?

The experience of engineering students entails more and more contact with sophisticated pieces of technology. For example, with the latest computer-aided design software for civil engineering it is possible not only to "build" structures such as bridges in the virtual space inside the computer, but also to test the integrity of a design against the effects of an earthquake. Underlying this computational power is a huge amount of invisible mathematics, and it is clear that technology is allowing students to use mathematics to an unprecedented degree-in the case of computeraided structural design, the most advanced numerical techniques for solving nonlinear equations become available at the press of a button, and with barely a mathematical equation in sight. In these circumstances, the future role of mathematics teaching for engineers is uncertain, especially since most of the mathematical methods which form the staple diet of traditional mathematics service courses are now themselves available effortlessly in a computer mathematics system such as Mathematica [START_REF] Wolfram | The Mathematica Book[END_REF].

In this respect, it seems that the computer is, if anything, making the toolbox metaphor ever more appropriate. If solving a nonlinear equation is a question of pressing the right buttons, it is not inappropriate to think of it as similar to selecting the right spanners -and we don't seem to need much instruction about how spanners work (or are designed) to use one. On the contrary, we might be forgiven for asking how the connections between mathematics and its applications in engineering can be made more visible by using computer mathematics software, which, it is commonly acknowledged, hides mathematics inside general-purpose, black box functions for doing integration, equation solving and the like?

This role of technology, together with (in the UK at least) the well-documented decline in mathematical preparation of incoming students (e.g. [START_REF] Sutherland | The Changing Mathematical Background of Undergraduate Engineers: A review of the issues[END_REF], [START_REF] Lms | Tackling the Mathematics Problem[END_REF]), has led some to make a reasonable case for the downplaying of the role of mathematics in engineering (see, for example, [4, p. 264]). There is, undoubtedly, an argument that significant kinds of engineering can be done with mathematics which has already been done by someone else, and wrapped up into computational tools which the engineer needs only to use.

In fact, it seems that this is a very partial view. There remains a strong case for the inclusion of mathematics as more than a set of passive tools, catalysed by the computer in new ways (see, for example, [START_REF] Wilkinson | Mathematics for engineers: an industrial view[END_REF]). In this paper, we will outline our case that the computer, if appropriately conceived, affords an opportunity to make visible important parts of the mathematical agenda, rather than to relegate mathematics into a set of tools whose workings remain opaque.

Our position demands attention to epistemology rather than merely technology and its application. Of course, we will need to consider technology-focussed issues, such as what is possible with a piece of software such as Mathematica-what can be done with in the context of a given mathematical or engineering topic. But we want to focus on issues to do with the basic relationships between mathematical and engineering knowledge. These are fundamental in our attempt to rise to the challenge of designing and structuring activities which simultaneously lead students to use and understand the mathematics they are deploying in their computationally-based activities.

In the UK, students specialise in their degree subject from the start. For this reason, if no other, they meet some demanding mathematics as soon as they enter university, and are called upon to "apply" it almost immediately. The example student activities in this paper are drawn from a short (6 hour) introductory course in Mathematica, for first year undergraduates in the Civil Engineering department at Imperial College, which was designed and delivered, for the first time, by the METRIC Project 1 in early 1998.

Developing a "structural feel" for beams and bridges

An introductory Mathematica course for undergraduates can easily fall into the class of generic software training; there is so much that seems to need to be discussed (symbolic calculations, numerics, graphics, programming, using the document interface), that it is easy to spend the whole time exhibiting the functionality of the software. With the Civil Engineers, however, we wanted to use the course to present activities to the students where Mathematica is being applied in specific engineering contexts, and we enlisted the help of a colleague from Civil Engineering 2 to develop these contextual examples 3 . At the very least, seeing Mathematica applied to relevant situations in engineering is likely to be good motivation for students, but our hope is that the "bridging" effects can be more significant than this.

The two contexts we chose were both to do with structures. The idea of the engineers was that, by letting Mathematica take the mathematical strain, we could help students begin to get a "structural feel" for how structures behave (something that their present courses seem to be deficient on). After the first run, we can claim to have existence theorems for this; we will be seeking more substantial evidence when the course runs again in the current academic year.

The first Mathematica session for the students was a quick overview of the numeric, symbolic and graphical capabilities of the software. As a final exercise, the students were invited to "apply" their fresh knowledge to a typical loaded beam problem, such 1 The project team is Phillip Kent and Phil Ramsden. See http://metric.ma.ic.ac.uk/ .

2 Dr David Lloyd-Smith, to whom we express our thanks.

3 METRIC has developed a similarly "contextualised" approach for Chemical and Mechanical Engineers, and for Chemistry students; for details, see the web site already mentioned.

as they meet in their initial engineering course on structures. The "structural feel" idea prompted an emphasis not on the mathematics of the problem-which is given to them in full-but on estimating important structural quantities as a load (W, below) is varied, using whatever combination of graphical, numerical and symbolic methods they choose. These quantities include the point of maximum displacement along the beam, and the point of "contraflexure" (i.e. where the curvature changes sign).

The students were presented with the situation shown in Figure 1, where the uparrows denote fixed supports, there is a distributed load between x = 3 and x =6, and W is a variable point load: The students were given the solution for small deflections of the beam. This is conventionally written by engineers in the form:

y x W 2 kN/m 3m 3m 4m 4m
EI y W x W x x x W x = - + - - - + - + + - 5
The constant term EI is typically around 10 7 so the deflections are very small, of order mm for y when x is of order m (Figure 2). The y-vs-x equation is not quite as it seems, because the polynomial terms ( ) x-3 4 etc, are written in the normal way, but in fact represent piecewise-defined ramp functions, defined to be zero when x < 3 and to be ( ) x-3 4 when x ≥ 3 , etc. (Figure 3). While this notational shift is implicit in much of the spoken and written language used with the students, it is made explicit, indeed it must be made explicit, in the Mathematica expression for y which the students are given:

y = (1/EI)*(5(-22275+4W)x/3 + (33000-4W)x 3 /60 - (250/3)*If[x<3,0,(x-3) 4 ] + (250/3)*If[x<6,0,(x-6) 4 ] + (27000+14W)*If[x<10,0,(x-10) 3 ]/60)
From the engineer's point of view, there is nothing strange; as our colleague put it: "of course, these terms here are ramp functions...". But for us it was surprising to discover something new about polynomials: basic and boring mathematical objects, but when looked at in a certain (engineer's) way, they are "ramps". Moreover, this is true in a dual sense, both as a visual metaphor, and as an expression of the role that the functions serve in the structural analysis: that the terms containing ramp functions have zero effect on the beam displacement until x reaches a threshold value.

We don't want to overstate the importance of a small episode, but it does point up the fact that mathematics is not a passive agent. In use, mathematics becomes a means of making sense of the underlying engineering principles. But reciprocally, the mathematics itself is shaped by its application-it takes on meanings which are derived from the setting in which it is used.

Incidentally, we discovered one other curious (to us) phenomenon: in situations like this, the beam's weight is often negligible in magnitude relative to the other forces in the problem. So the weight is abstracted into the form of a distributed load (pressure).

In effect, the beam is abstracted to an "ideal structure" defined only by its geometry, flexibility and material strength.

The Rainbow Bridge

The second Mathematica activity for the Civil Engineers is based on second-year mathematics material, and it represents an instance of didactical inversion: using the capability of the technology to allow students to carry out some task using mathematics which the students don't know yet in order that they can focus on some conceptual points which the mathematics makes accessible. Two pre-written Mathematica functions generate animations of a test load moving across two different simple bridge structures; at each step in the animation, the colour of each of the struts in a bridge represents the magnitude of the force in that strut induced by the test load (Figure 4). Again in the interests of "structural feel", the students are not required to understand the Mathematica details of how the animation is generated, nor the mathematical details of how the strut forces are calculated (which involves the solution of systems of linear equations). What they do have control over are the magnitude of the test load and the "colour function" which maps a numerical force value onto a range of output colours. They are asked to consider how to design a colour function which yields the most useful information about what is going on in the bridge as the test load moves across it, and to design a function which would allow them to detect the maximum safe load that can cross a bridge given a maximum safe force for any strut.

It should be clear what are the engineering lessons from this activity: the students can get experience in how the patterns of forces vary in a loaded structure, and they are invited to consider, albeit for a toy example, a central engineering design question of determining what loads a given structure can safely support. The mathematical lessons may not be so obvious. Indeed, one might ask, where is the mathematics at all? Haven't we hidden all the relevant mathematics inside the Mathematica functions? In the most obvious sense, we have hidden the mathematics of the problem-the solving of systems of linear equations. But in fact, the mathematics can be made visible in two ways. First, the didactic inversion allows us to hide the details of the mathematical processes whilst keeping visible the very useful results of that process. Second, the students are invited to engage in an interplay between bottomlevel Mathematica programming-defining colour functions-and high-level visualisation. For example, the "maximum safe load" question demands some kind of piecewise-defined function along the lines of (taking 35000 N as the max. safe force in a strut): This process of making and criticising representations (i.e. the colour function mappings) is not conventionally recognised as mathematical work, at least not for non-expert, beginning students. But, we think that it reflects a kind of mathematical thinking that has a great deal to do with having a good structural feel. (And we think that it has a strong relationship with diSessa et al's [START_REF] Disessa | [END_REF] idea of "meta-representational competence").

These structures activities highlight a problem of visibility in design. They illustrate the complexity of the questions we asked at the outset: for now it should be clear that the question is not only whether or not to make the mathematics visible, but what mathematics, and in what form, to make visible? Designers must choose to make certain pieces of the mathematics visible and functional, and it is functionality which is the very real contribution of the technology. At the same time, students have to map the expression of the mathematics (in whatever form) into the results they see-and to try to re-represent those results in terms of the (Mathematica-based) mathematics. This is not simply a matter of multiple representation, it is a matter of construction.

Discussion: Designing for visibility in a mathematical software environment

In this section, we want to consider the visibility of mathematical calculations in three different software packages-Mathcad, Mathematica and Maple-as well as how, and to what degree, different teachers of engineering mathematics are choosing to make mathematics visible whilst using those packages with students. Insight into the latter was gained from the proceedings of a recent workshop on the use of mathematical software packages in undergraduate engineering education 4 . This happened to allow us some rather intriguing views on the relationship between epistemology and visibility-in other words, how the intentions of the designer are translated into the conceptual mathematical models developed by the user/learner. As we shall see, the relationship is not straightforward.

Clearly, all the various software manufacturers are interested in appealing to as large an audience of mathematics users as possible, and their "box top" slogans express this: Mathcad-"the worldwide standard for technical calculations"; Mathematica-"the world's only fully integrated technical computing system"; Maple-"complete mathematics and visualisation system". But if one looks inside the box, we think that the different epistemologies of engineers and (applied) mathematicians can be made out in the software designs.

Mathematica and

Maple are examples of "computer algebra systems", and represent what an applied mathematician might expect of "computer mathematics": comprehensive sets of symbolic, numerical and graphical functions, expressed in a precise, extensible mathematics-like programming language. (It should be said that neither are much up to doing formal proofs automatically, but then again they don't claim to be replacing human mathematical thought-fortunately!-but to provide an environment to support it.)

Mathcad is a package very much designed for, and commercially targeted at, engineers. It works rather like a "sketchpad" combination of word processor, spreadsheet and mathematical (mostly numerical) toolbox: inputs and results can be placed quite freely on the screen/page, but they are causally connected behind the scenes. It is interesting to trace the evolution of Mathcad over its past three or four versions. As usual, more functions, menus and palettes have appeared, but a couple of developments seem more indicative of a particular design philosophy: the first is the way that Mathcad's developers are acquiring the electronic rights to many of the standard engineering data books, and making them available as $200 "electronic library" add-ons to the basic system. Mathcad is data-oriented, and proud of it.

Second, there is a "symbolic toolbox", that performs a selection of symbolic algorithms, which has grown in mathematical coverage with each new version (in fact, it is a portion of the Maple "mathematics engine" running in the background). We presume that this growth depends not least on the fact that enough users have requested a particular symbolic function to be added. Also, presumably, the developers of Mathcad have to pay the developers of Maple more to use more symbolic functions, which implies a certain conservatism on the part of the former. Now, from the point of view of a typical Mathcad user, this growth process must seem quite natural. An engineer is faced with a problem to solve, and needs to apply mathematical techniques to solve it; chances are it won't yield to a symbolic technique anyway (few mathematical equations of practical use do possess exact analytic solutions outside of special cases), but having an improved package of symbolic techniques to hand is going to turn out to be useful some of the time.

However, from the viewpoint of a mathematician, this haphazard growth process could seem pretty worrying. Mathematicians make strong distinctions between symbolic and numerical procedures, so much so that the latter are often treated as a separate field of the discipline (i.e. numerical analysis). This is natural, too; surely it's the business of mathematicians to make such distinctions?

The message that we take from this comparison of perspectives is that visibility is not a simple issue: it depends on what designers, and users, think "mathematics" is. While it is no surprise that engineers and mathematicians see the function of mathematics differently, it is perhaps more surprising that they may not be thinking about the same mathematics: if that is true, it raises some difficult questions about the nature of applied mathematics itself, and surely indicates that the metaphor of application is, at best, limited.

Should it be the business of engineers to make the same distinctions as mathematicians-to work with the same mathematical epistemology? In particular, do engineering students get the most appropriate mathematical training by following traditional mathematics courses which give pride of place to symbolic techniques, and relegate numerical methods to second place? For example, a speaker at the workshop, a teacher of civil engineering students, declared that a particular bugbear of his is having to re-orientate students who have been taught in school that integration is firstly about backward differentiation (a symbol-oriented view), and secondly (if at all) about determining areas under curves.

In principle (institutional finances and academic politics notwithstanding), a lecturer in engineering mathematics can choose between offering students a package like Mathcad, or a package like Mathematica (or Maple). Mathcad has the advantage of being a tool tuned for engineers, whilst Mathematica may be less so (it is certainly difficult to get to grips with it if you're not willing to think in explicit mathematical terms).

Either way, the educator has to come to terms with the design choices that the software developers have made: which mathematical aspects are visible enough, which need to be made more visible, or indeed less visible (as in our rainbow bridge example above).

It was interesting to note the high degree of unanimity at the workshop-amongst users of Mathcad, Mathematica, Maple, Derive and other systems-that the most important issue is how to develop competence in the use of appropriate mathematical technology, which all agreed was grounded in a combination of traditional and technologically-based mathematical understanding. One of the Mathcad users put it like this: I let [the students] work in the way that an engineer would: they tend to work backwards to the mathematical questions from Mathcad's answers: what do these things mean? If the answer isn't what you expected, did you expect a positive number, or a real number? Can you find by hand a bracketing answer? Then you've got confidence in the computer result.

And this from a Maple user:

What we need to ensure is that our students become failsafe users of these powerful mathematical tools. The problem is very similar to (but more complex than) that which we have with the uncritical use of calculators to compute unlikely numerical answers. … With computer algebra systems the problem of validation is much greater and the techniques of validation are more varied and more sophisticated. …. What rules or heuristics do you teach to your students to ensure that they become discriminating users of computer algebra systems and failsafe engineers? However, as we already pointed out in the introduction, there is certainly not unanimity in the engineering community as a whole about the way to do, and to learn, mathematical modelling:

We encourage students not to be linear toward the answer, and to go back around the loop to understand things better, but when the students graduate into industry they're usually told that "academic approaches" can't be tolerated, "here in industry we're results-oriented", it's a case of getting the answer and then moving onto the next problem. The point is to know when you've done enough to be sure of your answer.

Summing up, it is clear that the epistemological decisions built into software design far from determine the user's activities. Epistemological structures shape and are shaped by what the user does, but these are not straightforwardly linked.

Conclusions

We have presented our examples of mathematics teaching and educators' discussion with the aim of challenging the traditional view of mathematics: that it is either studied in its own right or must inevitably be viewed as a succession of recipes, preferably wrapped in computational dressing. The former view may be attractive to mathematicians, but it has consistently failed engineering and science students. On the other hand, it seems increasingly likely that the latter view will render invisible crucial parts of the scientific and technological endeavour, in ways which relegate mathematics only to the privileged few who design the programs. This is, we think, an increasingly problematic issue, and one which is facing all those whose work involves-implicitly or explicitly-mathematical knowledge and techniques.

Michael Clayton [START_REF] Clayton | Industrial applied mathematics is changing as technology advances: What skills does mathematics education need to provide?[END_REF], a mathematician working in the multidisciplinary environment of the telecommunications industry, has pointed to the "bridging" effects of technology on the relationships between mathematicians and engineers in industrial practice, overturning the traditional roles of mathematicians as makers of models, and other people in the design and production process as consumers of models: General-purpose IT tools such as spreadsheets, and mathematically based environments and workbenches such as Mathcad and Matlab have made it easier for engineers, dealers, salesmen, managers and others to construct their own models and refine them for specific applications. When time is of the essence, the value of these tools lies in the rapid prototyping they allow: initial modelling ideas can be investigated by the potential users, and the resulting interaction often leads to an improved match between the model and the users' requirements. Modern graphical user interfaces … can be designed to make even the most sophisticated special-purpose models accessible to the people who need to use them, helping to remove the "ivory tower" and "back room" images that have sometimes been attached to mathematicians in the past. [7, p. 25] Clayton's insight may be crucial for effective university mathematics teaching in the future. The perennial debate over whether mathematics should be taught by mathematicians or by engineers looks increasingly anachronistic in the light of technological change, and modern industrial working practice. We think it is more instructive to examine the potential of mathematical technology to change the relationships between mathematicians and engineers, and to connect both people, and the knowledge domains in which they work, in new ways.

The tools we use, as much as the activities we design, shape the kinds of understandings our students construct. Moreover, the mathematical models for phenomena-however straightforward they are to mathematicians-are not straightforward at all: they depend critically on the settings in which they are used, and the tools with which they are expressed. Provided we are explicit (at least to each other, and perhaps to our students) we see this as a mathematical opportunity: in contrast, leaving this issue (and the mathematics) invisible must, we think, be a source of difficulty.

Figure 1 :

 1 Figure 1: The loaded beam problem in Civil Engineering.

Figure 2 :

 2 Figure 2: Graphs of beam displacement (for W ranging from 0 to 4 kN).

Figure 3 :

 3 Figure 3: The ramp function " (

Figure 4 :

 4 Figure 4: Frame from a "bridge movie"; a test load (black disc) moves across a (2-D) bridge, whose struts change colour according to the magnitude of the force caused by the test load (the numbers just index the struts for reference).

  overloadFun[force_]:=If[Abs[force]>35000, GrayLevel[1], (*safe load exceeded -output white*) Hue[0.8*Abs[force]/35000] (*else output a colour spectrum*) ]

Organised by METRIC at Imperial College, June 1998. Proceedings are available from http://metric.ma.ic.ac.uk/symposium/ .


Koshi Hasatani

Mathieu Leocmach

Anthony J Genot

André Estévez-Torres

Teruo Fujii

Yannick Rondelez

High-throughput and long-term observation of compartmentalized biochemical oscillators

come

High-throughput and long-term observation of compartmentalized biochemical oscillators

Koshi Hasatani, 1 Mathieu Leocmach, 2 Anthony J. Genot, 1 André Estévez-Torres, 3 Teruo Fujii, 1 and Yannick Rondelez We report the splitting of an oscillating DNA circuit into ∼ 700 droplets with picoliter volumes. Upon incubation at constant temperature, the droplets display sustained oscillations that can be observed for more than a day. Superimposed to the bulk behaviour, we find two intriguing new phenoma -slow desynchronization between the compartments and kinematic spatial waves -and investigate their possible origins. This approach provides a route to study the influence of small volume effects in biology, and paves the way to technological applications of compartmentalized molecular programs controlling complex dynamics.

Recent progresses in molecular engineering have allowed the construction of synthetic biochemical systems that display subtle dynamics such as oscillations or multistability [1][2][3][START_REF] Padirac | Proc. Natl. Acad. Sci[END_REF]. The building of such molecular programs from scratch yields a unique glimpse of the challenges and locks that evolution had to overcome in its drift towards more and more complex life forms. Moreover, the inherent biocompatibility of these chemical circuits supports their use to monitor and control biological processes.

Compartmentalization of these man-made dynamic systems would offer tantalizing additional possibilities. Cell-sized compartments (∼ 1 -100 µm) may provide a model of biotic or prebiotic organization [START_REF] Szostak | [END_REF]. For example, stochastic effects arising from low numbers of molecules in small spaces are very important in biology [6,7], but hardly explored outside of living systems. Synchronization of many autonomous elements using quorum sensing provides another example [8].

Beside these basic biological motivations, microcompartmentalization is also the obvious way forward for the exploration of man-made molecular circuits. First, it permits the running of hundreds of circuits using volumes that would have yielded a single experiment otherwise. Such high-throughput will be required to tune the many parameters controlling the behavior of molecular assemblies. Even in the case where all the compartments possess an identical "program", compartmentalization may yield new analytical concepts such as digital PCR [START_REF] Vogelstein | Proc. Natl. Acad. Sci[END_REF]. Second, splitting -and later establishing controlled connections between -spatially distributed molecular circuits may unleash the computing power of molecular programs [START_REF] Toiya | [END_REF] by removing cross-talks and allowing the reuse of modules.

Yet, compartmentalization and long-term monitoring of biochemical reactions in micro-compartments has been difficult to achieve. Indeed, the reactions involved must run sustainably in closed systems and the potentially detrimental effects of large surface/volume ratios and leaks need to be tightly controlled.

While biochemical systems with a simple dynamic such as qPCR were monitored in droplets [11], out-ofequilibrium biochemical circuits pose distinct and additional challenges. PCR is driven sequentially by an external operator, uses only a single enzyme and lasts rarely more than 45 minutes. By contrast, our dynamic systems display oscillatory behaviour, are fully autonomous, require 3 enzymes, and last for days.

More importantly, sustainability of oscillations depends on a much more delicate balance between reagents, buffer and temperature than PCR. Most biochemical assays typically consist in a simple relaxation toward a stable steady state, but oscillations arise from the destabilization of all steady states, which requires tight control over reaction parameters.

Here we report the compartmentalization and day-long tracking of 700 oscillators into picoliter droplets. This represents an increase of throughput by two orders of magnitudes and a decrease of volume by 4 orders of magnitudes compared to the literature on synthetic DNA sys-tems1. We show that with an experimental setup minimizing evaporation, droplets act as independent chemical containers that satisfyingly reproduce bulk conditions. Moreover, compartmentalization also reveals two striking phenomena that may remain hidden in bulk: slow desynchronization and kinematic spatial waves.

We use here a recently reported synthetic biochemical oscillator [12]. While simple to prepare and well characterized, it involves subtle enzymatic kinetics, which should put a stringent test on the use of droplets as compartments. The biochemical system reproduces the ecological predator-prey mechanism: molecular preys catalyse their own replication, but also serve as fuel for the replication of their molecular predators. Both prey and predators are continuously degraded by an exonuclease. This system results in robust oscillations, which are monitored via the fluorescence intensity of a dye bound to the DNA template encoding the circuit [13].

Droplets were generated by a flow-focusing junction fabricated in a PDMS chip. For each run, we generated thousands of droplets from the same reactive mix. The droplets were then transferred to a chamber, which was formed between two glass slides sealed by araldite. The chamber was kept at constant temperature (45.5 • C) by placing it under the heated glass plate of a microscope stage equipped with a temperature controller (Fig. 1). A field of view contained about 1000 droplets (movies in ESI). While the oscillations lasted for 2 days, we restricted the analysis to one day in order to increase the number of trajectories successfully tracked [14][15][16].

Fig. 2a shows a representative field of view and time traces for the fluorescence of droplets. Most droplets are monodisperse, with a diameter of ∼ 100 µm and a dispersion of ∼ 15%. This corresponds to a volume of ∼ 500 picoliters, which is 4 orders of magnitude smaller than bulk reaction volume (20 µl). A few droplets are significantly larger, which may result from coalescence during the generation and transfer of droplets. The fluorescence of the majority of droplets oscillates, and the oscillation of each individual droplet is similar to the bulk control.

The collective behaviour was more surprising: spatial waves move across the chamber while slow desynchronization ultimately sets each droplet's fluorescence on a seemingly independent trajectory (Fig. 2b). Since the contents of all droplets derive from the same mix, we expect them to be initially synchronized. Indeed, noting T i the time required to reach the i th peak (Fig. 2a), we measure a standard deviation of 7 minutes for T 1 . After 15 hours, the standard deviation of T 7 had increased to 35 minutes (25% of the mean period). At this point no coherence can be visually detected in the droplet population. We tested if this desynchronization was specific to compartmentalization or also existed in bulk. In 20 µl tubes, we observed a similar albeit less pronounced desynchronization (S8, ESI). This common desynchronization points to the influence of an external factor such as temperature on the period of oscillations. We characterised the sensitivity of oscillator to temperature in bulk (S9, ESI) and found that an increase of temperature of 2 • C lengthens the period by 25%.

Desynchronization is not spatially random, but manifests itself initially as travelling waves (movies M1 and Fig. 2a). One may misinterpret spatial waves as the result of diffusive transport of reagents. However, pseudo (also called kinematic) waves may appear when properties determining phases or periods vary spatially [17]. In view of the temperature-sensitivity of the oscillator, a gradient of temperature could not only desynchronize the droplets, but also generates kinematic waves.

The kinematic wave model predicts in particular a wave speed inversely proportional to the wave number, which we observed experimentally (Fig. 1c,S11 ESI). The kinematic model also implies a thermal gradient of 0.3 K mm -1 , which agrees with the measured gradient of 0.15 K mm -1 (S10, ESI).

On the contrary, the experimentally observed speeds of the waves do not support a model where waves are caused by diffusion of DNA strands between droplets. This model would be consistent with a time for dropletto-droplet diffusion of ∼ 5 minutes (S11,ESI).

This timescale is 2 orders of magnitude larger than the diffusion of a smaller molecule (fluorescein) between droplets stabilized by a more permeable surfactant (Span-80) [18]. The timescale is also incompatible with the possibility of digital PCR in droplets (which lasts 30 minutes) [19].

To further test that diffusion between droplets was negligible, we prepared two populations of oscillating droplets with a delay between them. One population was pre-incubated at reaction temperature (45.5 • C) for approximately half a period, thus allowing the oscillation to progress. The remaining population was kept at room temperature where reactions are much slower. The two populations were then mixed and incubated together on stage. If diffusion between droplets is negligible, the delay between the two populations of droplets should persist over time. The movie M2 shows the sustained presence of two populations. Fig. 3b shows the distribution of time-to-peak for this experiment. The initial delay is clearly observable from the shape of T 1 , where the populations are easily distinguished. This bimodal distribution was maintained at least until T 6 despite the widening of both populations (overall standard deviation of 16 minutes at T 1 versus 51 minutes for T 6 , compared to 7 and 29 minutes for T1 and T 6 (in Fig. 3a).

This population stability supports the functioning of each droplet as an independent reactor. The kinematic waves suggest that droplets are essentially identical but placed in a thermal gradient. Yet, some droplet-specific desynchronization is also visible in the movies (Fig. 2d) and we cannot rule out that intrinsic factors such as stochasticity or partitioning noise play a role in the ob-served dynamic. Predator-prey oscillators are nonlinear and their autocatalytic loops may amplify minute differences or fluctuations between droplets. Noise-induced stochasticity is an important determinant of biological circuits dynamic.

In conclusion, the picoliter droplet-based compartments used here satisfactorily reproduce the behaviours observed in bulk, in spite of reducing volumes by a factor of 10 4 . Therefore this approach may open the highthroughput experimental exploration [20] of the range of dynamics displayed by synthetic molecular programs when varying their control parameters. Another interesting direction would be to engineer a form of controlled molecular communication between the droplets in order to trigger and tune the onset of collective behaviours4. A further shrinking in size would better match the dimensions of biological cells and provides a unique platform to study the impact of small molecule numbers on the dynamics of dissipative reaction networks.

Compartmentalized oscillators also offer an in vitro model to investigate pattern-formation in biology. Kinematic waves clearly illustrate that spatio-temporal patterns -ubiquitous in biological development-need not arise from a reaction-diffusion mechanism, but may also originate from a spatial gradient of parameters. Purposely tailoring spatial gradients of parameters will offer a route to create novel classes of patterns.

This research is supported by the CNRS (France) and the JSPS (Japan). We thank V. Taly, H. Guillou and A. Padirac for advice as well as B.J. Kim and R. Ueno for using their thermal camera. During the course of our research, E. Winfree and F. Simmel kindly informed us of

FIG. 1

 1 FIG. 1. a) Water-in-oil droplets are generated inside a PDMS chip with a flow-focusing junction. b) Droplets are then transferred to a glass chamber sealed with araldite for incubation and observation. c) The chamber is placed under a heated stage and observed by fluorescence microscopy.

FIG. 2

 2 FIG. 2. a) Microscopy images showing the first kinematic spatial wave. The arrow indicates the travelling direction (colours reflect the intensity of fluorescence). The temperature is set at 45.5 • C. b) Top. Superimposed fluorescence traces of all successfully tracked droplets (690 traces). The dashed box indicates the peak corresponding to the wave shown in a). Bottom, bulk fluorescence of the oscillator (set at the same temperature) measured in a 20 µl tube in a qPCR machine. c) Wave speed as a function of the wave number, fitted with a kinematic wave model (S11, ESI). d) Kymograph showing the desynchronization of a column of droplets lying on the same initial wavefront (column indicated by a vertical line in a)).

1 1

 1 LIMMS/CNRS-IIS, University of Tokyo, Komaba 4-6-2 Meguro-ku, Tokyo, Japan 2 Institut Lumière Matière, UMR5306 Université Claude Bernard Lyon 1 -CNRS, Université de Lyon, 69622 Villeurbanne, France 3 Laboratoire de photonique et de nanostructures, CNRS, route de Nozay, 91460 Marcoussis, France

their ongoing work on the effect of compartmentalization and stochasticity in biochemical oscillators.