Some materials, such as coffee, milk or marble, have a soft translucent aspect because of sub-surface scattering: light enters them, is scattered several times inside before leaving in a different place. A full representation of sub-surface scattering effects in illumination simulation is computationally expensive. The main difficulty comes from multiple scattering events: the high number of events increases the uncertainty on the result, forcing us to allocate more time for the computations. In this paper, we show that there is a strong correlation between the surface effects of multiple scattering inside the material and the effects after just two scatter events. This knowledge will help in accelerating multiple scattering effects. We exploit this knowledge to provide a model and implementation for fast computation of double-scattering events, using a precomputed density function stored in a compact way. 1. Material: σ t = (0.164, 0.26, 1.192), α = 0.95, η = 1.3, g = 0.3. 1

Introduction

Many materials, such as milk or wax, exhibit sub-surface scattering behaviour: the light that reaches them enters inside, and is scattered several times before leaving again. This effect gives the material a slightly translucent look. The overall aspect depends on how far the light can penetrate inside the material, and on how many times it will be scattered before the exit. This effect is visible on solids such as marble and jade, as well as liquids such as milk, coffee or orange juice. It results in a softer look for objects made with these materials.

We have to model and render this sub-surface scattering effect for photorealistic rendering of these materials. But taking it into account greatly increases the computational complexity of illumination simulation. In order to compute the outgoing light at a specific point, we now have to take into account the incoming light at all neighbouring points, from all directions, instead of just the light incoming on this point. This adds two dimensions to the sampling, increasing the computation time. We also have to store the behaviour of the material, a function that express the relationship between incoming and outgoing light. Because the relationship is both spatial and angular, this function, the BSSRDF, has 6 dimensions. With regular sampling on all directions, storage is prohibitive.

A full representation of sub-surface scattering effects is expensive both in storage and computation time. To reduce the complexity, we use approximations, such as the diffusion dipole. For materials with a small mean-free path, light encounters a large number of scattering events before exiting the material. As a consequence, outgoing light is diffuse, removing the angular dependency and reducing the number of dimensions. Although is is practical and efficient, this approximation is not sufficient to represent the full range of sub-surface scattering materials. In a seminal study, Donner et al. [DLR + 09], proved that materials can exhibit strong directional response, with complicated lobes. They also provided a detailed and compact representation for BSSRDF lobes in the angular dimension, but still sample regularly in the spatial dimension.

In this paper, we build upon the work of Donner et al. [DLR + 09]. We focus on multiple scattering effects (light being scattered several times before leaving the material). We reproduced the experimental study of Donner et al. [DLR + 09], but separate between paths with two scatter events and others. Our study shows a strong correlation between the light leaving after two scatter events and the light leaving after multiple scatter events, for a large range of materials. As the complexity of paths with two-scatter events is much lower than the complexity of multiple-scatter events, we can use this knowledge to speed up scattering computations. In a second contribution, we provide a fast, low-memory algorithm for computing two-events sub-surface scattering and extrapolate for multiple scattering.

Our paper is organised as follows: in the next section, we briefly review relevant work on sub-surface scattering. In section 3, we present our study, comparing two-scatter events and multiple scatter events. We show there is a strong correlation between two-scatter events and multiple-scatter events, except for strongly backward scattering materials. In section 4, we present our model for compact representation of two-events sub-surface scattering; the probability density of twice-scattered photons can be represented independently of the incoming direction. We then use this probability density to compute the outgoing light by integrating along the ray. Finally, in section 6, we conclude and present directions for future work.

Previous work

We can do illumination simulation with sub-surface scattering effects by tracing photons inside the material, taking into account every scatter event [START_REF] Kajiya | The rendering equation[END_REF]. This technique is highly expensive in computation time, resulting in many research for faster representation of scattering effects.

Single scattering events are complicated. When the camera and viewer are assumed to be at an infinite distance, Hanrahan and Krueger [START_REF] Hanrahan | Reflection from layered surfaces due to subsurface scattering, 20th annual conference on Computer graphics and interactive techniques[END_REF] showed that there single-scatter events can be represented as a BRDF. In the general configuration, single scatter events are highly directional and can result in caustics inside the material [START_REF] Walter | Single scattering in refractive media with triangle mesh boundaries[END_REF]. Subsequent research has focused on the multiple scattering problem: light being scattered several times before leaving the material.

When there is a large number of scattering events, exiting light has lost all directional information. It can be modelled using a diffuse representation, the dipole approximation, introduced by Jensen et al. [START_REF] Jensen | A practical model for subsurface light transport, 28th annual conference on Computer graphics and interactive techniques[END_REF]. When there is a small number of scattering events, the outgoing light keeps directional information. We need a more complex model to represent it. D'Eon and Irving [START_REF] Eon | A quantized-diffusion model for rendering translucent materials[END_REF] greatly increased the accuracy of the dipole model, using a better approximation for the diffusion inside the material. Their approach still loses directional information, while we keep it. The difference could be important for anisotropic materials.

Donner et al. [DLR + 09] conducted a very thorough study of these sub-surface scattering materials. They worked both on Monte-Carlo simulations and measurements from actual materials. They showed that the material response has a strong directional component and that the shape of the lobes can be complex. They also provided a model for the directional component of the material response, using elliptic coordinates. But they did not provide for the spatial response, and sampled in concentric circles around the point of impact. Their model is the most accurate to date, but it has a high memory cost (up to 250 MB for each material). It is also limited to a single lobe. This paper extends and completes their study. We also provide a model that is less accurate, but much more compact and can represent multiple lobes.

3

Analysis of subsurface scattering

Properties of translucent materials

As light enters a translucent material (see Figure 1), it is first refracted at the interface using Snell's law: the light beam changes direction, while staying in the plane defined by the incident light beam and the normal at the entry point. The angle of the refracted ray with the normal depends on the index of refraction of the materials:

η sin θ i = sin θ i
Where η is the ratio of the index of refraction. If the first material is empty or almost empty (such as air), η is equal to the index of refraction of the second material.

After being refracted at the interface, light travels inside the translucent material. Here, it can encounter particles, resulting in either absorption or scattering. If it is scattered, light changes directions and continues travelling. As we are interested in the outside appearance of the material, we keep track of light until it reaches the surface again. The effect of a translucent material on light can be modelled using a few parameters:

• The absorption coefficient, σ a express how the material absorbs light. It is the reciprocal of the average length travelled by a photon before being absorbed.

• Similarly, the scattering coefficient, σ s express how often the rays are scattered inside the material. Both coefficients have units in m -1 . • The intensity along a ray decreases with the distance travelled, d:

I(d) = I 0 e -(σ s +σ a )d
(1)

We call = (σ s + σ a ) -1 the mean-free path. It corresponds to the average length a photon can travel before getting scattered or absorbed. It is expressed in m. We call σ t = σ s + σ a = 1/ the transmission coefficient. • When a ray is scattered inside the material, the phase function describes the angular distribution of light after scattering. In this work, we used the Henyey-Greenstein phase function, which has rotational symmetry and depends on the angle between the direction of the ray before and after scattering:

p(θ, g) = 1 4π 1 -g 2 1 + g 2 -2g cos θ 3 2
(2)

The parameter g express the anisotropy of the phase function. For g = 0, the rays are scattered equally in all directions. For 0 < g < 1, the rays are scattered mainly in the same direction as the incoming ray (forward scattering). For -1 < g < 0, the rays are scattered mainly in the reverse direction (backward scattering).

The Henyey-Greenstein phase function was originally developed for radiations in the galaxy [START_REF] Henyey | Diffuse radiation in the galaxy[END_REF]. It can be used to approximate any phase function φ using the parameter g computed as:

g = Ω cos θφ(θ) dω (3)
We will use this important property in our work. • The albedo α = σ s σ s +σ a , express the relative importance of scattering and absorption. Materials with a high albedo (α ≈ 1) have mostly scattering effects, with little absorption. Conversely, in materials with a low albedo (α ≈ 0), absorption effects dominate. The albedo is a dimensionless quantity, between 0 and 1. As pointed out before (e.g. [DLR + 09]), the mean-free path gives the scale of the response for a given material. For two materials that differ only by , we can predict the response of the second material by scaling the response of the first material by the ratio of their mean-free-paths.

In the remainder of this paper, we parametrise using normalised coordinates. All spatial coordinates are divided by to give a dimensionless, generic parametrisation. The behaviour of translucent materials only depends on three parameters: α, g and η.

Since we're interested only in the outside appearance of translucent materials, we model their behaviour using a bi-directional sub-surface scattering reflectance distribution function S (BSSRDF). It defines the general transport of light between two points and directions as the ratio of outgoing radiance at point

