Nadine Le Bolay

Abdoulaye Seyni

Mehrdji Hemati

WATER UPTAKE AND DIFFUSION IN CO-GROUND STARCH-FILLED POLY (VINYL ACETATE) COMPOSITES

Keywords: Co-grinding, composite materials, degradable filler, water uptake, diffusion

The article presents results on the incorporation of starch in a poly (vinyl acetate) matrix by co-grinding. This treatment permits a good dispersion of the filler in the matrix and favors agglomeration phenomena which limit segregation and enhance interactions between the two materials. The composite materials offer interesting water resistance properties. Water absorption phenomena were studied and it was shown that water penetration in the materials is controlled by diffusion. The diffusion coefficients of water were determined for the matrix, the filler and the mixture and they are respectively equal to 3.3 10 -11 , 1.1 10 -10 and 5.5 10 -11 m 2 /s.

INTRODUCTION

Starch-filled composite materials have received considerable interest over the recent years due to the biodegradability, renewability and low cost of the filler. Among the different processes used to produce polymeric composite materials, co-grinding is the most recently developed. It allows a good homogeneity of the composite and good interactions between the two materials [START_REF] Zapata-Massot | Effect of the mineral filler on the surface properties of co-ground polymeric composites[END_REF], permitting to offer interesting properties. In a previous study [START_REF] Seyni | On the interest of using degradable fillers in co-ground composite materials[END_REF], it was shown that co-grinding makes it possible to enhance the mechanical properties of starch-filled poly (vinyl acetate) composites compared to blending.

Water uptake is generally the first step of biodegradation. It is strongly linked to the polymer hydrophilic behavior [START_REF] Verdu | Action de l'eau sur les plastiques[END_REF]. This one is influenced by the chemical structure. Indeed, it is low for carbonaceous or fluorinated groups, medium for moderately polar groups (like esters or ethers) and high for groups which are hydrogen bonds donors (acid, alcohol or amide). When a hydrophilic polymer is immersed in water, the liquid molecules diffuse in the polymer, leading to a material swelling occurring more or less rapidly, depending on the polymeric chains relaxation. Above the glass transition temperature, the chains relaxation is rapid and water diffusion can occur into the polymer network, following the Fick's laws. Below the glass transition temperature, chains relaxation can be slow and one can observe deviations from Fickian behavior. [START_REF] Alfrey | Diffusion in glassy polymers[END_REF] proposed three models to describe the transport phenomena into glassy polymers:

-Diffusion is Fickian if the water diffusion rate is lower than the relaxation rate. The mechanism is controlled by diffusion. -Water diffusion rate is higher than relaxation rate. The mechanism is controlled by relaxation. -The two rates are comparable. [START_REF] Friesch | Sorption and Transport in Glassy Polymers -A Review[END_REF] has expressed the water diffusion mechanism according to the following equation:

n i i t t k m m m m . = - - ∞ (1)
where m i is the initial mass of the polymeric sample, m t , its mass at immersion time t, m ∞ , its mass at equilibrium, n, the diffusional exponent and k, the system constant. This equation is considered to be valid for (m t -m i )/(m ∞ -m i ) < 0,6. The value of n depends on the diffusion mechanism. It is equal to 0.5 for Fickian diffusion, to 1 for a relaxation mechanism, and between 0.5 and 1 for the third case described above. Simplifying the second Fick's law, Crank (1975) has expressed equation (1) as follows:

t D e m m m m i i t         = - - ∞ π
In this paper, we present results on the influence of co-grinding poly (vinyl acetate) and starch on the behavior of the mixtures in water. Co-grinding was realised in a laboratory tumbling ball mill. We will first comment the size variation to explain how the composite materials are produced. Then, we will present the evolution of the water uptake of co-ground samples. From these results, we will study the influence of co-grinding on the water diffusion mechanism and coefficient.

II. APPARATUS, METHODS AND MATERIALS

The products were ground in a 5-liter ceramic tumbling ball mill containing ceramic balls (5.5, 9.5 and 12.5 mm). The balls occupied 20 vol. % of the chamber and the powder 10 vol. % of the void space between the balls. The rotating speed of the chamber was fixed at 75 % of the critical speed, i.e. 100 rpm.

The matrix was poly (vinyl acetate) (PVAc). The initial particles were sieved between 80 and 100 µm. As it can be observed in Figure 1.a, the particles, produced by emulsion polymerization followed by drying, are spheroids. They are composed of small grains. The PVAc density is 1783 kg.m -3 . The filler was waxy maize starch (Figure 1.b) with a particles size between 4 and 30 µm. Individual particles and agglomerates compose the initial powder. The starch density is 1330 kg.m -3 . The PVAc -starch mixtures contained 25 wt. % starch. Indeed, it has been shown in a previous study that this filler rate permits to improve significantly the mechanical properties of the composite [START_REF] Zapata | Production of small composite particles by cogrinding in a media mill -characterization of the granulometric and the mechanical properties[END_REF]. a: Poly (vinyl acetate) particles b : starch particles Fig. 1. SEM micrographs of poly (vinyl acetate) and starch particles Small powder samples were taken at various grinding times in different parts of the mill and analysed by means of a laser diffraction granulometer Malvern Mastersizer 2000 to determine their size distribution. The average size, d50, corresponding to a cumulated volume fraction of 50 %, was then calculated. Some particles were analysed by scanning electron microscopy (Leo 435 VP). Powder discs, with a mass around 0.6 g ± 10 -4 g, were formed by pressing the particles in a Carver press. The mould diameter was 8 mm and the applied force was 30 kN. The discs were introduced in distilled water contained in a flask and weighted at different times after being removed from water and rapidly dried on a blotting paper. Three discs were used for each series of experiments. The water uptake, WU, was calculated from the equation:

100 * i i t m m m WU - = (3)

III. RESULTS AND DISCUSSION

Size and morphology evolution

The analysis of the size and morphology evolution permits to understand the behavior of the products during grinding and to study the fragmentation and agglomeration phenomena occurring during the treatment. Results have been presented in details in a previous article [START_REF] Seyni | On the interest of using degradable fillers in co-ground composite materials[END_REF], and only the main results on this point are reported here, since they are useful to understand water uptake phenomena. PVAc and starch particles were ground separately first and together then. Figure 2 presents the variation of the mean sizes during grinding.

The mean size of PVAc particles decreases first due to a fragmentation of the particles, and after a limit size of 58 µm is reached, a competition between fragmentation and agglomeration occurs. The mean size of starch particles evolves differently, since it increases from the treatment beginning due to agglomeration, to oscillate then around 30 µm. These results are in agreement with those presented by [START_REF] Schönert | Size reduction (Fundamentals) -Chap. 1[END_REF] and may be explained by the fact that the limit size that can be reached when grinding polymer particles, which are macromolecules, is several tens of µm. Consequently, since the initial mean size of PVAc is about 86 µm, a little fragmentation may occur before agglomeration phenomena appear. On the contrary, initial starch particles are small (13 µm approximately, i.e. lower than the limit grinding size). So, the mean size increases due to immediate agglomeration. As for the co-ground mixtures containing 25 wt. % of starch, their size distributions are gathered in Figure 3. Before grinding, two peaks are observable, the right-hand side one corresponding to PVAc and the other to starch. The peak of PVAc shifts towards the left during the first 300 minutes, due to particle fragmentation already observed when grinding this polymer alone. During this period, the starch peak does not evolve significantly, but at 900 minutes, one can observe that it tends to disappear. Moreover, a peak appears around 350 µm. Finally, for longer co-grinding times, the main peak shifts to the right, while the proportion of the peak at 350 µm increases. Indeed, after PVAc particles are fragmented, small starch particles stick on the matrix fragments (left micrograph of Figure 4). Since filler particles are no more free in the mixture, they can not be detected by the granulometer. Furthermore, agglomeration continues as the processing time is prolonged, and some large agglomerates made up of the two materials are formed (right micrograph of Figure 4): PVAc-starch agglomerates tend to agglomerate between them. Thus, the filler is well dispersed in the matrix and no segregation may occur during handling.

The consequence of this is that the mean size of the mixture drawn in Figure 2 decreases strongly up to 300 minutes resulting from PVAc fragmentation, increases slightly up to 900 minutes because of an agglomeration of starch on PVAc, and increases significantly then due to an agglomeration of PVAc -starch agglomerates.

Finally a compaction of these big agglomerates, combined with a competition between fragmentation and agglomeration lead to oscillations of the mean size.

Behavior of the materials in water

PVAc discs produced by pressing particles ground different times (between 0 and 4200 minutes) have been immersed in water. Figure 5 shows how the water uptake evolves with the immersion time. The water uptake increases with the immersion time for values of this one lower than 1140 minutes (19 hours).

The affinity of the polymer with water may be attributed to the presence of acetate groups. Water is absorbed by the polymer. In this range of the immersion time, an increase of the grinding time leads to a decrease of WU, in particular for high grinding times. This may be the result of an amorphisation of the polymer or a degradation of the molecular chains. Between 1140 and 100 000 minutes (70 days), WU does not evolve significantly and equilibrium in the water uptake is reached. The influence of the grinding time on the values of WU at equilibrium is similar to that on the absorption rate in the first part of the curves. Above 70 days of immersion, a swelling of the discs and a degradation of their surface is observed. Consequently, WU increases again.