x o in direction ω o , L o (x o , ω o ) to the incident radiant flux Φ i (x i , ω i ) at point x i from direction ω i : S (x i , ω i , x o , ω o ) = dL o (x o , ω o ) dΦ i (x i , ω i ) (4)

The experimental framework

Our entire experimental setting follows closely Donner et al. [DLR + 09], for easier comparison with their results. We restrict ourselves to spatially uniform, homogeneous materials, whose boundary is an infinite plane. This reduces the number of dimensions as we only need the relative positions in space and angle. We express the spatial difference x o -x i in polar coordinates: r and θ s . We orient our axes so that the angular dependence of the BSS-RDF depends only on three angles (θ i , θ o , φ o ). For a given material, the BSSRDF depends on five parameters: S (r, θ s , θ i , θ o , φ o ). We compute the BSSRDF using Monte-Carlo simulation. We send photons along a collimated beam, hitting the translucent material with a direction ω i making angle θ i with the surface normal. These photons are refracted as they enter the material, and propagate a distance d before being scattered or absorbed:

d = -log ξ σ s + σ a ( 5
)
Where ξ is a random variable in [0, 1]. This expression for d corresponds to the exponential falloff for illumination. When a photon is scattered, we compute its new direction by importance-sampling the phase function [START_REF] Wann | Realistic image synthesis using photon mapping[END_REF] and modulate its intensity by the albedo. If a photon hits the interface, we compute its reflected direction and modify its power by applying the Fresnel term.

To increase the speed of our simulations, we select sample points on the surface of the material. For each sample point x o , for each photon in the simulation, we compute the power Φ e that this photon would contribute if it was scattered directly to this sample point and average the contributions of all the photons:

Φ e = α p( d • ω p , g) e -d F t ( d)Φ p (6)
• Φ p is the power of the photon, • d is the vector joining the position of the photon and the sample point,

• d = d/ d is the normalised version,
• ω p is the direction of propagation of the photon, • F t is the Fresnel term for transmission. If the angle between d and the normal n is larger than sin -1 1 η , then the photon is fully reflected at the material interface and does not contribute to this specific sample point. We sample the BSSRDF for the incoming directions: θ i ∈ {0, 15, 30, 45, 60, 70, 80, 88} at sample points defined by: r ∈ {0.01, 0.05, 0.1, 0.2, 0.4, 0.6, 0.8, 1, 2, 4, 8, 10}. Directions θ s are sampled regularly every 15 degrees. We store outgoing directions θ o and ϕ o in buckets of width 5 degrees. We used the same sampling strategy as [DLR + 09], for easier comparison and validation.

The whole simulation uses roughly 350 lines of code, and we provide it as supplemental material for researchers wishing to reproduce our results2 . We checked that we reproduce the results of [DLR + 09], even though we use a different method for sampling outgoing directions: we used regular samples in θ o and ϕ o while they used directions uniformly distributed on the sphere using the Healpix format.

Correlation between two-scatter events and multiple scatter events

In this paper, we only focus on multiple scattering events. Single scattering events are more complex and exhibit subtle spatial variations [START_REF] Walter | Single scattering in refractive media with triangle mesh boundaries[END_REF].

Among multiply scattered photons, we separate between photons scattered exactly twice and the others. We compute the outgoing light leaving the material after exactly Figure 2: Relationship between radiance leaving the material after two scatter events (x) and multiple scatter events (y). The two quantities are strongly correlated, with r = 0.998

(θ i = 0, α = 0.95, g = 0.3, η = 1.3).
two scatter events, along with the light leaving the material after multiple scatter events. Figure 2 shows the relationship between these two quantities, measured at multiple sample points and directions on the surface. Looking at the curve, you can identify two areas: Where there is a strong double scattering (L > 0.02), there is also a strong correlation between double and multiple scattering. When double scattering is negligible (L < 0.02), the relationship is not linear, although double scattering is still a good predictor of single scattering.

Figure 2 corresponds to a single material: η = 1.3, α = 0.95, g = 0.3, and a single input direction, θ i = 0. We conducted an extensive study for a complete range of materials, and found that this strong correlation always exists, except for materials with a strong backward scattering (g < -0.5). Note that Figure 2 corresponds to a material with a very high albedo, giving a strong importance to multiple scattering.

In our study, we tried all possible combinations of material parameters and input directions. For each material, we compute the correlation coefficient between outgoing light due to multiple scatter events and the light due to only two scatter events. For two distributions X and Y, the correlation coefficient r XY is defined as:

r XY = Cov(X, Y) σ X σ Y
r XY is equal to 1 for perfectly correlated distributions, and 0 if the distributions are independent from each other. We found that the relationship holds for all materials with g > -0.5, for all input directions. Thus light leaving after two scatter events is a strong predictor of multiply scattered light, except for materials with strong backward scattering properties. (c) Index of refraction η (α = 0.99, θ i = 0) Figure 3: Evolution of the correlation coefficient r between second-order scatter events and multiple scatter events, depending on the simulation parameters. Observe that we consistently observe a strong correlation for materials with g > -0.5. increasing the importance of two-event scatterings among multiple scattering, and increasing the correlation coefficient. In all figures except this one, we used a large value for the albedo (α = 0.99): the correlations we find can only get better with lower values of the albedo. For strongly backward scattering materials (g < -0.5), there is no correlation between second-order scattering and multiple scattering. Incoming direction: Figure 3(b) shows the variation of the correlation coefficient r as a function of g for different values of the input direction θ i . When there is a strong correlation between second-scatter and multiple-scatter, we see that this correlation remains the same, independently of the input direction: for g > -0.5, the curves follow similar patterns, even though there are slight variations. For strongly backward scattering materials, we see more variations in the correlation coefficient (there is more correlation at grazing angles), but as we are looking for a model that is independent from input direction this information is not very useful to us. More importantly, the correlation we find is also independent on θ i : for g = 0, α = 0.99, η = 1 (illustrated on Figure 2) we consistently get y ≈ 1.04x + 0.02, for all θ i . This key property means that we can replace the computation of multiple scattering by a computation of two-event scattering, and then deduce multiple scattering by interpolating. Note that in this case (high albedo, isotropic scattering), two-event scattering amounts for 96 % of multiple scattering. Index of refraction: Figure 3(c) shows the variation of the correlation coefficient r as a function of g for different values of the index of refraction η. As you can see, η has little influence on the quality of the correlation. The main difference is between η = 1 (no refraction at the interface) and η > 1. Materials with a index of refraction larger than 1 can experience total internal reflection: light reaching the surface above a certain angle does not exit and is reflected towards the material. This effect reduces both two-events scattering and multiple scattering, but apparently increases the correlation between the two. In all figures except this one, we placed ourselves in the worst case scenario: η = 1. Thus changing the index of refraction can only increase the correlations we found.

Conclusion and discussion

We have found a strong correlation between light exiting the material after two scattering events and multiply scattered light, for a large range of materials (forward scattering and moderately backward scattering, g > -0.5). We found no visible correlation for strongly backward scattering materials (g < -0.5). The strong correlation holds for all values of the albedo, incoming direction and index of refraction. We can use this correlation to speed-up the computation of multiple scattering in translucent materials: we compute two-event scattering and extrapolate to extract multiple scattering. The correlation coefficients depend on the actual material parameters (g, α and η) but can be precomputed. Figure 7 illustrates the results of this extrapolation. Where double scattering is low, we cannot use the affine correlation y = ax + b, as it would overestimate multiple scattering for low values of double scattering. We extended the affine correlation in a C 1 way so that its value is null at the origin (see Figure 2):

f (x, a, b, p) =      ax + b if x > x 0 b 1-p x x 0 p if x ≤ x 0 where x 0 = pb a(1-p)
We obtain the a, b and p parameters by fitting the computed values. The fact that multiply-scattered light is strongly correlated with two-event scattering explains a key re- Figure 4: We compute the density of two scatter events around the scattered ray inside the material, along with average direction. At rendering time, we find the contribution at sample point P using the density, average direction and phase function. sult of Donner et al. [DLR + 09]: they found that for most materials, lobes were aligned with the scattering plane (defined by the entry point and direction and the exit point). This is obviously true for light scattered twice. Since multiply-scattered light is correlated to that, we naturally get the same result.

We expect this result can be used with any light simulation algorithm such as path tracing, Monte-Carlo ray-tracing or Bi-Directional Path Tracing to speed up the computation of sub-surface scattering. In the next section, we provide a specific algorithm for low-memory computation of two-event scattering.

4

Fast computation of two-scatter events

In this section, we present an algorithm for fast computation of double-scattering events in translucent materials. When the light ray enters the material, it changes direction with refraction, but keeps propagating in a straight line. All the first scattering events take place on this line. As a consequence the distribution of second-scattering events has rotational symmetry around this line. We exploit this property to store photon density in a compact way.

We store the photon density in the medium, using the refracted incoming ray as the main axis for our frame. Since the photon density has rotational symmetry around this axis, we simply have to store a 2D function of ρ and z.

In the rendering step, for each outgoing ray, we compute the refracted ray inside the material and extract the contribution from the second scatter event on this ray (see Figure 4).

Precomputation: density and direction

Figure an incoming ray hitting the material, at an angle θ i with the normal. This ray is refracted when it hits the surface. The new ray makes an angle θ i with the surface normal, with θ i defined by Snell's law: η sin θ i = sin θ i . We express the incoming ray after refraction as R.

All first scatter events take place on R. As a consequence, the probability distribution of second scatter events has rotational symmetry around R (if we allow second scatter events to take place above the surface of the material). We introduce cylindrical coordinates based on R: (ρ, z) (see Figure 4). We compute three quantities:

• the probability density of second scatter events, d,

• the average direction for all photons scattered twice reaching that point v, stored as an angle in 2D, θ. • the anisotropy of the directions of photons reaching this point, expressed as the g parameter of a Henyey-Greenstein phase function, and computed using equation 3. We store all second scatter events, even those happening outside of the medium, and store them before multiplication by the albedo. Thus the quantities we compute and store only depend on the anisotropy of the phase function g, not on the incoming direction θ i , the index of refraction, or the albedo.

We do our precomputation in 2D-space for better efficiency. The density in 3D is connected to the density we computed in 2D by:

d 3D (ρ, z) = 1 2πρ d 2D (ρ, z) (7)
After precomputation, we store these functions into 2D textures. Using a texture resolution of 128 × 256, the cost for a three-channels, 4 byte-per-channel texture is 512 KB. The overall cost for storing 20 textures sampling g regularly from -0.9 to 0.95 is 7.5 MB.

To compute the main direction v and anisotropy, at each sample point, we compute the outgoing radiance in all directions, taking into account the phase function and exponential falloff at both scatter events. We sum the contributions for all primary scatter events, resulting in an anisotropic lobe (see Figure 5). We approximate this anisotropic lobe with a Henyey-Greenstein phase function, using equation 3, and store its main direction and value.

Sampling at runtime

At runtime, we have a sample point x o and an outgoing direction ω o and we need the outgoing light at this point. We takes ω o the refracted version of the outgoing ray (the ray that is such that η sin θ o = sin θ o ). We call S the refracted outgoing ray: (x o , ω o ).

If we consider a single point P along S , we obtain the contribution of double scattering on this point P by:

• express P in the frame of reference of the incoming ray R and divide coordinates by the mean-free-path. • extract the density d P , main direction v P and anisotropy g P at point P by bi-linear interpolation from the pre-computed values, • compute the phase function p v P • ω o , g P , • and finally multiplying by the density, albedo, Fresnel term and density volume element at this sampling point, dV:

L P = σ 2 s d P p v P • ω o , g P F t (ω o ) dV (8)
Depending on our application, we either want to use the contribution of a single point P (for example to connect light and eye paths in bi-directional path-tracing) or integrate the contributions of all points P along the ray S . For the latter, we sample at points P along the ray using the exponential falloff as importance function. We sample regularly for ξ ∈ [0, 1[ and take P = x olog(1 -ξ)ω o . This reduces the number of sample points while taking into account the exponential falloff with the distance.

Direct BSSRDF lobes comparison

Figure 6 shows a comparison between the lobes for two scatter events generated by our approximation and Monte-Carlo simulation, for different values of the material parameters. Notice how our approximation fits all the features of the lobes computed using Monte-Carlo approximation, both for size and shape. We are even able to capture the secondary lobe which happens for slanted incoming lighting (see second row).

Comparison with full Monte-Carlo simulation in a simple case

Figure 7 shows a side-by-side comparison between our algorithm and a reference, computed using Monte-Carlo simulation. We chose a simple setting so that all illumination effects would be visible: a collimated beam of light hitting a half-space with a planar surface, filled with a blue translucent material: σ s = (0.07, 0.53, 0.52), σ a = (0.93, 0.47, 0.48), g = 0. We used a refraction index of η = 1.0 to visualize only scattering effects.

The two pictures on the left display only single and double scattering. You can see that our algorithm provides a good approximation of double scattering. The main problem is that we underestimated the light propagated along the direction of the beam (visible on the right edge of the picture). We attribute this to under-sampling for small values of r. Even for isotropic material, most of the interesting events take place close to the original beam. Instead of regular sampling, we could use more samples for small values of r.

The two pictures on the right show all scattering events. We computed multiple scattering by extrapolation from double scattering in our algorithm, and compare it with the reference solution using Monte-Carlo. Our algorithm exhibits the right behaviour in all directions. We tend to underestimate illumination in the direction of the beam (this is a consequence of underestimating double scattering), and to overestimate it in areas at the periphery of the halo.

Integration in a ray-tracer

We ported our algorithm in the Mitsuba renderer [START_REF] Wenzel | Mitsuba renderer[END_REF], as an extension of its subsurface class, which uses Jensen and Buhler [JB02] method. Our implementation is an extension of [START_REF] Wann | A rapid hierarchical rendering technique for translucent materials[END_REF]: first, we precompute incoming light at randomly distributed sample points and directions on the surface of the object. In a second pass, we integrate the contributions from these directional samples. All rendered images in this paper were generated using this implementation. The main difference with [START_REF] Wann | A rapid hierarchical rendering technique for translucent materials[END_REF] is that we need to store the Figure 6: Comparison between two-scattering lobes computed using our approximation (left) and Monte-Carlo (1.2 million photons). Both columns are using the same scaling. Notice how we pick the right direction and shape for the lobes, including the presence of a secondary lobe on the second row.

Full model Reflections/Refractions

All Scatterings incoming radiance in both space and direction. This significantly increases the number of samples we need to store at runtime, depending on material density and anisotropy of the phase function.

Figure 1 displays a comparison between our method (center column), the classical dipole approximation and a reference solution. For the reference solution, we used bidirectional path-tracing and a large number of samples per pixel (4,000). All pictures were computed using the same material definition, inside the same renderer (Mitsuba). The computations were spread over three networked computers with Intel core 2 processors, for a total of 16 cores (two machines with 4 cores, one machine with 8 cores). We used 32 samples per pixel for the dipole approximation, for more accurate shadows and anti-aliasing, and 16 for our implementation.

The top row shows multiple scattering effects only: the result of the dipole approximation (left), extrapolating from double scattering for our method (center) and multiple scattering effects for the reference solution. To present only multiple scattering effects for the reference solution, we used the following method: first, we computed a full reference solution, with all effects included. Then, we subtracted from it a picture with only reflection, refraction and single scattering, computed by restricting the number of events in bi-directional path tracing.

The bottom row shows the picture with all lighting effects: reflection, refraction, translucency and single scattering combined with multiple scattering (for bi-directional path-tracing, it corresponds to the picture computed without any tweaks).

The dipole approximation is very fast (16 seconds), but differs in terms of color. Our implementation, using precomputed double scattering effects and extrapolating multiple scattering, provides the right color, while being 20 times faster than the reference solution, and only 6 times slower than the dipole approximation Figure 8 displays the relative contributions of all the elements used in rendering translucent materials: specular reflections on the surface, transmittance through the material, both of these combined, then multiple scattering (extrapolated from double scattering) and the full model. The effects of scattering are clearly visible on the central part of the statue.

Figure 9 shows the effects of varying material density, while keeping other material parameters constant. This is equivalent to changing the scale of the object compared to the mean-free-path. We show separately the effects of reflection and refraction (left column), multiple scattering extrapolated from double scattering (center column) and all of them combined together (right column). As the object gets larger (top row) or smaller (bottom row), scattering effects are less visible, compared to reflection and refraction effects. Scattering effects are more visible if the object has the right size compared to the meanfree-path, not too large nor too small.

Conclusion

Our first contribution is a study of sub-surface scattering effects in translucent materials. Except for materials with strong backward scattering (g < -0.5), our experiments show a strong correlation between multiple scattering and double scattering. Any light simulation algorithm can exploit this property to speed up the computations of multiple scattering in sub-surface scattering, by computing only the double scattering and extrapolating. What we have found is a an affine correlation: it does not mean that multiple scattering does not exist after the second event, it means that, on average, its distribution can be predicted by the results of the second event. Our study was only conducted on surface events; we haven't studied the distribution of light inside translucent materials. This would require a separate study.

Our second contribution is an approximate model to compute double-scattering effects in an efficient way. We store the average material response, independently from the incoming direction. We use these precomputed values to get the response from the second scatter event, and extrapolate to get multiple scattering. We integrated this model inside a sub-surface scattering plug-in, but this model could be combined with any illumination simulation algorithm. Our model is an approximation, but provides accurate colors and fair light distribution, while being very effective in terms of memory and CPU.

Figure 1 :

 1 Figure 1: Comparison between our algorithm and standard algorithms for homogenenous scattering material 1 . Top row: multiple scattering effects only. Bottom row: full light simulation (reflection, refraction single and multiple scattering). The dipole approximation (left) results in the wrong color. Reference solution (right), computed using bi-directional path tracing (BDPT), takes a long time. Multiple scattering effects extrapolated from double scattering computed using our algorithm (center) is 20 times faster than the reference solution and provides the right color.

Figure 3

 3 presents our results. It displays the variation of the correlation coefficient as a function of the input parameters of the simulation. A coefficient r > 0.95 corresponds to a very strong correlation; r > 0.8 corresponds to a strong correlation. Albedo: Figure 3(a) shows the variation of the correlation coefficient r as a function of g for different values of the albedo α. At each scatter event, the light is multiplied by α. Reducing the albedo reduces the amount of light being multiply scattered, Albedo α (η = 1, θ i = 0) Input angle θ i (α = 0.99, η = 1)

  Cylindrical coordinates (ρ, z), centered on the refracted ray

Figure 5 :

 5 Figure5: We approximate the contributions at each sample point with a single Henyey-Greenstein lobe, with direction v and anisotropy g.

  (a) Reference, single and double scattering (b) Our algorithm, single and double scattering (c) Reference, single and multiple scattering (d) Our algorithm, single and multiple scattering extrapolated from double scattering

Figure 7 :

 7 Figure 7: A collimated beam, making an angle of π/3 with the surface normal, hitting the surface of a half-space containing a blue translucent material. Comparison between values computed using our algorithm and full Monte-Carlo simulation. Material: σ t = 1, α = (0.07, 0.53, 0.52), η = 1, g = 0.

Figure 8 :

 8 Figure 8: Compositing the different components of the model. From left to right: all effects except scattering, all scattering effects (single, double and multiple), and the full model (reflections, refraction and scattering). Material: σ t = (0.44, 0.88, 0.55), α = 0.6, η = 1.3, g = 0.5.

Figure 9 :

 9 Figure 9: The effects of varying material density. Left column: only reflected and transmitted light. Central column: scattered light (single scattering plus multiple scattering, extrapolated from double scattering). Right column: complete lighting. Material: σ s = (0.2, 0.3, 0.6), α = 0.9, η = 1.3, g = 0.2.

See the IEEE digital library for this paper.

Acknowledgments

The Armadillo, Asian dragon and Thai statuette models appearing in the figures of this paper are kindly provided by the Stanford Computer Graphics Laboratory. All pictures were ray-traced using the Mitsuba renderer [Jak10], with our custom plug-in for double scattering effects.


Keywords: deb, starvation, stochastic feeding, rejuvenation, respiration, fish

Developmental patterns of yolk-sac larvae are well captured by the standard deb model: (i) when feeding is delayed post birth the size at which post-feeding growth begins is reduced but the rate of growth post-feeding is unaffected and (ii) maternal effects (initial energy in egg) show up as differences in condition at birth and maximum length of non fed individuals. We extended the standard deb model in two ways to account for starvation. (I): if somatic maintenance can no longer be paid structure is also mobilized to cover the costs, but at an extra cost-conversion efficiency of structure to energy. Death occurs if structure reaches a fraction of the maximum at the onset of shrinking. (II): if maturity maintenance can no longer be paid then maturity level decays exponentially (rejuvenation). Hazard due to rejuvenation is proportional to the difference between maturity and the maximum maturity at the onset of rejuvenation.

We performed Monte Carlo simulation studies which treat feeding as a random process to evaluate the contribution of the metabolic handling of starvation to early teleost life history. The simulations suggest that food density strongly impacts growth, energy reserves, mineral fluxes, hazard and mortality from shrinking. Environmental factors can soon override maternal induced differences between individuals. Moreover in the low food density, simulated individuals from eggs of lower caloric content experience mortality from shrinking earlier than their counterparts issued from higher energy eggs. Empirically observed patterns of real data, i.e. high scatter in respiration in combination with low scatter in lengths, can be expected when the metabolism is treated as a deterministic system while behaviourally controlled input is stochastic. At low food densities where mortality from shrinking reaches 10% almost all individuals experience hazard due to rejuvenation. This hazard is difficult to access experimentally but represents moments of heightened susceptibility to pathogens and toxicants and could be ecologically significant.

List of Figures

Modified from [START_REF] Kamler | Parent-egg-progeny relationships in teleost fishes: an energetics perspective[END_REF]. Egg size is measured in many different ways. Assumptions of auxiliary theory link these measurements to amounts of reserve and structure (egg size to initial amount of energy in reserve), which comes with some parameters (see table 2) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Impact of delays in first feeding on growth at reference temperature (20 • C). s H = 2; k′ J = 2.4; f = 1 (see table 2). (A) structural length (L) against age for E 0 = 1.7 (black lines) and E 0 = 1.1 (grey lines); arrowheads: birth; dashed lines: instantaneous death by shrinking. E 0 = 1.7: feeding is never initiated (negative growth till death by shrinking) and feeding is delayed 0, 4 and 8 days after birth. E 0 = 1.1: feeding feeding is never initiated and feeding is delayed 0 and 4 days after birth. Survival probability for E 0 = 1.67 (B) and E 0 = 1.1 (C). Survival for delays in initial feeding of 0 days after birth (*), 4 days after birth (**), 8 days after birth (***) and feeding is never initiated ( †). Dotted line: all organisms have died from shrinking. . . . . . . . . . . . . . . . . . . . . . . . . . 3

Points: data from [START_REF] Jardine | Direct yolk sac volume manipulation of zebrafish embryos and the relationship between offspring size and yolk sac volume[END_REF] (28 • C). (A) maximum observed length attained during development in the absence of feeding against yolk volume expressed in nanolitres (nl). Dashed line: model predictions of this maximum length for 1.1 ≤ E 0 ≤ 2.5 J. (B) observed length against age for eggs with all yolk volumes shown on left. Solid dark grey: model predictions supposing κ G = 0.67, solid light grey: model predictions supposing κ G = 1, dashed grey: death from shrinking in both cases, i.e. the length reached a fraction δ X = 0.75 of the maximum at the onset of negative growth. . . . . . . . . . . . . . . . . . 4

Model outputs computed at 20 • C. s H = 2 and k′ J = 0.024 day -1 . First row: an example of scaled reserve density (e), observed length (L f ), wet mass, rejuvenation hazard ( ḣJ ) for an individual sampled randomly from the Monte Carlo simulation study. Second row: an example of fluxes of CO 2 , H 2 O, O 2 and NH 3 against age since birth from the same individual sampled randomly from the Monte Carlo simulation study. Grey: stochastic mineral fluxes for that individual (x = 0.75, E 0 = 1.67 J, f alternates between 0 and 1) are compared to deterministic predictions (black smoth line) for f = x/(1 + x) and f = 0 (E 0 = 1.67 J). When f = 0 the deterministic mineral fluxes go to 0 and the dotted line is when structure reaches a fraction δ X = 0.75 of the maximum structure (instantaneous death is assumed at this point). Bottom row: distribution of each mineral flux at age 30 since birth for 5000 Monte Carlo trials (x = 0.75, E0 = 1.67 J, f alternates between 0 and 1). Please note that the x and y axis of fluxes of NH 3 differ. . . . . . . . . . . . . . . . . . . . . . . . . .

2

Output of 5000 Monte Carlo simulation studies at reference temperature (20 • C). s H = 2 and k′ J = 0.024 day -1 . Open symbols: x = 0.75, closed symbols x = 0.3, circles: E 0 = 1.8 J, squares: E 0 = 1.2 J.

List of Tables

Introduction

The ecology of teleost larvae has generated much interest in marine fisheries research. Survival through the early life history stages is thought to regulate recruitment and subsequent yearclass strength [START_REF] Hjort | Fluctuations in the great fisheries of Northern Europe viewed in the light of biological research[END_REF][START_REF] May | Larval mortality marine fishes and the critical period concept[END_REF][START_REF] Cushing | Marine ecology and fisheries[END_REF]Bailey and Houde, 1989, among others).

Starvation [START_REF] Hjort | Fluctuations in the great fisheries of Northern Europe viewed in the light of biological research[END_REF] and predation [START_REF] Bailey | Predation on the eggs of larvae of marine fishes and the recruitment problem[END_REF]) have been hypothesized as major factors controlling. Predation has been considered to be the primary agent of mortality during the yolk-sac stage [START_REF] Blaxter | The role of sensory systems of herring larvae in evading predatory fishes[END_REF][START_REF] Paradis | Vulnerability of fish eggs and larvae to predation: review of the influence of the relative size of prey and predator[END_REF] while starvation may only become important after the transition to exogenous feeding [START_REF] Leggett | Recruitment in marine fishes: is it regulated by starvation and predation in the egg and larval stages?[END_REF] as there are potential interactions of larval size [START_REF] Miller | Larval size and recruitment mechanisms in fishes: toward a conceptual framework[END_REF], developmental rate [START_REF] Houde | Early life dynamics and recruitment variability[END_REF], environmental effects [START_REF] Shepherd | The effect of development rate on swimming behaviour and escape response in American plaice (Hippoglosoides platessoides) larvae[END_REF], predation and starvation. To date, with much of the emphasis focused on how growth and survival relate to environmental conditions, few studies consider how larval metabolism handles starvation. The focus of our paper is (i) to formalise starvation rules which operate at the individual level and (ii) further investigate the importance of metabolism on starvation during teleost early life history stages. The study is conducted within the conceptual framework of Dynamic Energy Budget (deb) theory [START_REF] Sousa | From empirical patterns to theory: A formal metabolic theory of life[END_REF][START_REF] Kooijman | Dynamic Energy Budget Theory for metabolic organization[END_REF][START_REF] Sousa | Dynamic energy budget restores coherence in biology[END_REF]: a well-tested theory on the uptake and use of substrate by all organisms. Parameters for the standard deb model (see Kooijman, 2010, Chap.2) have been estimated for a number of animal species and can be found in the Add My Pet collection http://www.bio.vu.nl/thb/deb/deblab/add_my_ pet/ (Lika et al., 2011, this special issue). The standard deb model is extended in this paper to deal with starvation in more detail. The question we ask is how does the metabolic handling of starvation contribute to growth, survival, dioxygen consumption, and carbon dioxide, water and ammonia production (hereafter referred to as mineral fluxes) during the early life history of the individual when resources are fluctuating. It is practical to conduct this study using a deb model where the parameters of a teleost have already been estimated as well as with data which relate egg size to growth and survival. A full life cycle model for zebrafish, Danio rerio was previously developed [START_REF] Augustine | Developmental energetics of zebrafish, Danio rerio[END_REF] and [START_REF] Jardine | Direct yolk sac volume manipulation of zebrafish embryos and the relationship between offspring size and yolk sac volume[END_REF] have unique results where initial egg size was controlled by microinjection technique and lengths, yolk volumes and survival were observed for each individual. This makes zebrafish an ideal candidate for conducting theoretical studies on starvation.