Starch discs were submitted to the same treatment. For grinding times lower than 60 minutes, the discs disaggregate rapidly. This is due to the good affinity between starch and water, resulting from the presence of all the -OH groups of the glucose molecules constituting starch. For grinding times between 120 and 300 minutes, similar curves as those with PVAc were obtained for immersion times lower than 1140 minutes, probably because of a modification of the molecular chains. However, when the maximum water uptake is reached, WU does not level off as for PVAc, but it decreases, certainly because of a partial starch release in water. For grinding times higher than 600 minutes, one notes a degradation of the starch discs, all the more fast as starch was ground a long time. This behavior can result from a too important degradation of the molecular chains under the action of the balls.

Finally, the same experiments were realized with discs of co-ground mixtures (Figure 6). The water uptake increases as for the matrix ground alone in the first zone of the curves (immersion times < 1140 minutes) and levels off then when an equilibrium is reached. It can be seen that WU is slightly higher for the mixture than for the matrix, because of the strong hydrophilic behavior of the filler that permits to enhance water absorption. For co-grinding times lower than 600 minutes, the absorption rate increases with the treatment time, on the contrary of PVAc. This may be related to the size and SEM observations presented in the previous paragraph. Indeed, as co-grinding occurs, the dispersion of the filler in the matrix and its adhesion on PVAc fragments (limiting segregation) are favored. For co-grinding times between 600 and 1800 minutes, the water uptake is strongly reduced, as for the matrix alone, probably for the same reasons, but it does not evolve for longer co-grinding times, since an equilibrium in the composite formation and thus in the water absorption kinetics may be reached. It has to be noted that the water uptake equilibrium period is shorter for the mixture than for the matrix because of a starch release in water after 12 000 minutes immersion, and WU decreases. The starch release is low for long co-grinding time, because starch is protected by the matrix. Finally a disc swelling is observed after 70 days of immersion, in particular for the shorter co-grinding times.

Fig. 6. Influence of the immersion time and the grinding time on mixtures water uptake These results are interesting since they mean that, when co-grinding the mixtures during an adapted time, it is possible to partially replace the synthetic matrix by a degradable and cheap filler, and to have a good resistance of the composite to water.

Analysis of water diffusion in the materials.

As shown in Figures 5 and6, water uptake increases with the immersion time because of a diffusion of water in the materials. Equation (1) has been used to characterize the diffusion mechanism. In this way, we have plotted in Figure 7 the variation of ln [(m t -m i )/(m ∞ -m i )] versus ln (t immersion ) for the matrix and the mixture. The points plotted on the figures represent all the results obtained for all grinding and co-grinding times. Fig. 7. Variation of ln [(m t -m i )/(m ∞ -m i )] versus ln (t immersion ) for the matrix and the mixture These points are well distributed around a single line, whatever the processing time of the powder in the mill. The exponent n, defined by the slope, is equal to 0,57 for the polymer alone and 0,52 for the mixture. These two exponents are close to 0,5, indicating that the mechanism of penetration of water within the discs is controlled by diffusion. Thus, using equation (2), the diffusion coefficients of water in PVAc and in the composite have been determined. They are respectively equal to 3.3 10 -11 m 2 /s for PVAc and 5.5 10 -11 m 2 /s for the mixtures, and they are independent of the grinding or co-grinding time. The diffusion coefficient is higher for the mixture because of the presence of starch. It was not possible to realise the same work with the filler data for all the grinding times because of the disintegration of most of the discs when immersed in water. However, we have teated the data obtained for 300 minutes of grinding and we have determined a diffusion coefficient equal to 1.1 10 -10 m 2 /s, which is similar to the value obtained by [START_REF] Russo | A study of water diffusion into a high-amylose starch blend: the effect of moisture content and temperature[END_REF]. The diffusion coefficient of the mixture logically lies between those of the two materials.

Mixtures

In order to go further, we were interested in the variation of the discs mass at equilibrium (Figure 8). It is constant for PVAc and increases for the mixture during the first grinding times, and then it decreases to stabilize for long grinding times. To continue the analysis, we plotted in Figure 9 the variation of m ∞ with 1/d50. One can observe two phases, the first one for ascending values of 1/d50 (i.e. for the size reduction) and the second one for decreasing values of 1/d50 (i.e. for the agglomeration of the particles). In the first phase, the mass at equilibrium is constant in the case of PVAc, which could be expected since this parameter depends on the material which is only the polymer. On the other hand, the mass increases slightly in the case of the mixture. This can be due to a modification of the disc structure, resulting of a modification, with co-grinding, of the dispersion of the filler in the matrix. In the second phase, the mass at equilibrium decreases in both cases, due to an amorphisation of the polymers or a chain degradation.

IV. CONCLUSION

A study on co-grinding of poly (vinyl acetate) and starch was implemented in order to show the possibility to replace partially the synthetic polymeric matrix by a vegetable, degradable and largely available filler. A size and morphology study has permitted to explain how the composite particles are formed. Then, a study of the behavior of the materials in water has led to the conclusion that introducing starch in PVAc increases the water absorption and the starch release, but these phenomena can be reduced by applying an adapted co-grinding treatment that may avoid the use of a compatibiliser to favor interactions between the matrix and the filler. It was shown that the water uptake is controlled by diffusion and the diffusion coefficients were calculated. That of the mixture is higher than that of the matrix because of the presence of starch which is strongly hydrophilic.

Fig. 2 .

 2 Fig. 2. Variation of the mean size of the two products and the mixture with the grinding and co-grinding time

Fig. 3 .

 3 Fig. 3. Evolution of the size distributions of the co-ground mixtures The disappearance of the starch peak and the presence of the peak at 350 µm may be explained by agglomeration phenomena that are illustrated by the SEM micrographs of Figure 4.

Fig. 5 .

 5 Fig. 5. Influence of the immersion time and the grinding time on PVAc water uptake

Fig. 8 .

 8 Fig. 8. Variation of m ∞ with the grinding time Fig. 9. Variation of m ∞ with 1/d50

(2) where e is the polymer sample thickness and D the diffusion coefficient. Several authors[START_REF] Peppas | Anomalous transport of penetrants glassy polymers[END_REF][START_REF] Russo | A study of water diffusion into a high-amylose starch blend: the effect of moisture content and temperature[END_REF], Yu et al., 2008[START_REF] Tang | Water absorption and solubility of PHBHV/HA nanocomposites[END_REF] studied water diffusion in hydrophilic polymers, using this equation to determine the diffusion coefficient in their materials.


Rosana El Jurdi

Olivier Colliot

HOW PRECISE ARE PERFORMANCE ESTIMATES FOR TYPICAL MEDICAL IMAGE SEGMENTATION TASKS?

Keywords: Segmentation, Performance, Validation, Statistical analysis, Confidence interval, Standard error

An important issue in medical image processing is to be able to estimate not only the performances of algorithms but also the precision of the estimation of these performances. Reporting precision typically amounts to reporting standard-error of the mean (SEM) or equivalently confidence intervals. However, this is rarely done in medical image segmentation studies. In this paper, we aim to estimate what is the typical confidence that can be expected in such studies. To that end, we first perform experiments for Dice metric estimation using a standard deep learning model (U-net) and a classical task from the Medical Segmentation Decathlon. We extensively study precision estimation using both Gaussian assumption and bootstrapping (which does not require any assumption on the distribution). We then perform simulations for other test set sizes and performance spreads. Overall, our work shows that small test sets lead to wide confidence intervals (e.g. ∼8 points of Dice for 20 samples with σ ≃ 10).

INTRODUCTION

In medical imaging, it is not uncommon that sample sizes are in the order of dozens of subjects, at best hundreds or thousands. In 3D medical image segmentation, the size of the set used to evaluate the performance may be even smaller than for other medical imaging tasks as obtaining the ground truth requires voxel-wise annotation by trained raters.

Intuitively, the precision of the estimation of the performance depends on two factors: the variability of the performance among the test set (the more variable, the less precise) and the size of the test set (smaller sets will lead to lower precision and therefore larger confidence intervals). However, papers usually report the average performance for different metrics (e.g. average Dice) but not the precision 1 with which this average performance is estimated. Such precision can be provided in the form of confidence intervals or equivalently standard error of the mean (SEM) which are not often reported. What is more often reported is the empirical standard deviation over different folds of a cross-validation. While this may qualitatively characterize the variability of the learning procedure when the training and testing set change, it should never be used to compute the SEM, since here n would be the number of folds or splits, which is arbitrary and can be made as large as one wants, thereby making the confidence interval arbitrarily narrow. It is not even an unbiased estimate of the standard deviation of the performance metric [START_REF] Bengio | No unbiased estimator of the variance of K-fold cross-validation[END_REF].

Quantifying the precision of the estimation of the performances thus requires an independent test set, on which confidence intervals or SEM are reported. Since this is not typically done in medical image segmentation papers, one may ask the following question. What precision can be expected for a typical sample size? How trustworthy are the average performance estimates (for instance Dice coefficients) reported in medical image segmentation papers? Surprisingly, this question has been little studied in medical imaging. In the case of a different task, namely image classification, it is necessary to have large sample sizes for a precise estimation of the accuracy (typically 10,000 samples to achieve a 1%-wide confidence interval given an accuracy of about 90% -95%) [START_REF] Varoquaux | Cross-validation failure: Small sample sizes lead to large error bars[END_REF][START_REF] Varoquaux | Evaluating machine learning models and their diagnostic value[END_REF]. However, to the best of our knowledge, this is not widely known in the case of segmentation. We hypothesize that the test size needed to achieve a given precision is lower than for classification due to the continuous nature of performance measures [START_REF] Douglas | The cost of dichotomising continuous variables[END_REF].