Within the context of deb theory the life cycle is divided in three stages: embryo (no assimilation), juvenile (assimilation) and adult (allocation to reproduction but no longer to maturation). Two important events bound the embryonic stage: (i) age zero where reserve is maximal and starts to reduce and structure and maturity are close to zero and start to increase and (ii) age at birth where external feeding is initiated and the embryo switches to a juvenile. Hatch precedes birth for many teleost species since the mouth is not yet formed and the individual is still reabsorbing the yolk sac inherited from the mother. Birth is initiated the moment the oesophagus opens and the digestive tract is open to the environment. Birth corresponds to a specific developmental milestone and so it makes sense to compare ages and lengths at this point in development between different individuals to understand maternal effects. Age and length at hatch is a more commonly found endpoint, but has the drawback of preceding birth and not occurring at strictly the same stage of development for all individuals (e.g. Kimmel et al., 1995, for zebrafish). Furthermore, age at hatch can be influenced by rearing protocol: sterilisation procedures such as bleaching can delay hatching by hardening the chorion (e.g. [START_REF] Zhang | Depletion of zebrafish Tcap leads to muscular dystrophy via disrupting sarcomere-membrane interaction, not sarcomere assembly[END_REF]. Yolk consists of lipo-proteins, and during embryonic development it is rapidly converted into proteins, lipids and carbohydrates while hardly loosing energy or building blocks. These three products are stored in different places in the body, a process known as internalisation of yolk. These 'details' are not included in the standard model explicitly, and the yolk-sac is treated as a temporary organ with specialised functions, similar to adipose tissue of juveniles and adults in many taxa.

We first examine to what extent the deb model is representative of impacts of maternal effects (initial egg caloric content) and delays in initial feeding on early juvenile development. Predictions are compared to actual observations of length against age of individuals born of eggs with initial yolk volumes ranging from 126 to 628 nl (see Jardine and Litvak, 2003, for detailed protocol).

It has already been shown that stochastic food availability can explain the large variance in size of individuals in a same aquarium from a same brood due to amplification of difference by social interaction (Kooijman, 2009a). Large variation in larval zebrafish lengths (and high mortalities) were demonstrated to be a consequence of low food density, which was amplified by changes in food types, from small to larger food particles [START_REF] Eaton | Growth and the reduction of depensation of zebrafish, Brachydanio rerio, reared in the laboratory[END_REF]. To emulate a random encounter rate between food items and juveniles we consider a natural stochastic version of the feeding module and perform Monte Carlo simulation studies to generate lengths, mass, reserve, mineral fluxes and survival for 5000 individuals. We perform the studies at low and high food densities and with populations born of different initial egg energy content. The importance of including the metabolic handling of starvation when studying the ecology of juvenile teleosts is discussed. and juvenile) using 3 state variables: energy in reserve E (J), structural length L (cm) and maturity E H (quantified in J) (see table 1) and 7 parameters (see table 2). L = V 1/3 where V (cm 3 ) is the structural biovolume.

The assimilation flux ṗA (J day -1 ) is ṗA = 0 for embryos and ṗA = f { ṗAm }L 2 for juveniles with { ṗAm } the maximum (surface-area) specific assimilation rate, and f the scaled functional response, defined as the feeding rate on a particular food type as a fraction of the maximum possible one for an individual of that size. It can be specific for a particular diet. The standard model uses the Holling type II function response, f = x/(1 + x), with x the scaled food density.

A fraction κ of the flux ṗC (J d -1 ) that is mobilized from reserve is invested in somatic maintenance ṗM = [ ṗM ] V and growth ṗG = κ ṗC -ṗM . A fraction κ G (the conversion efficiency) of the flux allocated to growth is fixed in new structure. (Lika et al., 2011, this special issue), with [E G ] the cost per unit of structure and µ V , d V , w V the chemical potential, density and molar weight of structure (see table 2).

κ G = µ V d V ([E G ]w V ) -1
A fraction (1 -κ) ṗC is invested in maturity maintenance ṗJ = kJ E H and maturation ṗR = (1 -κ) ṗC -ṗJ . Maturity does not have mass or energy itself, but is quantified as the cumulative energy investment in maturation. Maturity maintenance can be conceived as all processes responsible for maintaining the current state of maturity, i.e. immune, hormonal and cellular defence systems. An analogy would be that considerable energy is invested in learning and if the knowledge level is not maintained it is forgotten.

The mobilization flux is given by ṗC = E ( v/L -ṙ), with v (cm d -1 ) the energy conductance and ṙ (d -1 ) the specific growth rate. Reserve dynamics are specified as follows: d dt E = ṗA -ṗC . Maintenance always has priority over investment in growth or maturation. Starvation occurs when energy mobilised no longer suffices to cover (I) the somatic maintenance cost, i.e. κ ṗC < [ ṗM ] V , with shrinking as result and (II) the maturity maintenance costs, i.e. (1-κ) ṗC < kJ E H , with rejuvenation as result. It seems reasonable to assume that shrinking has a maximum and death occurs instantaneously when it exceeds a fixed fraction δ X . This construct comes with the need to introduce a new state variable max L, i.e. the maximum length the individual once had. It seems likely that maturity maintenance is more optional than somatic maintenance in the sense that paying for the immune or the cellular defence systems might not be obligatory. Yet there is a penalty in the form of an increased hazard proportional to the fraction of maturity maintenance that is not paid (assuming implicitly that there is a constant need for defence). Also this construct comes with a need to introduce a new state variable max E H , i.e. the maximum maturity level the individual once had.

Both shrinking and rejuvenation have an extra parameter in the specification of their dynamics. On the assumption that the somatic maintenance costs expressed as energy flux is the same for growth and shrinking, the change in structural length is [ ṗM ] and the energy investment ratio g = [E G ] v κ{ ṗAm } are compound parameters of the standard deb model. Somatic maintenance might have a 'building block' aspect that might cause growth efficiency κ G to deviate from the value mentioned above.

Following the previous analogy we can imagine that forgetting follows a first order process with rate parameter k′ J , say:

dE H dt = (1 -κ) ṗC -kJ E H if positive, else dE H dt = -k′ J E H - (1 -κ) ṗC kJ
Hazard due to rejuvenation ḣJ is considered proportional to the difference between actual maturity level and maximum maturity level at the onset of rejuvenation max