Our objective is to study the precision that can be expected in 3D medical image segmentation for typical test set sizes. We first conduct experiments using a standard deep learning network applied to a classical segmentation task from the Segmentation Decathlon Challenge [START_REF] Antonelli | The medical segmentation decathlon[END_REF] in order to estimate confidence intervals which are obtained for variable test set sizes. We then perform simulations for other sizes and spreads. We insist that the aim of the present paper is not to propose a new segmentation methodology. Instead, the main aims are to provide information regarding the confidence intervals that can typically be expected in medical image segmentation and to raise awareness of the community on this important issue.

OVERVIEW

Our aim is to provide confidence intervals (or standard-errors) for the mean of a given performance metric for different test set sizes and different spreads (standard deviation) of the performance. If the performance metric follows a Gaussian distribution, one can obtain those values as follows:

SEM = σ √ n CI = [µ -1.96 × SEM, µ + 1.96 × SEM] (1)
where µ denotes the mean, σ the standard-deviation, n the test set size, SEM the standard error of the mean and CI the 95% confidence interval. In the following, for conciseness, we will also denote the width of CI = [a, b] as w = b -a.

One does not know a priori if using Equations 1 is valid for a given performance metric. On the other hand, in the absence of assumption on the distribution of the metric, one can use bootstrapping to estimate SEM and CI [START_REF] Efron | An introduction to the bootstrap[END_REF]. In Section 3, we will perform experiments in order: i) to verify the validity of using Equations 1 by comparing the estimates obtained using these equations to those using bootstrapping; ii) obtain empirical estimates of µ and σ; iii) study the effect of using subsamples of reduced size.

In Section 4, we will perform simulations using Equations 1 to study how the precision varies when varying the test set size and performance variability (as measured by σ).

EXPERIMENTS

Dataset and segmentation method

We used the Hippocampus dataset from the Medical Decathlon challenge [START_REF] Antonelli | The medical segmentation decathlon[END_REF], composed of 260 3D MR images. The task is to segment the anterior and posterior parts of the hippocampus. For evaluation, we merged the two regions (in prediction and ground-truth respectively) and considered the hippocampus as a whole. From the 260 samples, 100 patients were randomly selected for training, 50 for validation and the remaining 110 samples constituted the test set.

For the segmentation, we used a U-net type network [START_REF] Ronneberger | U-Net: Convolutional networks for biomedical image segmentation[END_REF]. Note that our aim is not to achieve the highest possible segmentation scores or to propose a novelty in the segmentation method but rather to obtain typical performances. We thus relied on a standard approach. We treat a 3D MRI as a sequence of 2D images and used a 2D architecture. At inference, we predict for each slice independently then stack the slices belonging to the same patient together to form a 3D-volume prediction. The architecture is a 3-stage structure composed of convolutional and de-convolutional blocks, bottleneck and skip connections. An ensemble of convolutional and batch normalization layers constitutes the encoder part. Each stage within the decoder path is composed of 2 consecutive convolutional blocks followed by an upsampling layer. The bottleneck is composed of 2 convolutional blocks separated by a residual block [START_REF] Zhang | Road extraction by deep residual U-Net[END_REF]. The architecture has been used in previous publications [START_REF] Kervadec | Boundary loss for highly unbalanced segmentation[END_REF][START_REF] Jurdi | A surprisingly effective perimeter-based loss for medical image segmentation[END_REF]. The optimizer was Adam, the learning rate was 0.001 and the batch size was 8. The learning rate was halved if the validation performances did not improve over 20 epochs as proposed by [START_REF] Kervadec | Boundary loss for highly unbalanced segmentation[END_REF]. Note that we used the standard generalized Dice loss [START_REF] Sudre | Generalised dice overlap as a deep learning loss function for highly unbalanced segmentations[END_REF]. The model was trained over 500 epochs. We performed a three-fold cross-validation with the 150 patients of the training/validation sets and selected the best model across these three folds. Note that the test set was left untouched and was never used at any stage for training, model selection or architecture/parameter optimization. To evaluate the performance, we computed the Dice coefficient in %2 for the whole hippocampus (thus merging anterior and posterior parts before computing the metric).

The code used to generate the results is available online3

Precision on the whole test set

We first studied the precision of performance estimates using the maximum test set size. The distribution of Dice values over the test set is shown in Figure 1. One can observe that the Gaussian assumption is not unreasonable despite underlying outliers and skewness. More importantly, we will now compare estimates based on this assumptions to corresponding non-parametric bootstrap estimates. We first compute µ, σ, SEM and w using Equations 1. We then compute their bootstrap counterparts µ * , SEM * , w * as follows 4

SEM * = 1 M M m=1 (µ * m -µ * ) 2 . The 95% confidence interval CI * = [a * , b * ]
is the set of values between the 2.5% and 97.5% percentiles of the sorted bootstrap means {µ * 1 , µ * 2 , . . . , µ * m , . . . µ * M }. We finally define the width of the confidence interval as w * = b * -a * . The estimates using the Gaussian assumption and the bootstrap are very close as can be seen on Table 1.

Precision on subsamples of size k ≤ n = 110

We now study experimentally the relationship existing between the test set size and the precision of the estimation of the segmentation performance. To that end, we draw subsamples of variable size k ∈ K = {10, 20, 30, 50, 100, 110} from the whole test set of size n. In order not to depend on a particular drawing (which may be lucky or unlucky), we repeat the procedure 100 times for each k. We denote the subsamples as (S k,j ) where k is the subsample size and j ∈ {1, . . . , 100} is the index of a particular drawing.

We then proceed with the computations of the different estimates based either on the Gaussian assumption or on the bootstrap.

For the Gaussian assumption, we denote as µ k,j =

1 k k l=1 D j,k,l and σ k,j = 1 k k l=1 (D j,k,l -µ k,j ) 2 the
empirical mean and standard deviation for the subsample S k,j where D j,k,l is the Dice coefficient of a given subject in the subsample S k,j . Similarly, we use the notations SEM k,j = σ k,j √ k and w k,j = 2 * 1.96 * SEM k,j . We can then study how these values vary across the 100 subsamples of a given size k. To that end, we compute the average and standard-deviation of µ k,j , σ k,j , SEM k,j and w k,j across the different subsamples S k,j for k fixed and j ∈ {1, . . . , 100}. This provides the following estimates

µ k ± σ µ k , σ k ± σ σ k , SEM k ± σ SEM k and w k ± σ w k .
The values are displayed in Table 2. One can gather that, as the sample size increases, the standard deviation and the standard error decrease.

Bootstrap estimations are performed as follows. For a given subsample S k,j of size k and index j, M = 15000 bootstrap samples of size k are drawn with replacement. We denote a given bootstrap sample as S * k,j,m and its mean as µ * k,j,m where m ∈ {1, . . . , M } is the index of the m th bootstrap sample of subsample S k,j . The bootstrap mean µ * k,j of S k,j is the mean of the bootstrap sample means µ * k,j,m :

µ * k,j = 1 M m=1 µ * k,j,m .
The standard error of the mean µ * k,j (denoted as SEM * k,j ) obtained via bootstrapping is the standard deviation of the means of all bootstrap samples of subsample S k,j : SEM * k,j =

1 M M m=1 µ * k,j,m -µ * k,j 2 .
The 95% confidence interval of the sample S k,j is denoted as [a * k,j , b * k,j ] and is the set of values between the 2.5% and 97.5% percentile of the sorted bootstrap means of subsample S k,j . We finally define the width of the confidence interval via bootstrapping as w * k,j = b * k,j -a * k,j . We study how these values vary across the 100 samples of a given size k. To that end, we compute the averages and the standard deviations of µ * k,j , SEM * k,j and w * k,j across the different subsamples S k,j for k fixed and j ∈ {1, ..., 100}. This provides the following estimates µ

* k ± σ µ * k , SEM * k ± σ SEM * k and w * k ± σ w * k . Results are shown in Table 2.
As for the whole test set, estimates using Equation 1 and bootstrapping are very close across different subsample sizes. As expected, precision decreases with the sample size.

SIMULATIONS WITH A GAUSSIAN DISTRIBUTION

In the previous section, we showed that estimates of SEM and confidence intervals computed using either Equations 1 or the bootstrap are in accordance. Given this, we now perform simulations using Equations 1 for more values of k (including in particular larger test sets) and for different values of σ. Note that this is independent of µ which, in itself, has no impact on the SEM nor on the width of the CI (even though it is usually observed that lower performing models, thus associated with a lower value of µ, also have a more variable performance and thus a larger value of σ). Results are displayed in Table 3.

Comparing the gray column in Table 3 to Table 1, one can observe that both the standard error and the confidence interval width are close to the previously obtained experimental values (for k ≤ 100)). Nevertheless, the experimental values are slightly lower than those of the simulation. As the value of k increases, the gap between experimental and simulated values decreases.

DISCUSSION

In this paper, we have provided elements regarding the precision with which segmentation performance can be estimated in typical medical imaging studies. As hypothesized, the test set size needed to obtain a given confidence interval is smaller than in image classification. Typically, with 10,000 samples for classification (for a high performing model with over 90% accuracy), one obtains a 1.2%-wide CI. For our segmentation experiments, such width is obtain with only about 1000 samples. For a 4%-wide CI, one needs about 1000 samples for the aforementioned classification and about 100 samples for segmentation. Of course, one needs to keep in mind that this depends not only on the test sample size but also on σ of the performance. For example for σ = 5, the 1%-wide CI would then require between 300 and 500 samples.

Confidence intervals are rarely reported in medical im-age segmentation papers. To illustrate this, we have conducted a search for papers published in 2022 in IEEE Transactions on Medical Imaging (IEEE TMI) dealing with segmentation of 3D images. We found 51 papers containing 107 experiments (a given paper may include several experiments). Only 71 (66%) of these experiments (corresponding to 21 papers) were done on an independent test set, the rest reporting cross-validation results on training/validation sets. This is problematic for two reasons. First, the performance on cross-validation results can lead to optimistically biased performances [START_REF] Wen | Convolutional neural networks for classification of alzheimer's disease: Overview and reproducible evaluation[END_REF]. Moreover, they cannot be used to estimate standard errors. For the experiments that had an independent test set, the median size of the test set was 21 (minimum: 4, mean: 61.5, maximum: 697). In light of our experiments, it is not unreasonable that the associated confidence intervals will often be wide (typically between 9 and 10 for a test set with 20 samples). Finally, of the 21 papers that adopted an independent test set, only 3 of them (amounting to 6% of the 51 surveyed papers) reported confidence intervals or standard-error [START_REF] Zhao | Segmentation of the true lumen of aorta dissection via morphology-constrained stepwise deep mesh regression[END_REF][START_REF] Chen | Beyond mutual information: Generative adversarial network for domain adaptation using information bottleneck constraint[END_REF][START_REF] Han | Deep symmetric adaptation network for crossmodality medical image segmentation[END_REF]. Note that 11 other papers used statistical testing to compare different approaches even though they did not provide confidence intervals [START_REF] Song | Global and local feature reconstruction for medical image segmentation[END_REF][START_REF] Chen | Attention-assisted adversarial model for cerebrovascular segmentation in 3D TOF-MRA volumes[END_REF][START_REF] Jiang | Unpaired cross-modality educed distillation (CMEDL) for medical image segmentation[END_REF][START_REF] Ling | CAT-Net: A cross-slice attention transformer model for prostate zonal segmentation in MRI[END_REF][START_REF] Mehta | Propagating uncertainty across cascaded medical imaging tasks for improved deep learning inference[END_REF][START_REF] Huang | 3D lightweight network for simultaneous registration and segmentation of organs-at-risk in CT images of head and neck cancer[END_REF][START_REF] Chen | Semi-supervised unpaired medical image segmentation through taskaffinity consistency[END_REF][START_REF] Maria | KiU-net: Overcomplete convolutional architectures for biomedical image and volumetric segmentation[END_REF][START_REF] Zheng | Automatic liver tumor segmentation on dynamic contrast enhanced MRI using 4D information: Deep learning model based on 3D convolution and convolutional LSTM[END_REF][START_REF] Wang | RECISTSup: Weakly-supervised lesion volume segmentation using RECIST measurement[END_REF][START_REF] Chen | Model-Driven deep learning method for pancreatic cancer segmentation based on spiral-transformation[END_REF].

Our study has the following limitations. First, we studied only one dataset. Second, we used only one segmentation model. Future work would need to assess other models and datasets, in particular to see which is the typical range of values that can be expected for σ. Second, it was restricted to the Dice performance metric and it would be interesting to study if the same observations hold for other metrics (e.g. Hausdorff distance, volume error. . . ).

Overall, the experiments presented in our paper show the importance of reporting confidence intervals on independent test sets and that, in general, studies with small test sets cannot claim accurate estimation of the performances. We believe that it is an important issue on which the community should focus.

COMPLIANCE WITH ETHICAL STANDARDS

This research study was conducted retrospectively using human subject data made available in open access by the Medical Segmentation Decathlon challenge [START_REF] Antonelli | The medical segmentation decathlon[END_REF][START_REF] Amber L Simpson | A large annotated medical image dataset for the development and evaluation of segmentation algorithms[END_REF]. Ethical approval was not required as confirmed by the license attached with the open access data.

Fig. 1 :

 1 Fig. 1: Histogram of Dice accuracy over the entire test set. It is shown together with a kernel density estimation (KDE) which smoothes the observations with a Gaussian kernel.

  σ

Table 1 :

 1 . Given a test set of size n, M = 15000 Results on the full test set (n = 110). µ and σ are the empirical mean and standard deviation of the Dice coefficient across all patients in the test set. SEM is the standard error of the mean and w is the width of the 95% confidence interval calculated via Equations 1. SEM * , µ * and w * are the values obtained via Bootstrapping.

		n	µ	σ	SEM	w	µ *	SEM *	w *
		n = 110 80.70 10.75	1.02	4.02 80.70	1.02	3.99
	Subsample size k	µ k	σ k		SEM k		w k		µ * k	SEM * k	w * k
	k = 10	81.01 ± 3.04 8.17 ± 4.75	2.58 ± 1.5 10.13 ± 5.88 81.01 ± 3.04 2.59 ± 1.51 9.86 ± 5.62
	k= 20	80.61 ± 2.16 9.96 ± 3.74 2.23 ± 0.84 8.73 ± 3.28 80.61 ± 2.16 2.23 ± 0.84 8.61 ± 3.21
	k= 30	80.63 ± 1.6 10.36 ± 2.94 1.89 ± 0.54	7.41 ± 2.1	80.64 ± 1.6 1.89 ± 0.54 7.34 ± 2.08
	k= 50	80.95 ± 1.14	9.99 ± 2.1	1.41 ± 0.3		5.54 ± 1.16 80.95 ± 1.14 1.41 ± 0.3 5.51 ± 1.15
	k= 100	80.64 ± 0.32 10.78 ± 0.53 1.08 ± 0.05 4.22 ± 0.21 80.64 ± 0.31 1.08 ± 0.05 4.21 ± 0.21
	k= 110	80.70 ± 0.0	10.75 ± 0.0	1.02 ± 0.0		4.02 ± 0.0	80.71 ± 0.01 1.02 ± 0.01 3.99 ± 0.03

Table 2 :

 2 Results on subsamples of size k ≤ 110. Results are shown as mean ± std where mean and std are the mean and standard-deviation over all the subsamples S k,j of a given size k (k is fixed and j ∈ {1, . . . , 100}). bootstrap samples of size n are drawn with replacement. We denote a given bootstrap sample as S * m and its mean as µ * m where m ∈ {1, . . . , M }. The bootstrap mean µ * is the mean of the bootstrap sample means µ * m . The standard error of the mean µ * obtained via bootstrapping (SEM * ) is the standard deviation of the means of all bootstrap samples:

Table 3 :

 3 Simulation of SEM k and the width w k for different sizes k of the test set and for different values of σ. The gray column with σ = 10.75 corresponds to the standard deviation found in the experimental section.

		2		5		8		10.75	12		15		18
		SEM w k	SEM w k	SEM w k	SEM	w k	SEM	w k	SEM	w k	SEM	w k
	k = 10	0.63 2.48 1.58	6.2	2.53 9.92	3.4	13.33 3.79 14.88 4.74 18.59 5.69 22.31
	k = 20	0.45 1.75 1.12 4.38 1.79 7.01	2.4	9.43	2.68 10.52 3.35 13.15 4.02 15.78
	k = 30	0.37 1.43 0.91 3.58 1.46 5.73 1.96	7.7	2.19	8.59	2.74 10.74 3.29 12.88
	k = 50	0.28 1.11 0.71 2.77 1.13 4.43 1.52	5.96	1.7	6.65	2.12	8.32	2.55	9.98
	k = 100	0.2	0.78	0.5	1.96	0.8	3.14 1.08	4.22	1.2	4.7	1.5	5.88	1.8	7.06
	k = 200	0.14 0.55 0.35 1.39 0.57 2.22 0.76	2.98	0.85	3.33	1.06	4.16	1.27	4.99
	k = 300	0.12 0.45 0.29 1.13 0.46 1.81 0.62	2.43	0.69	2.72	0.87	3.39	1.04	4.07
	k = 500	0.09 0.35 0.22 0.88 0.36	1.4	0.48	1.89	0.54	2.1	0.67	2.63	0.8	3.16
	k = 1000 0.06 0.25 0.16 0.62 0.25 0.99 0.34	1.33	0.38	1.49	0.47	1.86	0.57	2.23

Throughout the paper, precision means how precise are the estimates of the performance. It has nothing to do with the performance metric Precision also known as Positive Predictive Value.

Performance metric was computed using this code: https:// github.com/deepmind/surface-distance

https://github.com/rosanajurdi/SegVal/tree/ ISBI2023

Throughout the paper, the bootstrap estimate of a given x is always denoted as x *

ACKNOWLEDGMENTS