E H : ḣJ = kM max E H -E H E h H and E h H = (1 -κ) [E G ] L 3 m κ s H
where s H is the rejuvenation stress coefficient. The result is that ḣJ contributes to survival as:

dS dt = -S( ḣJ + ḣa + ḣacc + ḣsh )
In this study we consider that hazard from ageing ḣa and from accidents ḣacc is negligible. We consider that hazard from shrinking of structure ḣsh is infinite with instantaneous death as a result when the shrinking threshold is passed.

The starvation module has, in summary, two new state variables, max L and max E H , and four new parameters κ G , δ X , k′ J and s H . Growth efficiency also occurs under non-starvation conditions, but there it does not play a role in the dynamics of the state variables.

Embryo and juvenile mineral fluxes are specified in the Appendix. For a full discussion we refer to Kooijman (2010, Chap.4).

Linking state variables to measurements

Structural length L is taken proportional to some well-chosen physical (observed) length L f : L = δ M L f , with δ M the (constant) shape correction factor. [START_REF] Augustine | Developmental energetics of zebrafish, Danio rerio[END_REF] show that the shape of early juvenile zebrafish changes during the first month of development. This is most likely the case for many teleost species who undergo metamorphosis after early juvenile development. We exclude this detail from our study by considering that δ M is constant.

Freshly spawned eggs are taken to be composed uniquely of reserve inherited from the mother. Maternal effects are expressed as differences in the initial energy content of the egg E 0 (J). We relate observed egg volume to initial energy content as follows:

V Y f = E 0 δ 3 Y w E µ E d E
δ Y is an empirical shape correction factor. w E , d E and µ E are the molar weight, density and chemical potential of reserve (table 2). We assume an Arrhenius relationship between metabolic rates and temperature (Eqn. 1.2 [START_REF] Kooijman | Dynamic Energy Budget Theory for metabolic organization[END_REF]. Computations are made at reference temperature of 20 • C. Predictions which are compared to observations by [START_REF] Jardine | Direct yolk sac volume manipulation of zebrafish embryos and the relationship between offspring size and yolk sac volume[END_REF] are temperature corrected to 28 • C. Dry biomass W d has contributions from both reserve and structure and is expressed as:

W d = d V L 3 + w E µ E E
with d V the density of structure (table 2). We make the simplified assumption that wet biomass W w is equal to 6 W d on the basis of observations by [START_REF] Craig | Growth and mortality of zebra fish, Brachydanio rerio (Hamilton Buchanan), maintained at two temperatures and on two diets[END_REF]; [START_REF] Bagatto | Growth and metabolism of larval zebrafish: effects of swim training[END_REF].

Stochastic feeding module

Suppose that food density, in terms of particles per volume, is constant, such as in a polyculture (e.g. [START_REF] Best | A novel method for rearing first-feeding larval zebrafish: polyculture with Type L saltwater rotifers (Brachionus plicatilis)[END_REF]. We simplify the system by considering all food items to be of constant chemical composition with mass M X (µmol). Each fish is either searching for food (i.e. f = 0) or handling food (i.e. f = 1). Time interval spent handling food t h is inversely proportional to squared length (of the individual), so larger individuals take less time to handle a food item. { JXAm } is the surface-area specific maximum intake rate (in C-mol food per day per surface area) and relates to surface-area maximum assimilation rate { ṗAm } as: { JXAm } = { ṗAm }y XE /µ E . y XE is the yield of food on reserve (see table 2).

t h = M X { JXAm } L 2
Food searching intervals t s are considered as an exponentially distributed random variable with mean t h /x. Mean feeding rate is defined as ḣX = (t s + t h ) -1 = f /t h . This natural introduction of stochasticity does not come with any new parameters.

Two populations are created from two pools of initial energy contents: E 1 0 = 1.18 J and E 2 0 = 1.68 J. Populations E 1 0 and E 2 0 both develop in two different environments: x 1 = 0.3 and x 2 = 0.75. This makes 4 experimental conditions: (i) E 1 0 , x 1 , (ii) E 1 0 , x 2 , (iii) E 2 0 , x 1 , and (iv) E 2 0 , x 2 . We run 5000 Monte Carlo simulation for each condition and represent results as the distributions of lengths, mass, scaled reserve density e, and respiration at 30 days since birth. In addition, we look at the distribution of age at death from shrinking for the entire study and the number of live individuals at 30 days since birth who have experienced rejuvenation hazard. Stochastic simulation routines (traject.m and traject m.m) are part of the DEBtool software [START_REF] Kooijman | From food-dependent statistics to metabolic parameters, a practical guide to the use of dynamic energy budget theory[END_REF] and are freely downloadable at http://www.bio.vu.nl/thb/ research/.

Results and Discussion

Constant food

Model simulations show that when initial feeding is delayed post birth, size at which postfeeding growth begins is reduced but the rate of growth post-feeding is unaffected (figure 2A). This is in line with general observations by [START_REF] Guillaume | Nutrition et alimentation des poissons et crustacés[END_REF]. Individuals from the lower E 0 die from shrinking if feeding is delayed beyond 4 days. The initial caloric content of the egg influences post-birth lifespan when feeding is never initiated which is again in agreement with general patterns noted by [START_REF] Kamler | Parent-egg-progeny relationships in teleost fishes: an energetics perspective[END_REF]. Rejuvenation hazard contributes to reducing survival probability (figure 2B andC). The time an individual can handle delays in initial feeding before rejuvenation hazard starts contributing to survival relates to initial energy in the egg.

The model predicts that maternal effects (initial energy in egg) influence condition at birth in terms of amount of reserve per unit structure and maximum length of non fed individuals. Birth does not occur if E 0 < 1.1 J because the cumulated energy invested in maturity never attains the threshold value E b H (see table 2). Hence eggs are considered non viable when E 0 < 1.1 J. It would be interesting to verify experimentally the realism of this lower boundary value of E 0 for which an egg is viable. Should such a lower boundary exist then maternal effects can also contribute to embryonic survival.

The maternal effect rule (Kooijman, 2009b), implemented by default into the standard deb model, stipulates that the reserve density of the mother equals the reserve density of her offspring at birth. E 0 = 1.67 J for mothers at f = 1 (for this parameter combination). Between the narrow caloric range of E 0 = 1.1 to 1.67 J, the maximum structure for a non fed juvenile increases by 10%, and the maximum maturity attained increases by 40%. This explains the increased hazard experienced by lower E 0 individuals during starvation (figure 2C). Caloric values of egg dry mass for teleost species are generally very conserved ranging from 20 -30 J mg -1 [START_REF] Kamler | Parent-egg-progeny relationships in teleost fishes: an energetics perspective[END_REF]. If we consider the dry mass of an egg as W 0 = E 0 w E / µ E , then the caloric value of E 0 = 1.6 J is 22.9 J mg -1 which falls within the boundaries suggested by [START_REF] Kamler | Parent-egg-progeny relationships in teleost fishes: an energetics perspective[END_REF].

We do not propose a maximum egg size in this study. Links between egg size and the mother's condition might be more complex and depending on the applications deviations from ing on growth at reference temperature (20 • C). s H = 2; k′ J = 2.4; f = 1 (see table 2). (A) structural length (L) against age for E 0 = 1.7 (black lines) and E 0 = 1.1 (grey lines); arrowheads: birth; dashed lines: instantaneous death by shrinking. E 0 = 1.7: feeding is never initiated (negative growth till death by shrinking) and feeding is delayed 0, 4 and 8 days after birth. E 0 = 1.1: feeding feeding is never initiated and feeding is delayed 0 and 4 days after birth. Survival probability for E 0 = 1.67 (B) and E 0 = 1.1 (C). Survival for delays in initial feeding of 0 days after birth (*), 4 days after birth (**), 8 days after birth (***) and feeding is never initiated ( †). Dotted line: all organisms have died from shrinking.

the maternal effect rule should be considered. The maximum surface-area specific assimilation rate { ṗAm } is sensitive to the type of food. Consequently, maximum assimilation rate which corresponds to f = 1 for one type of food might correspond to f < 1 for another. Since the maternal effect rule implies that the e of the embryo at birth equals the value of f of the mother at spawning and that diets differ between individuals, using E 0 > 1.67 J as starting values to calculate life history traits is coherent with model assumptions.

In figure 3A we computed lengths of non fed juvenile zebrafish for 1.1 ≤ E 0 ≤ 2.5 J as suggested by the large range of observed yolk volumes V Y f . Predicted maximum length as a function of V Y f is compared with the actual observations from [START_REF] Jardine | Direct yolk sac volume manipulation of zebrafish embryos and the relationship between offspring size and yolk sac volume[END_REF]. Maximum size attained seems to saturate for initial yolk volumes above 314 nl. This points to an uneasy relationship between yolk volume and water content. Large eggs tend to have higher caloric value than small eggs, but this is not always the case [START_REF] Kamler | Parent-egg-progeny relationships in teleost fishes: an energetics perspective[END_REF]. For Sardina pilchardus, [START_REF] Riveiro | Importance of parental effects on larval survival in Sardina pilchardus[END_REF] found that that time to yolk absorption was not a function of egg size but a function of egg biochemical composition (protein content) which illustrates the point we make in figure 1. We converted E 0 to V Y f with the implicit assumption that the fraction of water per volume is constant. There is no solid evidence that this is the case. Still, the data are overall supportive of a positive relationship between E 0 and maximum length of a non fed individual.

We further compare predicted lengths of individuals issued from eggs with energy content 1.1 ≤ E 0 ≤ 2.5 J to observed lengths against age for each V Y f (figure 3B). We show two cases: (i) paying maintenance from structure increases somatic maintenance costs by a factor κ -1 G and (ii) there is no extra cost involved (κ G = 1) which is not thermodynamically coherent. In both cases predicted negative growth rate at the onset of maximum structure is overestimated. The discrepancy might be explained by parameters values. On the other hand, does an organism shrink isomorphically? [START_REF] Yin | Morphological changes during growth and starvation of larval cod (Gadus morhua L.) and flounder (Platichthys flesus L.)[END_REF] found that it was the ratios of morphometric measurements which were sensitive to the state of starvation for larval Gadus morhua L. and Platichthys flesus L.. A length measurement following the spinal cord, as is the case here, might be relatively insensitive to shrinking. Nonetheless, the model reflects the observed relationships between V Y f and maximum lengths while lengths at starvation are underestimated. The predicted effects of initial energy in egg on growth and metabolic handling of starvation lends support to current views where maternal effects contribute to resistance to starvation. Yet it remains difficult to assess the importance of this maternal effect in real life situations. We tackle the question theoretically using Monte Carlo simulations with stochastic food input in the following section. during development in the absence of feeding against yolk volume expressed in nanolitres (nl). Dashed line: model predictions of this maximum length for 1.1 ≤ E 0 ≤ 2.5 J. (B) observed length against age (hours post fertilization hpf) for eggs with all yolk volumes shown on left. Solid dark grey: model predictions supposing κ G = 0.67, solid light grey: model predictions supposing κ G = 1, dashed grey: death from shrinking in both cases, i.e. the length reached a fraction δ X = 0.75 of the maximum at the onset of negative growth.

Stochastic searching

All stochasticity that we propose is in food searching, not in food handling. This is natural because we assume a fixed particle size for simplicity's reasons, consistent with the idea that the standard deb model pushes simplicity into the extreme. In many natural situations, food particles will show a particle-size distribution, which introduces stochasticity in handling as well, but also comes with new parameters and the results become dependent on the particle size distribution, which itself will be affected by many factors.

The values of e, L f , W w, ḣJ and mineral fluxes against age since birth of an individual experiencing a stochastic encounter with food particles in environment E 1 0 , x 1 is compared to deterministic predictions for an individual at f = x 1 /(1 + x 1 ) (expected mean) and a starved individual (f = 0). See figure 4, first and second row.

Stochastic mineral fluxes show rapid small scale variation since f is fluctuating continuously between 0 and 1. The distribution of daily mineral production (CO 2 , H 2 O, NH 3 figure 4I, J and L) or consumption (O 2 figure 4K) of 5000 individuals at the end of the simulation time shows considerable deviation from the mean. Future research however could include digestion (see e.g. Kooijman, 2010, Chap.7) which would make assimilation a smoother process at small time scales. Differences in prior food history will still accumulate however and contribute to inter-individual differences in lengths, mass, reserve and mineral production and consumption. The mineral flux curves all show the same morphology, but this is a result of choice of parameter values and not a model property. The mineral fluxes are closely tied to energetics and go to zero when the individual is starved. No individuals in the x = 0.75, E 1 0 environment experience rejuvenation hazard. Figure 4D shows rejuvenation hazard when f = 0 and the organism is shrinking.

The Monte Carlo simulation studies reveal the predicted scatter in lengths, weights, survival and respiration at different time points when food availability fluctuates. Monte Carlo results of simulations for each of the 500 individuals at day 30 since birth are presented in figure 5. Interestingly the largest scatter is observed in the predicted mass and respiration data and the lowest in the predicted length data: the coefficient of variation (CV) of mass and O 2 consumption is almost four time times the CV of length in all cases. In the x 2 , E 1 0 environment for example the CV is 9% for L f , 34 % for W w and 38% for O 2 consumption. Monte Carlo simulations of deb models with stochastic food input seem to inherently capture patterns of scatter observed in real data sets where high scatter in measured mass or respiration combines with much lower scatter in observed lengths.

The CV increases twofold between the high and the low food density environments for all predicted observations regardless of the initial energy content of the egg. For example, the CV of the L f distribution for x 1 , E 1 0 is 4% wheareas it is 9% for x 2 , E 1 0 (figure 5A). This links up beautifully with the results of [START_REF] Eaton | Growth and the reduction of depensation of zebrafish, Brachydanio rerio, reared in the laboratory[END_REF] where lengths of individuals sampled after 16 days of development in containers with 2 free-swimming prey per mL had a CV of 75, E 0 = 1.67 J, f alternates between 0 and 1) are compared to deterministic predictions (smooth lines) for f = x/(1+x) and f = 0 (E 0 = 1.67 J). When f = 0 the deterministic mineral fluxes go to 0 and the dotted line is when structure reaches a fraction δ X = 0.75 of the maximum structure (instantaneous death is assumed at this point). Bottom row: distribution of each mineral flux at age 30 since birth for 5000 Monte Carlo trials (x = 0.75, E0 = 1.67 J, f alternates between 0 and 1). Please note that the x and y axis of fluxes of NH 3 differ. 16.8% and lengths of individuals sampled after 15 days of development in containers holding 200-300 free-swimming prey per mL had a CV of 3.8%.

Populations E 1 0 and E 2 0 evolve in a very similar manner in terms of lengths, weights, respiration and survival at both food densities (x = 0.3 and x = 0.75). Thus the present simulations hint at a negligible contribution of maternal effects to measurable endpoints after a month of development regardless of the food density. Maternal effects can be diluted by the individuals (chance) encounter rate with food items.

The most sensitive difference between E 1 0 and E 2 0 at x = 0.3 is in the distribution of age at death (figure 5E). Population E 2 0 shows higher mortalities from age 8 to 20 since birth than the population E 1 0 . After 20 days, age at death distributions are quasi identical between both populations. At x = 0.3 about 10% of both populations die from structure shrinking. At x = 0.2 about 50% of both populations die from shrinking (not shown).

Almost all individuals in both populations in the low food density environment (x = 0.3) experience hazard due to rejuvenation. This is shown in figure 5F where 79% of live individuals experienced hazard in E 1 0 , x = 0.3 and 85% of live individuals experienced hazard in E 2 0 , x = 0.3. Rejuvenation hazard hardly occurred for individuals in x 1 , E 1 0 and x 2 , E 1 0 . This hazard ( ḣJ ) can be translated biologically as moments of increased susceptibility to pathogens and stressors and contributes to survival probability. It is difficult to evaluate ḣJ experimentally, but the present simulations suggest that it could be a relevant factor regulating survival should additional environmental stressors and small periods of starvation occur simultaneously.

Concluding remarks

There exists strong empirical support for maternal effects as a general pattern but our simulations show how environmental conditions soon override such effects. Inter-individual differences in prior food history accumulate and can be responsible for the observed scatter in the data.

The problem of inherent variability in biological data must be addressed in all fields of biology. The present study shows that empirically observed patterns of real data, i.e. high scatter in respiration in combination with low scatter in lengths, can be expected when the metabolism is treated as a deterministic system while behaviourally controlled input is stochastic. This is consistent with general patterns in observations.

Monte Carlo simulation studies are powerful tools to understand how stochasticity of food searching translates into variability in measured endpoints such as respiration, growth or reproductive output. This insight can be useful for designing experiments: e.g. feeding protocols, types and frequencies of measurements and sampling.

The concept of rejuvenation and rejuvenation hazard follows directly from the structure of the DEB model but further theoretical work needs to be done to understand how to estimate it from data. The concept links up neatly with empirical support that starved organism show a heightened susceptibility to disease and environmental stressors such as toxicants [START_REF] Rougier | Copper and zinc exposure to zebrafish, Brachydanio rerio (hamilton-buchaman): effects in experimental listeria infection[END_REF].

The extension we propose of the standard DEB model specifying shrinking and rejuvenation is not species-specific and applies to the entire life-cycle. Further dedicated research will test to what degree these rules are inherent to the functioning of all organisms and further our understanding of which adaptations are species specific and/ or specific to a particular life stage of the organism.

  (A) Distribution of lengths of live individuals after 30 days . (B) Distribution of wet mass of live individuals after 30 days. (C) Distribution of scaled reserve density (e) of live individuals after 30 days. (D) Distribution of respiration (nmol O 2 day -1 ) of live individuals 30 days. (E) Age at death by shrinking for each individual. Black: E 1 0 , grey: E 2 0 . Note that no deaths occurred for x = 0.75. (F) Percentage of live individuals who experience hazard due to rejuvenation during the simulation interval. Black: E 1 0 , grey: E 2 0 . No hazard is experienced for individuals in E 1 0 , x = 0.75. Percentage of population in E 2 0 , x = 0.75 is barely perceptible on the graph. . . . . . . . .

  d dt L = L ṙ/3 with specific growth rate ṙ = kM g e/l -1 e + g if positive, else ṙ = kM g e/l -1 e + κ G g where l = L/L m is scaled length and e = E v L 3 { ṗAm } is scaled reserve density. The somatic maintenance rate coefficient kM = [ ṗM ]/[E G ], the ultimate structural length L m = κ { ṗAm }

Figure 2 :

 2 Figure 2: Impact of delays in first feed-

Figure 3 :

 3 Figure 3: Points: data from Jardine and Litvak (2003) (28 • C). (A) maximum observed length attained

Figure 4 :

 4 Figure4: Model outputs computed at 20 • C. s H = 2 and k′ J = 0.024 day -1 . First row: an example of scaled reserve density (e), observed length (L f ), wet mass, rejuvenation hazard ( ḣJ ) for an individual sampled randomly from the Monte Carlo simulation study. Second row: an example of fluxes of CO 2 , H 2 O, O 2 and NH 3 against age since birth from the same individual sampled randomly from the Monte Carlo simulation study. Black: stochastic mineral fluxes for that individual (x = 0.75, E 0 = 1.67 J, f alternates between 0 and 1) are compared to deterministic predictions (smooth lines) for f = x/(1+x) and f = 0 (E 0 = 1.67 J). When f = 0 the deterministic mineral fluxes go to 0 and the dotted line is when structure reaches a fraction δ X = 0.75 of the maximum structure (instantaneous death is assumed at this point). Bottom row: distribution of each mineral flux at age 30 since birth for 5000 Monte Carlo trials (x = 0.75, E0 = 1.67 J, f alternates between 0 and 1). Please note that the x and y axis of fluxes of NH 3 differ.

Figure 5 :

 5 Figure 5: Output of 5000 Monte Carlo simulation studies at reference temperature (20 • C). s H = 2 and k′ J = 0.024 day -1 . Open symbols: x = 0.75, closed symbols x = 0.3, circles: E 0 = 1.8 J, squares: E 0 = 1.2 J. (A) Distribution of lengths of live individuals after 30 days . (B) Distribution of wet mass of live individuals after 30 days. (C) Distribution of scaled reserve density (e) of live individuals after 30 days. (D) Distribution of respiration (nmol O 2 day -1 ) for live individuals at day 30. (E) Age at death by shrinking for each individual. Black: E 1 0 , grey: E 2 0 . Note that no deaths occurred for x = 0.75. (F) Percentage of live individuals who experience hazard due to rejuvenation during the simulation interval. Black: E 1 0 , grey: E 2 0 . No hazard is experienced for individuals in E 1 0 , x = 0.75. Percentage of population in E 2 0 , x = 0.75 is barely perceptible on the graph.

Table 1 :

 1 State variables of the deb model.

	Variable Unit Name
	Standard deb model
	E	J	reserve
	L	cm	structural length
	E H	J	cumulated energy invested in maturity
	Forcing variable
	x	-	scaled food density
	Additional state variables
	max L	cm	maximum structural length before shrinking
	max E H	J	maximum E H before rejuvenation

Table 2 :

 2 Model parameters (affecting state variables at 20 • C) .

	Parameter Value	Unit	Name
	deb model parameters (values taken from Augustine et al., 2011)
	{ ṗAm }	246.3	J day -1 cm -2	maximum surface-area assimilation
	v	0.0287	cm day -1	energy conductance
	[ ṗM ]	500.9	J day -1 cm -3	volume-linked somatic maintenance costs
	[E G ] kJ	4652 0.0166	J cm -3 day-1	cost of a unit of structure maturity maintenance rate
	κ	0.44	-	allocation fraction to soma
	E b H	0.54	J	cumulated energy invested in maturity at birth
	Temperature correction module (values taken from Augustine et al., 2011)
	T A	3000	K	Arrhenius temperature
	Auxiliary parameters (values taken from Lika et al., 2011, this special issue)
	w E	23.9	g C-mol -1	molar weight of reserve
	w V	23.9	g C-mol -1	molar weight of structure
	d E	0.15	g cm -3	density of reserve
	d V	0.15	g cm -3	density of structure
	µ E µ V	555 500	kJ C-mol -1 kJ C-mol -1	chemical potential of reserve chemical potential of structure
	Parameters linked to starvation (this study)
	k′ J	depends on study day -1	specific maturity decay
	s H	depends on study -	rejuvenation stress coefficient
	δ X	0.75	-	maximum shrinking fraction
	Parameters linked to respiration(this study)
	y XE	1.25	C-mol X/ C-mol E yield of food on reserve
	y P X	0.1	C-mol P/ C-mol X yield of faeces on food
	Shape coefficients used for comparing model output to data from Jardine and Litvak (2003)
	δ Y	1.2521	-	yolk shape coefficient
	δ M	0.1245	-	shape coefficient
	Stochastic feeding module
	M X	0.06	µmol	Mass of food particle

Information Number of DEB auxiliary theory assumptions

  Modified from[START_REF] Kamler | Parent-egg-progeny relationships in teleost fishes: an energetics perspective[END_REF]. Egg size is measured in many different ways. Assumptions of auxiliary theory link these measurements to amounts of reserve and structure (egg size to initial amount of energy in reserve), which comes with some parameters (see table 2, auxiliary parameters).

	Diameter Volume	Wet mass	Dry mass	Energy	Content of a key substance
	Figure 1:

Model and methods2.1 Formulation of starvation rulesWe use embryo parameter values for zebrafish provided in[START_REF] Augustine | Developmental energetics of zebrafish, Danio rerio[END_REF] while simplifying their model by neglecting deviations from isomorphy and effects of ageing. We also ignore any surface-area linked maintenance costs. As mentioned before, we focus on the early juvenile that follow the standard deb model[START_REF] Kooijman | Dynamic Energy Budget Theory for metabolic organization[END_REF] and have two life stages (embryo

Acknowledgements

This work is part of the ENVIRHOM research program supported by the Institute for Radioprotection and Nuclear Safety and the Provence Alpes Côte d'Azur region. Béatrice Gagnaire and Christelle Adam-Guillermin are gratefully acknowledged for the stimulating discussions and useful comments. We would like to thank two anonymous referees for their constructive comments on the manuscript

Appendix

Computations of fluxes of dioxygen ( JO ), carbon dioxide ( JC ), water ( JH ) and ammonia ( JN ) comes with the specification of the stoichiometry of elements ( matrix n M of chemical indices) and the ratios of C, H, O, N to C in the generalized molecule which makes up the four organic compartments: X (food) , V (structure), E (reserve), P (faeces) which are collected in matrix n O .

Chemical indices n O and n M (see Lika et al., 2011, this special issue) are as follows:

n O =

1.00 1.00 1.00 1.00 1.80 1.80 1.80 1.80 0.50 0.50 0.50 0.50 0.10 0.10 0.10 0.10

Organic fluxes are a weighted sum of the three basic powers in deb theory: assimilation ṗA , dissipation ṗD and growth ṗG . Assimilation and growth powers are defined in the main text (section 2.1). Somatic maintenance, maturity maintenance and maturation all contribute to the dissipation power such that ṗD = [ ṗM ] L 3 + (1 -κ) ṗC . Mineral fluxes are a weighted sum of the organic fluxes. Faeces production is linked only to assimilation so is zero for embryos.

The four organic fluxes, JX flux of food (C-mol X day -1 ), JV structure flux (C-mol V day -1 ), JE reserve flux (C-mol E day -1 ), and JP feces flux (C-mol P day -1 ) are specified as follows :

and f = 0 for embryos. The four mineral fluxes JC (mol C day -1 ), JH (mol H day -1 ), JO (mol O day -1 ), and JN (mol N day -1 ) are a weighted sum of organic fluxes (see Eqn.4.37 Kooijman, 2010, pp.139