The research leading to these results has received funding from the French government under management of Agence Nationale de la Recherche as part of the "Investissements d'avenir" program, reference ANR-19-P3IA-0001 (PRAIRIE 3IA Institute) and reference ANR-10-IAIHU-06 (Agence Nationale de la Recherche-10-IA Institut Hospitalo-Universitaire-6).


Frédéric Courteille

Florent Bourgeois

Michael Clifton

Martine Meireles

3D RECONSTRUCTION OF ULTRAFILTRATION CAKES FROM BINARISED IMAGES

Keywords: Ultrafiltration, Cake microstructure, 3D reconstruction

Transmission electron microscopy (TEM) imaging is a common technique for studying microstructure of ultrafiltration cakes. It yields local microstructural information in the form of two-dimensional gray-level images of slices a few particle diameters in thickness. This work presents a simple yet powerful particle-by-particle reconstruction scheme for simulating ultrafiltration cake microstructure from TEM images. The scheme uses binarised TEM images, thereby permitting use of lesser quality TEM images. Moreover, the reconstruction scheme embeds the multi-scale nature of ultrafiltration cake microstructure using a novel concept that consists in matching the morphology of simulated and measured structures at a number of scales identifiable in the microstructure. In the end, reconstructed microstructures are intended for bettering our understanding of the relationships between cake morphology, ultrafiltration performance and operating conditions.

Problem statement

Water purification using ultrafiltration membranes is an industrial-scale process that replaces both coagulation and sand filtration. At present, stakes are high for increasing the cost-effectiveness of the process. One key dimension to this problem is the correlation between ultrafiltration efficiency (production rate and filtrate quality) and the microstructure of the cake that forms on the membrane surface during the process. As in other filtration processes, cake microstructure is a crucial element in this problem. But ultrafiltration differs from other filtration processes in that direct observation of the morphology of an ultrafiltration cake is particularly difficult. TEM images are used to capture local information about ultrafiltration cake microstructure. The gray-level information of TEM images is potentially valuable for spatial reconstruction. However, it turns out that the quality of TEM images, when imaging nano-structures is often insufficient for making quantitative use of the gray-level information. Hence, there is value in developing a reconstruction scheme that uses binarised TEM images only. Such images correspond to the projection of the particles captured by TEM images onto a plane situated anywhere inside the cake slice captured by the TEM scan. A typical TEM image, with its binarised transform is shown in Fig. 1 for an ultrafiltration cake obtained using 26 nm monodispersed spherical particles. The thickness of the cake slice captured in this TEM image is 80 nm, which is typical of what can be expected from current TEM machines.

3D reconstruction scheme

In order to reconstruct ultrafiltration cakes from cross-section images, several options can be selected. At one end of the spectrum of reconstruction solutions, one finds statistical models such as the Cox process used by [START_REF] Jeulin | Multi-scale simulation of random spheres aggregates -Application to nanocomposites[END_REF]. With such an approach, Jeulin and Moreau were able to extract the parameters of several Boolean processes, which they combined to reconstruct microstructures with aggregates. In general, statistical models can be used provided particles are arranged according to some process whose spatial statistical properties can be extracted from morphological functions measured on 2D images. This approach is elegant and attractive both computationally and analytically. At the other end of the spectrum, one may consider reconstruction schemes which we shall refer to as global. For example, [START_REF] Oberdisse | Structure of interacting aggregates of silica nanoparticles in a polymer matrix: small-angle scattering and reverse Monte Carlo simulations[END_REF] have used the reverse Monte-Carlo (RMC) technique for simulating fractal aggregates. Although they did not start from images, as they used small-angle neutron scattering (SANS) information instead, their procedure starts from a randomly formed aggregate, which they rearrange using the RMC algorithm until the SANS spectrum of the reconstructed aggregate matches the measured one. In principle, global techniques can be applied to structures larger than single aggregates; however the task becomes more complex as the number of particles increases. Moreover, choosing the number of unknowns is a critical step with such reconstruction schemes, as the reconstructed volume must be filled with particles before the start of the rearrangement process. The particle rearrangement algorithm is computer intensive; however there is no reason why such an approach could not be made to work satisfactorily. The procedure developed in this paper sits somewhere in between purely statistical and global approaches. The former was not selected on the basis that there is no evidence that ultrafiltration cakes conform to any particular statistical configuration. The authors avoided using the latter approach because of foreseeable difficulties with simulating large volumes. The proposed reconstruction scheme is an incremental particle-by-particle reconstruction of the cake slice that corresponds to the volume captured by the TEM image. Adding particles one by one avoids having to specify the number of particles a priori and circumvents the difficult choice of this unknown. At every iteration, i.e. every time a new particle is added into the system, the point P(x,y,z) at which it is implanted is chosen such that it minimizes the difference between the specific surface area of the real microstructure and the specific surface area of the reconstructed microstructure at the current iteration with the new particle added at point P. In practice, the new particle can be positioned at random locations of two types. If the targeted specific surface area is less than that measured at the end of the previous particle addition step (see Fig. 2), the particle is positioned at a randomly chosen location such that its projected area does not overlap with that of other particles; hence it contributes its whole perimeter to the specific surface area. Otherwise, as shown in Fig. 3, it is positioned at a randomly chosen location such that its projected area overlaps with that of other particles, thereby contributing only marginally to the change in specific surface area. The specific surface area can be calculated directly from the image by counting pixels that are located at the interface between the porous and solid phases. However, the function that is used in this work is the classical two-point covariance function S [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF], which is defined over the set of chords of length h. Not only S is a reliable means for calculating specific surface area from binary images [START_REF] Berryman | Measurement of spatial correlation functions using image processing techniques[END_REF], but it also captures a fine morphological signature of the texture. The authors have found that, despite the fact that the proposed reconstruction algorithm relies on matching the specific surface area only, the two-point covariance function of the textures reconstructed with the proposed method matches that of the original textures over significant chord lengths. The reason for this interesting finding is believed to be the combination of the multi-scale reconstruction scheme that is presented later, and the specified particle shape, which constrains the covariance function significantly. 2. The targeted specific surface area is less than that measured at iteration n, hence particle P n+1 is positioned at iteration n+1 at random inside the reconstructed volume such that its projected area does not overlap with that of other particles. S and n S ɶ correspond to the two-point covariance function of the targeted image (solid line) and the reconstructed structure (dotted line) after n iterations respectively.

S(h) ( ) n S h ɶ chord length h Figure
S(h) ( ) n S h ɶ chord length h
Figure 3. The targeted specific surface area is higher than that measured at iteration n, hence particle P n+1 is positioned at iteration n+1 at a position such that its projected area overlaps with that of other particles. S and n S ɶ correspond to the two-point covariance function of the targeted image (solid line) and the reconstructed structure (dotted line) after n iterations respectively.

The next issue consists in capturing the multi-scale nature of the texture from the image of the texture, from single particles to large aggregates. Quantification of the multi-scale information of the texture is achieved in two successive steps, as described in Fig. 4: step one, a set of images of individual scales is generated by successive application of an erosion operator1 ; step two, the 2-point covariance function is computed for each scale. The multi-scale texture simulation algorithm then simply consists in matching the two-point covariance function at every individual scale, starting from the upper scale (rightmost image in Fig. 4) and moving to the left of Fig. 4 towards the lower scales. Addition of particles reduces the porosity, hence the two-point covariance function S(h) moves upward every time a particle is added to the system. Once the porosity of the current scale is matched, the algorithm uses the slope of the two-point covariance function of the next scale as the next target for the particle-by-particle addition scheme presented earlier.

Results and analysis

Firstly, the simulation algorithm was tested against images obtained from slices through an artificial cake generated using the physico-chemical model developed by [START_REF] Madeline | Restructuring of colloidal cakes during dewatering[END_REF]. This is an ideal situation, in that particles in the images are spherical and every spatial property is known for the particles inside the cake, making it possible to develop and fully test the validity of the reconstruction scheme. Fig. 5 shows the results obtained when simulating with the proposed algorithm the texture whose frame is bolded. Fig. 5 was deliberately presented as a children's "Spot the difference" game, with the original cake slice image placed amongst 8 simulations. Should the bold frame be removed, it then becomes rather difficult to spot the original image in Fig. 5, thereby confirming the potential of the proposed 3D simulation technique for ultrafiltration cake slice simulation. We note however some differences with the reconstruction of the aggregates. The cause is an insufficient sampling of the larger particle aggregates in the original image, due to the fact that the technique used to generate the artificial ultrafiltration cake could only produce small cake volumes.

Figure 5. Projected image of the artificial ultrafiltration cake slice (bold frame) and 8 projected images of simulated slices using the proposed reconstruction method.

The second application of the proposed method is the reconstruction of the binarised TEM image of a real ultrafiltration cake, which was presented in Fig. 1. Contrary to the artificial cake reconstructed in Fig. 5, the sampling of the aggregates is quite satisfactory, as the field of view is large relative to the size of the largest aggregates. One reconstruction result is shown in Fig. 6, directly to the right of the original image. A notable difference with the previous ideal example is that particles no longer appear as spheres in the raw image. This is due to both the size of the particles (26 nm nominal diameter) and the binarisation step applied to the gray-scale TEM image. Nevertheless, the projected image of the simulated slice gives a satisfactory match to the original image. Moreover, as we hinted earlier, Fig. 6 confirms the match between the simulated and measured 2-point covariance function for a number of intermediate scales. Visually, all reconstructed images are satisfying. One observes that the simulation scheme recovers correctly both the void space and the state of particle aggregation. The state of aggregation is better recovered than with the artificial cake used earlier as the raw image provides a good sampling of the morphological features that correlate over long distances.

S(h) h (in pixels)

Quantitatively, the quality of the convergence of the distance watershed objective function can be assessed using chi-square statistics, which shows equality between measured and reconstructed objective functions with a confidence greater than 95%. In the end, the interfacial area is precisely matched, as confirmed by the match between the slopes of the 2-point covariance functions near the origin. Some improvement to the particle-by-particle addition scheme could be made in order to match the covariance functions over a longer range of chord length h. Indeed, because particles of given shape and size distribution are added into the system, the shape of the resulting covariance function is strongly constrained. This explains why the simple particle-by-particle addition method embedded into the proposed method yields in the end a good match of the measured 2-point covariance function even at relatively long chords. The reader should also notice that the simulated covariance function is not expected to match the measured one at long chord lengths since the uncertainty of the measurements is directly proportional to the chord length h. Hence, the authors believe that a further improvement would be required should the technique be applied to textures whose scales cover several orders of magnitude. Fig. 7 gives a 3D view of one of the cake slices whose projection was shown in Fig. 5. The reconstruction scheme described here uses binarised TEM images only, which do not contain information that can be used directly to derive particle position in the vertical direction. Without additional constraint, the reconstruction algorithm yields a particle network that is looser than the real one. This problem is successfully overcome by using the abovedescribed algorithm to position particles in the X-Y plane, and then selecting vertical positions that favour inter-particle contact in 3D when possible (i.e. when the particle's projection is not isolated in the projected image of the slice under reconstruction). This constraint embeds the physics of particle deposition in an ultrafiltration cake, which dictates that particles must form a connected network. The authors have checked from artificial cakes (See earlier example illustrated in Fig. 5) that the number of particles that are implanted in the reconstructed volume is very close to the number of particles inside the original cake slice. One explanation is that ultrafiltration cakes are very dilute systems, so almost every particle inside a slice a few particle diameters in thickness will contribute something to the projected binary image of that slice. In other words, only very few particles can possibly hide directly behind others inside the thickness of a cake slice captured by TEM. Nevertheless, the reconstruction scheme presented in this paper does not yield the exact number of particles inside the slice, nor can it possibly decide what the best vertical position might be for each particle. Before concluding, it should be stressed that a major advantage of reconstructing 3D textures is that it permits generation of many different random realisations, thereby providing confidence intervals for any property that can be derived from microstructure simulations. This is a particularly important practical point as it is often the case that only one single ultrafiltration cake image is available for a given set of experimental conditions.

Conclusions and perspectives

This article has presented a new scheme for simulating ultrafiltration cakes from binarised TEM images. Such images can be obtained routinely in practice even from poor quality gray-level images: this makes the approach particularly valuable for microstructure studies of nanostructures.

The originality of the proposed 3D simulation scheme is twofold. Firstly, the scheme uses a particle-by-particle addition process. Hence, the 3D simulation does not rely upon any assumption about the spatial statistics of the system, and simulated volumes are not limited in size as with global numerical reconstruction schemes. Secondly, the proposed simulation scheme simulates multi-scale textures successfully by matching the morphology of simulated and measured textures at every individual scale. The morphological properties of the texture at every scale are quantified using the 2-point covariance function. Because the scheme adds particles of given shape and size distribution into the reconstruction volume, it is found that it is sufficient to match the yintercept and the slope near the origin of the 2-point covariance function in order to obtain satisfactory 3D simulations. Hence, the proposed technique is quite straightforward to implement. Convincing simulation results were obtained with both ideal and real ultrafiltration cakes.

In a forthcoming publication, the authors will present a generalized version of this 3D simulation algorithm that solves the problem of reconstructing a whole cake from one or more TEM images. The principle consists in using a moving reconstruction volume whose thickness equals that of the TEM image, centred on the vertical position of the particle that is added during the particle-by-particle addition reconstruction scheme. Hence, a volume larger than the simulated slice is built in the end, in such a way that every slice centred on every particle inside the system has a projected image whose 2-point covariance function matches the measured values. Using an elevation-invariant covariance function, or an elevation variant one, the technique is able to reconstruct not just one single cake slice, but a whole ultrafiltration cake. Moreover, the technique that is developed here with spheres can readily be applied using other particle shapes (e.g. polyhedra) and distributions of particle size. Hence, with such a technique, which is quite simple to implement, one is able to simulate a wide variety of multi-scale two-phase materials. It is worth noting that the 3D simulation technique presented here is also extendable to multiphase systems, using inter-phase 2-point covariance functions.

Figure 1 .

 1 Figure 1. Original and binarised TEM image of an ultrafiltration cake (image size ≈ 1 µm 2 )

Figure 4 .

 4 Figure 4. Illustration of the multi-scale image analysis scheme.

Figure 6 .

 6 3D simulation of a real ultrafiltration cake slice (original binarised TEM image to the left; simulated image to the right). Measured 2-point covariance functions are represented as a solid line, whereas simulated ones are shown as dotted lines.Confidence intervals are given for 300 simulations.

Figure 7 .

 7 Figure 7. Example of a 3D ultrafiltration cake slice simulated with the proposed method.

Erosion is one of the two elementary operators in the area of mathematical morphology, the other being dilation. The basic effect of the operator is to erode away the boundaries of regions of foreground pixels (black pixels).


José Pérez-Escobedo

Luis

Morales-Mendoza Fernando

Aguilar-Lasserre Alberto

Azzaro-Pantel Catherine

Pibouleau Luc

Analysis of new pharmaceutical product mix selection with discrete-event simulation

Keywords: New Product Development, Discrete-Event Simulation, Pharmaceutical, Portfolio management

This paper is devoted to the presentation of a discrete-event (DES) simulation model to represent the drug development process and to help decision-making. The presentation is mainly focused on the development of the simulator initially developed for scheduling, planning and design purposes. A typical example involving 6 drugs is used as a test bench to illustrate the proposed methodology. Monte Carlo trials are used for each portfolio of drugs. Each trial can be viewed as a simulated path of the potential product along the probabilistic network. Different scenarios are analyzed, reflecting the pessimist or optimistic viewpoint of the decision's maker.

Introduction

A fundamental challenge in managing a pharmaceutical or biotechnology company is identifying the optimal allocation of finite resources across the infinite constellation of available investment opportunities. In that context, the optimal management of the new product pipeline has emerged at the forefront of all strategic planning initiatives of a company. This issue is traditionally identified as a complex one since it integrates various areas such as product development, manufacturing, accounting and marketing, with a large variety of parameters and decisionmaking levels. A strategic investment plan should simultaneously address and evaluate in a proper manner the following four main issues: product management, clinical trials uncertainty, capacity management and trading structure. It is also generally viewed as a multistage stochastic portfolio optimization problem. This paper is devoted to the presentation of an architecture that combines an optimization procedure based on a genetic algorithm and a discrete-event (DES) simulation to assess the uncertainty present in the pipeline and to help decision-making. The presentation is mainly focused on the development of the simulator initially developed for scheduling, planning and design purposes. The objective of the simulation section is to evaluate some performance indicators that may be useful for the decision's maker, for instance the expected net present value (NPV) computation. A typical example involving 6 drugs is used as a test bench to illustrate the proposed methodology. Monte Carlo trials are used for each portfolio of drugs. Each trial can be viewed as a simulated path of the potential product along the probabilistic network. Different scenarios are analyzed, reflecting the pessimist or optimistic viewpoint of the decision's maker.

System description and model formulation

Basically, three stages are involved in the life cycle of a pharmaceutical product: discovery, development and commercialization. More precisely, the tasks involved during the New Produc Development (NPD) to market (e.g. development, manufacture, clinical trials) require resources (e.g. capital, in-house capacity, personnel, contract manufacturing capacity). Our idea was to use the potential of discrete-event simulation to model the series of decision points along the drug development pathway. For example, at the end of each phase of clinical trials the probability of clinical success resulted in go/no-go decisions. Further decision points include those made at each manufacturing stage which could be carried out inhouse or by a contract manufacturer; this was determined as a function of the available capacity at a particular time. Figure 1 shows the main decisions and some examples of the options available for a drug process development, established from the work of [START_REF] Rajapakse | Integrated approach to improving the value potential of biopharmaceutical R&D portfolios while mitigating risk[END_REF]). A set of potential products are given, each of which must undergo a set of testing tasks. Each task has an associated duration, cost, and probability of success. Given the income for each product as a function of the time of product introduction, the problem is to model tasks scheduling while computing the traditional Net Present Value (NPV) economic criterion.

From batch plant design to project management by Discrete-Event simulation

A discrete-event approach was then adopted to model the pharmaceutical development pathway. This technique has been identified as an efficient method to study the impact of business and process decisions considered during the management of pharmaceutical product development on the portfolio profitability. In a DES, a process is described as it evolves with time and changes take place only a finite number of times, i.e. event occurrence date. The DES was developed using C++ object-oriented language, keeping the approach proposed by [START_REF] Bérard | Towards an incremental development of discrete-event simulators for batch plants : use of object-oriented concepts[END_REF]. The tool seeks to integrate various aspects, including resource management and the development of manufacturing activities required. A conceptual model is derived from the generic product development pathway presented in Figure 1. In [START_REF] Bérard | Towards an incremental development of discrete-event simulators for batch plants : use of object-oriented concepts[END_REF], a four layer framework was proposed based on the following items {engine, event, object, supervisor}, the aim being the development of a standard library for the simulator classes that are general to any case, thus minimizing the task of treating different study cases or the variants of a given one (i.e. design or scheduling objectives). In this approach, at the lowest level, the common engine can be found. Initially, the events in the next level are generic events common to all batch plant simulations: in this case, the definition must be adapted since we have consider the whole life cycle of a project related to a product. In the same way, the objects taken into account present some similarities but differ in their appreciation: for instance, in batch plant scheduling problems (BPS), material resources are constituted by equipment whereas in NPD problems, resources may be viewed more globally. In fact, the main differences at this step occur from a terminology point of view and this can be easily transposed in the NPD formulation (see Table 1). The core of the simulator is the Engine, which has two functions: the former is to order the Events in its Calendar by their occurrence date whereas the latter is to activate them if the necessary resources are available; if not, it reports the Event to a next Date. As previously said, an Event represents a change of the real system at a given time. The class Event is a basis class from which the different Events must be defined. If resources for this activity are available, the Event is activated; conversely, if resources are not available, the activity will be scheduled later. An Event is characterized by its occurrence date, its action over the system and a type that enables to give priorities when two or more Events have the same occurrence date. As a general rule, Events which release resources have priority over the others, and when Events have the same type, the classical FIFO rule (First In First Out) is applied. This will be useful when different projects compete for the same resources. The Event Class previously developed was generic enough to embed the NPD formulation. An important class is the so-called Stage Class that involves the activity of a product at a given stage requiring a given amount of resource and time. The net present value (NPV) of each product candidate is generated as a function of the year and computed at several steps. Constrained resource availability is considered in the model : it means that if two tasks compete for the same resource with a limited level, a decision rule must be applied. The traditional FIFO (First In-First Out) rule is applied first. In case of equality, priority is given to the drug with the highest market value. It is also assumed since the duration of each task is never shorter than one year that the capital resource is renewable (yearly budget allocation). The probability of a candidate failing a task at the so-called Probability stage Class results in cancellation of the development process for that candidate.

Figure 2 Relationship between classes -system activities

Case study

The case study is an adaptation of the example presented in [START_REF] Rajapakse | Integrated approach to improving the value potential of biopharmaceutical R&D portfolios while mitigating risk[END_REF]. It involves a portfolio of six potential monoclonal antibodies (MAbs) that are ready for clinical development. Given the finite level of resources available, the simultaneous development of the six products can not be carried out. The interest of the simulation approach is to provide decision-makers with an explicit view of the best product mix in order to satisfy performance indicators. Probabilities were assigned to reflect the level of uncertainty related to the considered drugs (see Table 1). They are representative of the kind of information available from experts in the pharmaceutical industry. The probability distributions assigned to key risk factors for drugs relevant with differing levels of uncertainty are given in Table 2. The phase transition probabilities reported in [START_REF] Reichert | Monoclonal antibodies in the clinic[END_REF] from a large data collection and analysis for 186 antibodies entering clinical studies were considered for this case study and are shown in Table 3. The chance of failure due to technical reasons was set again according to the uncertainty level assigned to a particular drug candidate with those deemed to be high risk having a higher chance of failure due to technical reasons (Table 4). Drug B was modelled as a new indication for an existing drug and hence there was no technical uncertainty associated with it. Three scenarios are considered based on level of resources for drugs development (see Table 5). Scenario 1 considers all six drugs in the portfolio and unlimited levels of resources, thus establishing a base case against which to compare the results from resource-constrained simulations. Next, two resource levels of $500M and $750M were defined as constraints with a lower number of drugs to be developed.

Simulation implementation and results analysis

To take into account the stochastic nature of the problem, the simulation is repeated 350 times selecting random sampling values. The computation of the Net Present Value (NPV) is cumulated along the evolution of a project for a given product or a sequence of products. The considered economic criterion is then based on the expected value (ENPV) for the 350 simulations relative to each sequence. The risk associated to a sequence is evaluated by two indicators: the former, (Risk 1 ), is based on the number of times that the NPV is negative for the products in a given sequence and the latter, (Risk 2 ), considers the NPV standard deviation. The simulation was performed for the set of 36 sequences as it was defined in [START_REF] Rajapakse | Integrated approach to improving the value potential of biopharmaceutical R&D portfolios while mitigating risk[END_REF]. The results displayed in Figures 3 and4 represent the ENP versus (Risk i ). It can clearly be seen that two sequences are particularly interesting for the decision maker. Both involve a sequence of three products (ACD for portfolio 5 and ABD for portfolio 2). The results thus show that diversifying the portfolio is not recommended to generate profits: the 6-product sequence exhibits low performances. This behaviour can be explained by the high number of failures of products at trial phases as shown in Figure 5, representing the number of products reaching the commercialization phase for the 350 simulations. This trend is confirmed by the computation of the corresponding NPV (see Figure 6, for instance, 700 times). It must be yet observed that the NPV standard deviation is relatively high, which is due to the typical adopted cost values relative to each step of the process. These trends are confirmed by the results presented in Table 6, in which the confidence limits for each sequence are computed from the ENP viewpoint.

Conclusion

This study has presented a simulation approach to evaluate product portfolio selection in the development phase of pharmaceutical companied. The results show that simulations can help decision-aid to identify clearly where the company portfolio is positioned with regard to risk-return characteristics of alternative product portfolios and hence make appropriate investments based on limited resources. The risk is evaluated by two criteria: the former is a ratio measuring the number of sequences leading to negative NPV values relative to the total number of the performed simulations for a given sequence; the latter involves the standard deviation between the NPV values for an evaluated sequence. The flexibility of the simulation tool allows differences in product development to be captured easily.

Figure 5

 5 Figure 5 Frequency according to the number of products for the sequence of six products

Table 1

 1 General information about drugs in the portfolio.

	Drug candidate	MAb type	Market share	Uncertainty
	A	Chimeric	High	Medium
	B	Humanised	Low	Very low
	C	Chimeric	Medium	Medium
	D	Murine	Very high	High
	E	Chimeric	Low	Medium

Table 2 .

 2 Risk factors and probability distributions for drugs with differing levels of uncertainty

		High uncertainty	Medium uncertainty	Low uncertainty
		drug			drug	drug
		Value Probability Value	Probability	Value Probability
		12	20	9	20	6	10
	Ph I development time (Months)	18 24	30 50	12 15	50 30	9 12	70 20
		12	20	9	20	6	10
	Ph II Manufacturing time (Months)	18 24	40 30	30 50	30 50	9 12	70 20
		18	20	12	20	12	10
	Ph III Clinical trial time (Months)	24	30	18	50	15	70
		36	50	24	30	18	20
		750	10	1125	30	1200	20
	Market manufacturing cost (COG)	1000	20	1500	50	1600	70
	($ g-1)	1250	70	1875	20	2000	10
	Product yield (%)	30	40	30	30	30	10
		50	30	50	50	50	20
		70	30	70	20	70	70
	Delay in material delivery (Months)	0	20	0	50	0	70
		3	30	3	30	3	20
		6	50	6	20	6	10

Table 3

 3 General information about drugs in the portfolio.

	Monoclonal antibody	Phase I to II	Phase II to III	Phase III to	Review to
	type	(%)	(%)	review (%)	approval (%)
	Murine MAbs	77	52	45	33
	Chimeric MAbs	86	4	80	100
	Humanised MAbs	84	725	75	100

Table 4

 4 General information about drugs in the portfolio.

	Drug	Probability of failure due to	Market value
	candidate	technical reasons at Phase III (%)	($M)
	A	20	523
	B	0	427
	C	10	410
	D	30	1029
	E	5	31
	F	15	49

Table 5

 5 Resources levels and drugs in the portfolio

	Resource level ($M) Number of drugs per portfolio Number of possible portfolios
	Unconstrained	6	1
	500	3	20
	700	4	15


Mary Hanhoun

Catherine Azzaro-Pantel

Beatrice Biscans

Michele Freche

Ludovic Montastruc

Luc Pibouleau

Serge Donenech

REMOVAL OF PHOSPHATE FROM SYNTHETIC WASTEWATER BY STRUVITE PRECIPITATION IN A STIRRED REACTOR

Keywords: Struvite, Precipitation, P-recovery, Stirred reactor, Modelling

The cost of environmental protection and pollution prevention is increasing, above all because of stringent effluent quality standard. In that context, phosphorus impact on water pollution plays a major role because it accelerates eutrophication. This work addresses the problem of phosphorus recovery from synthetic wastewater by struvite precipitation. The general objective is to develop a methodology for an optimal design of an automated pilot reactor for the treatment of wastewater. In this perspective, preliminary precipitation experiments were carried out at laboratory scale to examine the effects of final pH, magnesium and phosphate concentrations in the influent. A thermodynamic model for struvite precipitation is proposed involving Davies activity coefficients.

INTRODUCTION

Wastewater discharges of nitrogen and phosphorus to the environment have produced an increase in water pollution because these nutrients accelerate eutrophication, producing grave consequences for aquatic life as well as for water supply for industrial and domestic uses. One proposed solution to this problem is the recovery of nutrients using precipitation. Two major precipitations have been developed for phosphorus recovery from wastewater, involving either the so-called calcium phosphate (CP) or magnesium ammonium phosphate (MAP), i.e., struvite. Struvite is a white crystalline substance consisting of magnesium, ammonium and phosphorus in equal molar concentrations (MgNH 4 PO 4 -6H 2 O). Struvite forms according to the general reaction shown below [START_REF] Bouropoulos | Spontaneous precipitation of struvite from aqueous solutions[END_REF]:

Mg 2+ +NH 4 + +PO 4 3-+6H 2 O MgNH  4 PO 4 , 6H2O
Struvite crystallisation from wastewater effluents is seen as an alternative to traditional biological and chemical phosphorus removal processes widely used in the wastewater treatment industry. It presents the advantage of not only removing phosphorus but also generating a compound that could be reused as a commercial fertilizer [START_REF] Durrant | Review of the feasibility of recovering phosphate from wastewater for use as a raw material by the phosphate industry[END_REF][START_REF] Nelson | Struvite precipitation in anaerobic swine liquid: effect oh pH and Mg:P ratio and determination of rateconstant[END_REF]. Struvite formation occurs relatively quickly because of the presence of excess supersaturation in the liquid, as a result of the chemical reaction of magnesium with phosphate in the presence of ammonium. Supersaturation may be developed by increasing the aqueous medium content in ammonium, magnesium or orthophosphate and/or pH. For this study, batch experiments were conducted to determine the efficiency of struvite precipitation in synthetic wastewater. More particularly, the influence of the following parameters was investigated: pH, magnesium and phosphate concentration. The experimental results are compared with those obtained by a thermochemical model, enough representative of struvite precipitation using a Davies activity coefficient modeling [START_REF] Hanhoun | A thermochemical approach for struvite precipitation modelling from wastewater[END_REF]. This is a first step towards the development of a model taking into account kinetics of struvite precipitation.

MATERIALS AND METHODS

A series of batch investigations was conducted to study the influence of species concentration and pH on struvite precipitation. Experiments were carried out at 25°C in a 500 ml stirred vessel. Stock solutions of magnesium and ammonium phosphate were prepared from corresponding crystalline solids MgCl 2 .6H 2 O and NH 4 H 2 PO 4 . Deionised water was used to prepare the synthetic wastewater solution. The supersaturated solutions (corresponding to a final phosphorus concentration of 4 mM with an Mg/NH 4 /PO 4 molar ratio equal to 1) were prepared by a rapid agitation of NH 4 H 2 PO 4 and pH was then adjusted (by the addition of the appropriate amount of a standard solution of sodium hydroxide), followed by the addition of the appropriate volume of stock magnesium chloride solution. Finally, pH was adjusted again as needed. Experiments were performed at several pH values [7-11.5] and keep constant all over the period of the reaction (1 hour). The homogeneity of the solution was performed by a magnetic stirrer. Other experimental runs were carried out in the same way with a final phosphorus concentration of 20 mM. The formed precipitate was collected by filtration through a 0.2 mm membrane filter and dried at ambient temperature. Residual concentrations of Mg 2+ (respectively PO 4 3-

) ions in solution after filtration were analysed by atomic absorption (AA) (respectively by spectrophotometry). X-Ray Diffraction (XRD) was used to determine the composition of the precipitate produced. Further samples of precipitate were dissolved in nitric acid: these samples were analysed for magnesium and phosphorus concentrations. The theoretical concentrations were calculated from the dissolved mass assuming that the formed precipitate is pure struvite.

Results

Figure 1 shows the results of an X-Ray Diffraction analysis performed on the formed precipitate. The presence of struvite is indicated by the location of the intensity peaks, corresponding to the standard database lines for struvite. The results show that only the peaks of struvite were found in the investigated pH range [7-10.5]. They were further confirmed by comparison with solid analysis of phosphate/magnesium. It must be yet pointed out that other precipitates may occur as reported in [START_REF] Bouropoulos | Spontaneous precipitation of struvite from aqueous solutions[END_REF] with pH values above 10.5, where a visible floculation of crystals is observed with the precipitation of Mg 3 (PO 4 )2.22H 2 O. Identical phenomena were confirmed within the same pH range in other studies [START_REF] Lee | Removal of nitrogen and phosphate from wastewater by addition of bittern[END_REF][START_REF] Corre | Impact of Reactor Operation on Success of Struvite Precipitation from Synthetic Liquors[END_REF].

CHEMICAL EQUILIBRIUM MODEL FOR STRUVITE PRECIPITATION

In this work, a thermodynamic model for struvite precipitation is proposed, based on Davies model of activity coefficients taking into account various phases which can precipitate in the range of pH to be considered. The model determines the phosphate conversion rate and is based on numerical equilibrium prediction of the study system Mg-NH 4 -PO 4 -6H 2 O. The parameters include the solubility product of struvite. The model involves, on the one hand, the mass balances for magnesium, ammonium and phosphate as a function of conversion rate and the electroneutrality equation, on the other hand [START_REF] Montastruc | A thermochemical approach for calcium phosphate precipitation modeling in a pellet reactor[END_REF]. As reported elsewhere, the thermodynamic model [START_REF] Hanhoun | A thermochemical approach for struvite precipitation modelling from wastewater[END_REF] uses the following equilibrium constants (see Table 1)

RESULTS AND DISCUSSION

The experimental results have been compared with the values predicted by the thermodynamical model for pH several values [7-10.5], as displayed in Figures 2 and3. It can be seen that the experimental values are in agreement with the computed results. Optimum experimental pH for struvite precipitation is 9.8, which corresponds to a phosphorus removal of 80% (respectively 94.8%) for an initial concentration of 4 mM (respectively 20mM). This experimental result is identical to the value reported in the literature [START_REF] Booker | Struvite formation in wastewater treatment plants: opportunities for nutrient recovery[END_REF].

It can be thus concluded that the developed model allows to predict thermodynamic equilibrium of struvite production for a pH range between 7 and 10.5. It will be used as a part of the thermodynamic core model for the studied precipitation process. The process design yet requires knowledge of struvite crystal growth and an establishment of kinetic models. For this reason, another experiment was carried out in a stirred reactor (1L), with an initial concentration of PO 4 , Mg, NH 4 of 4 mM, to study the effect of initial pH on struvite precipitation kinetics. During this experimental run, pH was not set at a fixed value as previously: only initial pH was controlled [8.5-10] and its evolution with time was measured (see figure 4). The measurement of Mg and P final concentrations exhibits a molar ratio approximately equal to unity, thus demonstrating the exclusive formation of struvite for each experiment. The phosphorus residual content was also compared with the value computed by the model as displayed in Figure 4. It must be highlighted that the observed induction time was affected by variations in initial pH. For example, for an initial pH value equal to 9, the solution became cloudy after 6 minutes, whereas the first particles were visible only after 20 seconds at pH 10. The result presented here showed that pH is an efficient measurement that can be used for the prediction of the kinetics of struvite precipitation. The pH drop that was typically observed during struvite precipitation reflects the rates at which struvite crystals may nucleate.

CONCLUSIONS AND PERSPECTIVES

This research investigated the potential of struvite precipitation as a method to recover phosphorus from wastewater, while concurrently removing ammonium and magnesium. For high pH values [7-10.5], a domain in which struvite precipitation can occur is observed. For this purpose, this study has presented a thermodynamic model for struvite precipitation: the obtained results have shown that an increase in pH enhances precipitation efficiency, with a predicted maximum precipitation rate at a pH value equal to 9.8. The model validity was checked for a pH range between 7 and 10.5, with a strong hypothesis based on a struvite exclusive formation. This assumption is validated by XRD analysis for pH-constant experiments, and for free pH experiments by Mg and PO 4 measurements.

The thermodynamic model developed in this study will now be embedded in a modelling framework for determination of process operating conditions for struvite precipitation in a stirred reactor.

Figure 1 .

 1 Figure 1. X-ray diffraction analysis results of formed precipitate from the synthetic solutions at different pH values and comparison with values from struvite database.

Figure 2 .Figure 3 .

 23 Figure 2. Comparison between the experimental values and the results of the thermochemical model. (Low concentration 4 mM)

Figure 4 .

 4 Figure 4. Average pH versus time for initial pH values

Table 1 .

 1 Equilibrium constants for complexes presented in the equilibrium equations.[START_REF] Hanhoun | A thermochemical approach for struvite precipitation modelling from wastewater[END_REF].

	Equilibrium		Value
	H 2 PO 4 HPO 4 2--+H + ↔H 3 PO 4 +H + ↔H 2 PO 4 -PO 4 3-+H + ↔HPO 4 2-H 2 PO 4 -+Mg 2+ ↔MgH 2 PO 4 HPO 4 2-+Mg 2+ ↔MgHPO 4 PO 4 3-+Mg 2+ ↔MgPO 4 -OH -+Mg 2+ ↔MgOH + H + +NH 3 ↔NH4 + H 2 O↔H + +OH -	+	-2.15 10 -7.20 10 -12.35 10 -0.45 10 -2.91 10 -4.80 10 -2.56 10 -9.25 10 10 -14
	NH4 + +PO 4	3-+Mg 2+ ↔MgNH 4 PO 4 .6H2O 10 -13.26