H Z Randrianandraina

Y Cressault

A Gleizes

Improvements of radiative transfer calculation for SF 6 thermal plasmas

Keywords: thermal plasmas, radiation, radiative transfer, SF 6

We present a comparison between an exact calculation of radiative transfer in SF 6 thermal plasma based on a fine description of the spectrum with 300 000 spectral points for each temperature value, for simplified conditions (1D and 2D geometries with imposed symmetrical temperature profiles and Local Thermodynamic Equilibrium) and two kinds of approximated calculations. The first is the classical Net Emission Coefficient largely used in arc modelling. The second one is based on a very simplified spectral description with only seven intervals assuming a grey body condition within each interval and using the Planck and the Rosseland averaging for deducing the mean absorption coefficient. At high temperature the use of the Rosseland averaging is not satisfactory. The other two approximated methods (Net emission and the Planck averaging) are acceptable but the radiative flux is in general not very accurate. The radiative transfer calculation can be improved following three ways: a better knowledge of the basic processes and in particular of the absorption coefficients for diatomic and polyatomic molecules; the use of different definitions of mean absorption coefficients (Planck averaging at high temperature and mean natural value at low temperature); a careful choice of the spectral intervals.

-Introduction

Radiation plays an important role on the behaviour of electric arcs and of thermal plasmas and is often the main mechanism of energy transport within the plasma and between the plasma and the surrounding materials. In the case of SF 6 plasmas, radiative transfer has been calculated for more than 30 years because this gas is widely used in High Voltage Circuit Breaker (HVCB) where radiation has three important effects: the radiation losses tend to cool the plasma during the high current phase; the partial absorption of radiation by the cold surrounding SF 6 gas heats the gas and thus increases the pressure within the closed apparatus which has a big influence on the interruption capability; finally the radiative flux is highly responsible for the nozzle ablation which is used in modern HVCB to increase the efficiency of the interruption. The difficulty of radiative transfer calculation is related both to the geometry and to the spectral dependence of the radiation intensity. In practical conditions, exact radiative transfer calculation is impossible in a general modelling process with a lot of iterations. This fact obliges to make approximated calculations and the choice of these methods depends on the objectives.

The first objective is to compute the radiation escaping the hottest regions of the plasma which allows calculating with a good accuracy the temperature in these regions and then to have a realistic description of the plasma. The simplest and nevertheless quite accurate method is the previous calculation of the Net Emission Coefficient (NEC) followed by the direct use of this coefficient in thermal plasma modelling. We will present the bases of this method in part 3. This method is of interest because the complicated calculation of the net emission (net meaning emission minus absorption, which is mathematically represented by the divergence of the radiative flux) integrated over all the spectrum (VUV, UV, visible and IR parts) is previously calculated for isothermal conditions and then injected in the energy conservation equation by a simple source term. This NEC was initially proposed by Lowke [1] and first calculated for SF 6 by Liebermann and Lowke [2]. Since this pioneer work, other teams have been calculating the NEC for SF 6 [3][4] and several works demonstrated the interest in and accuracy of this coefficient for calculating the temperature field in arcs [5] introducing a multiplying factor [6][7].

However, by definition, the NEC cannot give an estimation of the radiation absorption in the plasma and cold gas surrounding the hot core of the arc. The first idea to take into account this absorption was to modify the NEC introducing an empirical absorption adjusted with comparison with some experimental results [8]. A more sophisticated approach was proposed by Sevastyanenko introducing the Partial Characteristics Method (PCM) explained in [START_REF] Soloukhin | Radiative Heat Transfer in High Temperature Gases[END_REF]. This method is based on a previous and complicated calculation of source and sink terms (called Som and ΔSim) computed for segments of given lengths and functions of the two extreme points temperatures, assuming a linear variation of temperature along one segment.

Then these data can be used for calculating the radiative transfer in a general modelling along various directions. For SF 6 the Som and ΔSim terms have been computed by Aubrecht and Lowke [START_REF] Aubrecht | [END_REF] and by Raynal et al [6; 11].

The NEC and PCM are based on a fine description of the spectrum needing a very high number of spectral points at each temperature value. On the contrary, there exist classical methods of radiative transfer calculation that are based on a simplified description of the spectrum (P1 model, discrete ordinates method, see for example the book of Modest [START_REF] Modest | Radiative Heat Transfer[END_REF]).

The most popular simple spectral description for SF 6 plasma is to define a few spectral intervals (less than 10 in general) where the absorption coefficient is constant for a given temperature (grey body approximation by spectral intervals or Mean Absorption Coefficient (MAC)). It is well-known that this description may lead to very large errors because the line transfer is not well treated. Nevertheless, Gleizes et al proposed an original way to compute the MAC by calculating the Planck average coefficient and by using an escape factor for the lines [START_REF] Gleizes | XXIII ICPIG Int. Conf. Phen. Ionized Gases Toulouse[END_REF]. For SF 6 the spectrum was divided in five intervals and the values of the MAC can be found in [START_REF] Erraki | Etude du transfert radiatif dans les plasmas thermiques : application au SF6 et au mélange argon-fer[END_REF]; they have been used by Eby et al in a HVCB SF 6 arc modelling [START_REF] Eby | [END_REF].

Having the possibility to use several methods, it is interesting to test them and to define the best compromise between the accuracy of the results and the difficulty (including CPU times) of the methods. In general, for SF 6 plasmas, the validation of the methods is realised either by comparing the modelling results (temperature field, variations of global pressure) obtained using various methods and/or by comparing these results with some experimental data. It is to be noted that a direct validation of the total flux or of the divergence of the flux is not possible in practical cases because a large part of the radiation occurs in the VUV part and is mainly caused by resonance lines. Only partial experimental validations on the limited part of the spectrum [16] and on the absorption by cold SF 6 [17] have been performed. Among the works including a comparison of simplified methods for radiative transfer calculation, we may refer to the papers of Verite et al [18] and Dixon et al [19], both including the NEC and PCM methods (MAC also in [19]). The indirect validation with temperature, voltage and pressure measurements, allows to define an uncertainty of the order of a factor of 2 for the divergence of the flux in the hottest regions (comparison with temperature in the hottest regions) and for the absorbed radiation (comparison with pressure). The accuracy on the flux reaching the walls is difficult to evaluate.

Finally there are only few works on the direct validation of the calculation of radiative transfer obtained by simplified methods by comparison with an exact treatment of this calculation, for thermal plasmas. Indeed it is possible to perform an exact calculation taking into account a very fine description of the spectrum, and using simplified geometry conditions, i.e. imposing a temperature profile (1D calculation) or a temperature field (2D or 3D). For SF 6 Aubrecht and Protasevisch [START_REF] Aubrecht | Radiative Transport of Energy in SF 6 Arc Plasmas Review[END_REF] have shown a comparison between the PCM and exact solution, whereas the recent work of Nordborg and Iordanidis [START_REF] Nordborg | [END_REF] proposed an interesting analysis of various methods.

We have started a general project rather similar to the work of Nordborg and Iordanidis [START_REF] Nordborg | [END_REF] i.e. a comparison between the results obtained by an exact treatment of the radiative transfer and those obtained by simplified methods, in SF 6 plasmas in Local Thermodynamic Equilibrium (LTE). Our main differences with [START_REF] Nordborg | [END_REF] are the following:

-We calculate ourselves the fine spectral variations of the absorption coefficient and we plan to study the influence of various physical mechanisms such as the absorption due to the presence of diatomic and polyatomic molecules. We also can choice more refined description of the spectrum.

-We would like to optimise the number of spectral intervals for averaging the absorption coefficient. Nordborg and Iordanidis proposed a 28 bands model which is too important for plasma modelling.

-We also want to validate the radiative flux and not only the divergence of the flux, because it is very difficult to have simultaneous good values of the radiative flux and of its divergence, using simplified methods. This paper presents a first part of this general project and is devoted to five points: a short presentation of the SF 6 plasma radiation as a function of frequency and temperature with a very fine spectral sampling, accounting for the mechanisms of radiation emission and absorption (part 2); the calculation of NEC showing the role of line overlapping that was not taken into account in our previous calculations (part 3); the exact calculation of radiative transfer in 1D and 2D, with imposed temperature profiles (part 4); the presentation of the simplified description of the spectrum using MAC (with various definitions of the mean absorption), and some improvements from our previous calculations (part 5); the comparison between the exact solution and the approximate results obtained using the NEC and the MAC methods for simple imposed temperature profiles (part 6).

-Basic data for the radiative spectrum of SF 6 plasmas

We consider that the plasma is in LTE and we use the equilibrium composition required for the radiation calculation. Furthermore, assuming that the plasma state is in LTE, the radiation can be described using either the emission coefficient ε ν , or the absorption coefficient K ν because both are related by the Kirchhoff's law (in fact we use ' ν K taking into account induced emission):

ν ν ν ε ' K B = (1)
where B v is the blackbody intensity. Both line and continuum were considered in absorption coefficient calculation. The details of the studied continuum mechanisms (molecular dissociation and ionisation, Bremsstrahlung, electron-ion recombination, electron attachment) can be found in previous papers [3]. In SF 6 at atmospheric pressure the molecules are almost completely dissociated when the temperature is higher than 5000K. In these conditions the most important contribution to the continuum is due to the radiative recombination and we have calculated the Biberman-Schlüter coefficient using the work of Gongassian and Schlüter [22] based on semi-empirical solution of Schrödinger equation based on a scaled Thomas-Fermi potential. At low temperature (T<5000 K) the role of the molecules becomes important in the continuum part of the spectrum, because of photodissociation and photoionisation of these molecules. Nevertheless, in our calculation we do not take into account the radiation due to the rotational and vibrational bands emitted by the diatomic and polyatomic molecules, considering that in SF 6 arc plasmas the central temperature is much higher than 5000 K and that the radiation emitted in the external regions (where T < 5000 K) is negligible in front of the radiation emitted from the hottest regions. So the warm and cold gases play only a role of radiation absorption.

In thermal plasmas the atomic lines emitted from excited levels of neutral and ionised atoms represent in general the main contribution of the total radiation in a large range of temperature. In the case of SF 6 plasma, in the range of temperature between 300 and 30000K, we take into account the following atomic species S, S + , S ++ , S +++ , F, F + , F ++ and F +++ . The atomic structures used are taken from NIST data tables [23]. We consider all the lines tabulated by the NIST [24] and by the calculation of Kurucs and Peytermann [START_REF] Smith | Atomic spectral line database from CD-ROM 23 of R. L. Kurucz[END_REF]; the oscillator strengths (or Einstein coefficients) of the lines are first chosen from NIST compilation [24] and then completed by the values of Kurucs and Peytermann. For SF 6 we consider 8000 lines. One difficult part of the study concerns the determination of line profile.

In previous works we considered only the broadening mechanisms (line shifting was not needed because each line was studied separately for a given temperature), but in this development, line shift must be taken into account. Basic data for line shift and broadening were taken from [START_REF] Traving | Plasma diagnostics Ed Lochte-Holtgreven North[END_REF][START_REF] Griem | Spectral line broadening by Plasmas[END_REF]. The total line profile is then given by:

dy y a y a p T P D D ∫ ∞ ∞ - Δ - + ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ Δ - Δ = )² ( exp ) , , ( 2 ν ν π π ν ν (2)
where: Total absorption coefficient is clearly due to the continuum contribution (with several thresholds corresponding to the electronic structures of the components i.e. atoms and ions of fluorine and sulphur) and to the superposed lines. One of the numerical problems is the optimization of the spectral step to have a good description of the total spectrum. Decreasing this step tends to improve the accuracy but tends to exaggeratedly increase the number of spectral points.

2 L a ω = δ ν ν ν - - = Δ
The choice of the spectral step has been optimized following the following procedure:

-the spectral step Δν, is not constant over all the spectrum. We have first divided the total spectrum in five intervals considering the wavelengths of the resonance lines. Then we have calculated for a given Δν for each interval (ν 1 -ν 2 ) the total emission corresponding to:

ν ν ν ν ν d K B I ' 2 1 ∫ = (3)
-this emission can also be calculated by integrating the continuum and summing the intensities of the lines located in this spectral interval.

∑ ∫ + = l c I d K B I ν ν ν ν ' 2 1 (4)
-the integration (3) is computed for various values of the integration step Δν (with decreasing values of Δν) and compared to the value of expression (4) independent of the integration step.

We consider that the integration step is correct when the relative difference between both values of the total emission is less than 1%. For all the temperature range between 5000 and 30000K at atmospheric pressure (at higher pressure the lines are more broadened than at atmospheric pressure). Tables 1 and2 With this number of points the spectral resolution is fine enough to well describe the line profile and the possible line overlapping as illustrated in figure 2. Finally we would like to focus on one possible weakness of the spectral description of SF 6 plasma. We will see that the absorption of radiation can be important on the edges of the plasma because of the presence of molecules that may absorb radiation by photodissociation and photoionisation.

We have plotted in figure 3 the photo-absorption cross sections of SF 6 , SF 4 and F 2 , given by literature and used in our calculation, versus the wavelength of the absorbed photon. For SF 6 we have made a compilation of 4 papers [START_REF] Hitchcock | [END_REF][29][30][31] for representing the total cross section including dissociation and ionisation. The first dissociation threshold of SF 6 (

F SF SF + → 5 6
) corresponds to 3.995 eV (310.35 nm) but we can see from figure 3 that this dissociation phenomenon is not efficient for photoabsorption near the threshold. The same kind of remark can be noted for SF 4 (data from [31]) for which the first dissociation threshold

( F SF SF + → 3 4
) occurs at 3.477 eV (356.59 nm) and the second threshold (

F SF SF + → )
at 4.79 eV (258.79 nm) and for F 2 [32][START_REF] Okabe | The photochemistry of small molecules[END_REF] for which the dissociation threshold occurs at 1.602 eV (773.41 nm).

So the literature data cover only a part of the spectrum and should be extrapolated. The extrapolation of SF 4 data to the lower wavelength (i.e. higher energy) is not problematic: a constant value around 10 -17 cm 2 can be chosen. The extrapolation to the higher wavelengths is much more problematic. It is not correct to consider that the cross section becomes nil after the last point. Considering the dissociation energies of these 3 molecules we can perform an extrapolation but the shape of the curve is not known and the choice of this extrapolation may have a rather important influence on the radiation absorption at high pressure (p>10 bars), the condition which is usual in high voltage circuit-breaker. Apart from this aspect, we have no available data on photoabsorption cross section for other molecules as SF, S 2 and SF 2 and the presence of these molecules has probably an impact on the radiative transfer.

-Net Emission Coefficient: influence of line overlapping

As it has been written in the introduction, the NEC corresponds to an acceptable estimation of the radiation emitted from the hottest regions of an arc. This method, first proposed by Lowke [1] allows calculating the divergence of the flux at the centre of an isothermal spherical plasma:

q N ∇ = π ε 4 1 (5)
such as:

ν ε ν ν ν d R K K B P N ) ' exp( ' 0 - = ∫ ∞ (6)
where R p is the radius of the sphere. NEC incorporates directly self-absorption. In previous calculations made by our research group, the lines were treated independently, i.e. neglecting line overlapping. For each line an escape factor Λ was introduced to represent the ratio between a real emitted power and the emitted power without absorption. This factor is equal to 1 when the line is not absorbed and tends to zero when self-absorption increases. This assumption of negligible line overlapping simplifies the calculation because for each line the escape factor is previously computed as a function of the number density of the excited atomic level (deduced from the equilibrium composition) and of the broadening parameters, following the method of Drawin and Emard [START_REF] Drawin | [END_REF]. So the integration of the NEC along the frequency (Eq. 6) is reduced to the integration of the continuum absorption coefficient and the summation of the contributions of the lines. This assumption was made in order to reduce CPU time, knowing that the corresponding error on the NEC at atmospheric pressure was acceptable. For given values of pressure and temperature, it is important to note that the escape factor depends on the plasma size, and thus the NEC depends on R p . This parameter that is the radius of the plasma sphere in the initial calculation, should be estimated in practical calculations (in general the plasma is not spherical) as the conduction radius of the plasma (in fact the distance from the axis where the temperature represents around 70-80% of the axis temperature).

Using the very fine description of the spectrum for each temperature, as illustrated in figure 1, with 300 000 frequency points, it is possible to calculate the NEC (Eq. 6) by direct integration. The results of the NEC function of temperature are plotted in figure 4 for two pressure values. In the fictitious case R p =0 (in fact without any self-absorption) there is no difference between the values computed with the escape factors or without (in this latter case, that means integration over around 300 000 points in wavelength or frequency). This is normal and shows that the calculation with the total integration is correct. Then, figure 4 shows that about 90% of radiation is absorbed over the first millimetre of plasma, which is a well known result for the majority of thermal plasmas.

What is interesting in our results is the difference between the two kinds of results (with and without escape factors). We can see in figure 4 that the influence of line overlapping, weak at atmospheric pressure, is not totally negligible at higher pressure and tends to increase when temperature and pressure increase. This is due to the fact that several lines are very close (specifically triplets) and that the increase of temperature and pressure tends to broaden the lines (Stark and pressure effects) and then to favour the overlapping. The influence of overlapping may reach 40% of the NEC at a pressure of 8 bars.

-Exact calculation

Radiative transport equation is given by:

) ( ' ) , ( . ˆν ν ν ν I B K s r I s - = ∇ (7)
where I ν is the spectral radiative intensity; r, s ˆ the position and the direction, respectively.

Following one direction Ox, we obtain an ordinary differential equation:

( ) ν ν ν ν I B K dx dI - = ' ( 8
)
with the solution :

dl e l B l K e I x I x l x d K x dl l ∫ + ∫ = - - ∫ η η ν ν κ ν ν ν ν ) ( ' 0 ) ( ) ( ) ( ' ) 0 ( ) ( 0 ' (9)
In 1D configuration, knowing the temperature profile, the equation of radiative transfer can be numerically solved following Eq. ( 9). With a high number of mesh points, the solution can be considered as exact. In our calculation, we have imposed two kinds of temperature profiles, shown in figure 5, corresponding to the profiles along a radius. In 1D, we perform the calculation along a profile corresponding to a diameter (maximum temperature on the middle). In 2D, we calculate the radiative transfer in a plane with a circular geometry (section of cylinder plasmas, maximum on the centre) using discrete ordinates. Nordborg and Iordanidis [START_REF] Nordborg | [END_REF] showed that this method gives very similar solutions as the direct integration given in their paper, provided the number of directions is high enough. we present it only on half a segment (1D) of along the radius (2D). Among the results we must notice that for the considered case the divergence becomes negative for a radial position higher than 2 mm which corresponds approximately to T<15 000 K (see figure 5) and that a strong absorption occurs on the edge of the plasma when the temperature becomes lower than 3000K. The latter effect is mainly due to the absorption of radiation by the molecules which correspond to photodissociation and photoionisation of these molecules. From these preliminary tests, the general results presented in the last part were obtained using 200 spatial points in one direction, and using 36 directions for 2D computation.

-Mean Absorption Coefficients

The use of the fine spectral description (100 000 points or more for each temperature) is not possible in thermal plasma modelling where the calculation of temperature field is performed iteratively for each spatial point. On the other hand, the NEC is unable to calculate the absorption in the cold regions or to estimate the radiative flux. In these conditions, a simplification of the spectral description is needed to perform the iterative calculation of the radiative and temperature fields. The simplest approximation is to consider that the total spectrum can be separated in a limited number of intervals and that the plasma can be considered as a grey body in each interval, i.e. for a given temperature the absorption coefficient is constant within each interval. Then these variations of the mean absorption coefficients versus temperature can be used in the general modelling by means of classical methods of radiative transport (P1 method for example). However it is to be noted that for each interval a partial derivative conservation equation can be solved and in order to limit the cpu time the number of intervals should be very limited (less than 10 for 3D models for example).

The procedure is the following:

-from the evolutions of the absorption coefficient K' versus frequency it is rather easy to select between 5 and 8 intervals where K' varies slowly within each interval for all the temperature range.

-For each interval we may calculate a MAC for the continuum. We may choose 3 kinds of definition for the MAC and we will discuss on them just below.

-Finally for each interval we must add the contribution of the lines.

We present in figure 8, the variations of the continuum absorption coefficient versus temperature at atmospheric pressure for various temperature values. As explained in previous papers (see [3] for example) the sudden variations of K' are due to the electronic structure of the atoms (neutral and ions for T>5000 K) or of the molecules (T<3000 K). We have also plotted in this figure the arbitrary limits of 7 spectral intervals (vertical black dotted lines), considering that the absorption coefficient is rather constant for a given temperature (grey body approximation). The total range of studied frequency is between 2.10 14 Hz and 10 16 Hz (wavelength range between 30 nm and 1.5 μm). In fact, for a frequency interval between ν i and ν i+1 we may define three kinds of average. The Planck average is defined as:

ν ν ν ν ν d K B W K i i i i ' 1 ' 1 ∫ + = 〉 〈 (10)
with

ν ν ν ν d B W i i i ∫ + = 1 (11)
The Rosseland average is defined as: So in a general modelling we may first use the NEC (very fast) to have a good estimation of the temperature field and then engage a calculation of radiative transfer (P1 or discrete ordinates for example) using our values of MAC, without destabilizing the numerical code. A comparison of the values of the mean absorption coefficients for 7 bands in pure SF 6 at atmospheric pressure versus temperature, is given in figure 9 for R p =2 mm. For bands 1 to 3 (high energy) the curves are superimposed at low temperature and the differences occur at high temperature (over around 15000K) due to the presence of ionic lines that have a very weak contribution in the Rosseland averaging. For band 4 the difference is very marked and is due to intense resonance lines emitted by excited neutral atoms (we have verified that the Planck averaging of the continuum is superimposed to the Rosseland averaging for the total radiation). Finally for bands 5-7 the differences between the two kinds of results (Planck and Rosseland) are mainly, but not totally, due to the lines.

-Comparison of methods and discussions

This general comparison will be made using the temperature profiles defined in figure 5. As the profile is considered symmetrical, the divergence of the flux is also symmetrical and will be presented between x= 0 and x=R; R is the profile radius equal to 5 mm in case 1 and to 2 mm in case 2. The calculation has been performed in 1D and 2D as explained previously. The flux is not symmetrical because it depends on the direction and will be given only in 1D geometry between x=-R and x=R, corresponding to the diameter (the calculations were also performed for all the cords of the cylindrical section but will not be presented here because they gave similar results).

Net emission coefficient

Figure 10 presents the comparison of the divergence of the radiative flux obtained either by the net emission coefficient or by the exact method, for the case 1 (results very similar to case

2) and for a pressure of 1 bar. The NEC depends also on the parameter R p which is the radius of the isothermal sphere used for computing the NEC, and it corresponds in practical case to an equivalent plasma radius. In our case we show the results for two R p values, 1 and 2 mm. The results are classical. The NEC gives a rather satisfactory but not very good agreement in the hottest regions, whereas by definition the NEC is unable to compute a negative value of the divergence and thus is unable to treat the absorption of radiation in tedious and cold regions. It is also interesting to note another classical result showing the difference between 1D and 2D exact calculations: the 1D calculation performed following a diameter, always shows a more pronounced absorption in the outer regions than the 2D condition that averages this absorption along all the directions.

Mean Absorption coefficient: first results

Using the seven spectral intervals defined in figure 8 and the values of the Planck and Rosseland averaged mean absorption coefficients plotted in figure 9, we may calculate the radiative transfer for the given temperature profiles (defined in figure 5) and compare the results with the exact solutions in 1D and 2D conditions. This comparison on the divergence of the radiative flux (in 1D and 2D) is presented in figures 11a-b. From these figures it is clear that the use of the Rosseland averaging is not satisfactory because of a big and unacceptable underestimation of the divergence of the flux in the hottest regions. On the contrary we observe a qualitative agreement between the exact values (with almost 300 000 spectral points) and the approximated values using the Planck averaging (grey body with 7 spectral intervals): the differences are acceptable in the hottest regions whereas both methods show absorption of radiation in tedious and cold regions. Nevertheless we can see in figure 12 that in 1D condition the radiative flux near the wall (i.e. for the radial position equals to 5mm in case 1) depends rather strongly on the choice of R p when using the Planck averaging method and can be different from the exact solution. Using the value of R p giving the best agreement on the flux divergence near the centre (R p =1 mm, see figure 12) the approximated value of the flux at the edge is twice the exact value. This result (overestimation of the flux) means also that the total absorption of radiation is underestimated by the approximated method based on Planck MAC which may have a strong influence on the pressure rise within a SF 6 HVCB. We also may note from figure 12 the errors induced by the use of the Rosseland MAC.

-Conclusion

The improvement of the calculation of radiative transfer in thermal plasmas established in molecular gases (here SF 6 but the method is general), can be made using two correlated ways.

The first one that is physically complex consists in well treating the basic phenomena important in thermal plasma radiative transfer. The first kind of big difficulties, with rather high uncertainty, corresponds to the determination of atomic line profiles, in particular for resonance lines, responsible for a large part of the emission in spite of a strong selfabsorption. The second kind of difficulties concerns the photodissociation and photoionisation of molecules, responsible for a large part of radiation absorption. In the case of SF 6 we have two main problems that are not yet solved and that should be studied in the future: the role of the presence of diatomic molecules S 2 and SF that is not taken into account; the knowledge of the absorption cross sections in all the UV part of the spectrum because data are limited to the VUV range in general. This unknown part may play a rather important role in the radiation absorption at high pressure (in HVCB pressure may reach several tens of bars).

The second way to improve the radiative transfer is to optimise the spectral simplified description. Nordborg and Iordanidis [START_REF] Nordborg | [END_REF] presented a good approach and here we have proposed a method based on a lower number of spectral intervals that give satisfactory results using the Planck and the natural averaging for the mean absorption coefficients, whereas the use of the Rosseland averaging is not suitable. Other efforts should be developed to check the method in more general conditions (3D and higher pressures) and using other usual methods (P1 and discrete ordinates) based on the simplified description of the radiation.

ω

  represent the Doppler and Lorentz widths respectively, δ is the line shift.Expression (2) corresponds to a Voigt profile, convolution of Gaussian and Lorentzian profiles due to the broadening phenomena: Doppler, resonance, Van Der Waals and Stark effects. For an SF 6 plasma an example of spectral distribution of the absorption coefficient is given in figure1for two high values of temperature.

Figure 1 :

 1 Figure 1: Absorption coefficient of SF 6 at atmospheric pressure for temperatures of 10 000 K and 20 000 K respectively.

Figure 2 :

 2 Figure 2: Example of line overlapping with the fine spectral description (dashed line corresponds to the continuum).

Figure 3 :

 3 Figure 3: Photoabsorption cross section for 3 molecules

Figure 4 :

 4 Figure 4: Influence of line overlapping on the net emission coefficient for SF 6 plasma, for 2 pressure values

Figure 5 :

 5 Figure 5: Imposed temperature profiles In 1D, the accuracy depends on the number of spatial points N, and in 2D it depends also on the number of directions n. Figures 6 and 7 show the influence of these parameters on the divergence of the flux. For geometrical reason the divergence of the flux is symmetrical and

Figure 6 :Figure 7 :

 67 Figure 6: Calculation in 1D condition, of the divergence of the flux: influence of N, the number of spatial points

1 Figure 8 :

 18 Figure 8: Evolutions with frequency of the continuum absorption coefficient for various temperature values. The limits of spectral intervals for the grey body approximation are plotted in vertical black dotted lines and the numbering of the intervals (called also bands) are indicated in the figure (Band 1 corresponds lo higher frequency or photon energy)

Figure 9 :

 9 Figure 9: Evolutions of the mean absorption coefficient (the Planck and Rosseland averages) for the 7 bands, versus temperature and for R p = 2 mm for the Planck averaging

Figure 10 :

 10 Figure 10: Comparison of the NEC with the exact solution for the radiative flux divergence

1 )Figure 11 :

 111 Figure 11: Comparison of the values of the divergence of the radiative flux between the exact solution and the approximated method based on 7 bands with the mean absorption coefficients. a) 1D calculation; b) 2D calculation

Figure 15 :

 15 Figure 15: Evolutions of the absorption coefficient in Band 4 and comparison with the Planck, Rosseland and normal averaged values Another improvement can be made if we observe on figure 8 the variations of the absorption coefficient in band 5 at 2000K. The partial maximum is due to the presence of SF 4 and was not taken into account in the previous choice of the intervals. Following the results we can group bands 2 and 3 and split into two new bands the old band 5. With this new choice of spectral intervals and replacing the Planck average by the normal average for T<6000 K, we have calculated again the radiative transfer using the new 7 bands and one characteristic result is presented in figure 16 equivalent to figure 12. Obviously the final radiative flux at d = 10 mm obtained by the approximated method is now much better than the result obtained by the previous choice. The same kind of improvement has been noticed with the second temperature profile (case 2).

Figure 16 :

 16 Figure 16: Evolution of the radiative flux in 1D condition; comparison between exact solution and approximated values obtained using 7 new bands with the Planck mean absorption coefficients for T>6000K and normal mean absorption coefficients for T<6000 K.

Table 1 :

 1 give the main results, i.e. the absolute values of the emission computed by (3) and compared to the values computed by (4) (table1) and the relative difference (table2), for various numbers of spectral integration points, and for 6 values of temperature. The final choice corresponds to 286802 points. Values of the intensities given by relations (4) (upper line) and(3). The first column indicates the number of spectral points used for computing ralation(3)

	Nb of points	5000	10000	15000	20000	25000	30000
	Relation (4) 5.7242D+05 1.1431D+09 5.3413D+09 2.1459D+10 9.2400D+10 2.4105D+11
	160 802	5.6482D+05 1.1205D+09 5.2453D+09 2.1738D+10 -	-
	286 802	5.7084D+05 1.1442D+09 5.3792D+09 2.1558D+10 9.2417D+10 2.4211D+11
	538 801	5.7240D+05 1.1457D+09 5.3649D+09 2.1550D+10 9.2787D+10 2.4256D+11
	Nb of points	5000	10000	15000	20000	25000	30000
	160 802	-1.32D-02	-1.9D-02	1.79D-02	1.30D-02	-	-
	286 802	-2.76D-03	9.62D-04	7.09D-03	4.61D-03	1.83D-04	4.39D-03
	538 801	-3.49D-05	2.27D-03	4.41D-03	4.24D-03	4.18D-03	6.26D-03

Table 2 :

 2 Relative difference between the intensities computed by relations (3) and(4)

The natural average is:

Of course when the absorption coefficient is really constant within a spectral range, all the definitions give the same mean value. The problem occurs when ν ' K is not constant and in particular when strong lines appear in the spectral range. It is well known that the Planck average is a good approximation when ν ' K is rather low, precisely when the optical depth is low so that the intensity as defined in relation ( 9) is almost reduced to the local emission

. The Rosseland averaging coefficient is in general more precise for high values of ν ' K , and presents the great advantage that the central part of the strong lines (that is totally reabsorbed in the medium) has no influence on the mean value because of the definition in 1/K. Finally the natural average can be useful if we consider only the role of absorption (special conditions such as the local emission is much lower than the local absorption).

The previous calculation of radiation (for example through the NEC) have shown that the lines play in important role in thermal plasma radiation and it is well known that in these conditions the assumption of grey body intervals can lead to severe and unacceptable errors using the Planck averaging absorption coefficient [START_REF] Modest | Radiative Heat Transfer[END_REF][START_REF] Nordborg | [END_REF]. This is due to the fact that the most intense lines are strongly self-absorbed which is not taken into account if we sum directly the contribution of the lines in the interval [ν i -ν i+1 ]. In our previous works [START_REF] Gleizes | XXIII ICPIG Int. Conf. Phen. Ionized Gases Toulouse[END_REF][START_REF] Erraki | Etude du transfert radiatif dans les plasmas thermiques : application au SF6 et au mélange argon-fer[END_REF] we suggested to use the line escape factors, Λ l , introduced in part 3 or the NEC. So in relation [START_REF] Aubrecht | [END_REF] the contribution to the integral, of a line centred at 0 ν ν = is written as:

The escape factor depends on the temperature and on the size of the plasma characterized by the radius R p of the isothermal sphere used to compute the NEC. Thus our total MAC for a given interval depends on temperature and on R p . The best interest in this approach is that for the hottest regions we obtain a good consistency between the NEC and the MAC: All the previous results were obtained using one temperature profile (case 1 in figure 5)

corresponding to a maximum temperature of 20000K). Using the other temperature profile (case 2) the results are similar as it can be seen in figure 13, showing the divergence of the flux computed in 2D geometry..

Improvements on the description by spectral intervals

In order to evaluate the contribution of the various spectral intervals to the total radiative flux we have plotted in figures 14a-b the contributions of the seven bands in the Planck approximation (figure 14a) and the same corresponding contributions obtained by the exact solution (figure 14b). Intervals (or bands) 1-3 have a large contribution to the radiation, but the resulting flux for these high energy intervals at r = 5mm is nil, because of the strong absorption in the edges. Band 6 contributes strongly to the final flux and is well treated with the approximated method. On the contrary the approximated method underestimates the absorption of radiation in bands 4 and 5. So the improvement of the calculation of the radiative flux on the plasma edges needs a better treatment for these two bands. We have explained that the use of the Planck averaged absorption coefficient allows to have a rather good estimation of the net emission. But if we consider results presented in figures 11

we can see that the divergence of the flux is negative (i.e local absorption is greater than local emission) in regions where the temperature is lower than or of the order of 10000K, for case 1 of temperature profile. In any case when the central temperature is higher than 10000K, which is always the case in SF 6 arc, the regions where T<6000K are mainly absorbing regions. In this case is the Planck averaging useful? An answer to this question can be found analysing figure 15 representing the evolution of the absorption coefficient for Band 4 (zoom of figure 8 for this band or interval) at P = 1 bar and T = 300 K, and indicating three values of mean absorption previously defined by relations [START_REF] Aubrecht | [END_REF][START_REF] Modest | Radiative Heat Transfer[END_REF][START_REF] Gleizes | XXIII ICPIG Int. Conf. Phen. Ionized Gases Toulouse[END_REF].


Jeremie Guedj
email: jeremie.guedj@gmail.com

Caroline Bazzoli

Avidan U Neuman

France Mentré

Design evaluation and optimization for models of hepatitis C viral dynamics

Keywords: Design Evaluation, Design Optimization, Hepatitis C, Non-linear mixed effect models, Ordinary Differential Equations, Viral kinetics

Mathematical modeling of hepatitis C viral (HCV) kinetics is widely used for understanding viral pathogenesis and predicting treatment outcome. The standard model is based on a system of five non-linear ordinary differential equations (ODE) that describe both viral kinetics and changes in drug concentration after treatment initiation. In such complex models parameter estimation is challenging and requires frequent sampling measurements on each individual. By borrowing information between study subjects, non-linear mixed effect models can deal with sparser sampling from each individual. However, the search for optimal designs in this context has been limited by the numerical difficulty of evaluating the Fisher information matrix (FIM). Using the software PFIM, we show that a linearization of the statistical model avoids most of the computational burden, while providing a good approximation to the FIM. We then compare the expected precision of the parameters that can be expected using five study designs from the literature. We illustrate the usefulness of rationalizing data sampling by showing that, for a given level of precision, optimal design could reduce the total number of measurements by up 50%. Our approach can be used by a statistician or a clinician aiming at designing an HCV viral kinetics study.

Introduction

Chronic hepatitis C virus (HCV) infection is one of the most common causes of chronic liver disease, with as many as 170 million people infected worldwide [START_REF] Shepard | Global epidemiology of hepatitis C virus infection[END_REF]. HCV is considered as eradicated when a sustained virological response (SVR), defined by the persistent absence of serum HCV RNA for 6 months after therapy, is achieved. The current standard of care combines weekly injections of pegylated interferon (Peg-IFN) and daily oral ribavirin and leads to an SVR rate of 50% in HCV genotype 1 patients (2), which is the largest group of patients in western countries.

By mathematically describing early HCV RNA (viral load) decay after initiation of IFN-based antiviral therapy, crucial parameters of the in vivo viral kinetics have been estimated, such as the rates of production and clearance of free virus, and the rate of loss of infected cells [START_REF] Neumann | Hepatitis C viral dynamics in vivo and the antiviral efficacy of interferon-alpha therapy[END_REF]. Furthermore, by suggesting mechanisms of action for IFN and ribavirin, mathematical modeling has provided a means for evaluating and optimizing treatment strategies [START_REF] Mihm | Review article: predicting response in hepatitis C virus therapy[END_REF]. The standard approach is based on five non-linear ordinary differential equations (ODE) that describe the interaction between the virus and its target cells, and the changes in antiviral drug effectiveness with varying levels of Peg-IFN and ribavirin [START_REF] Talal | Pharmacodynamics of PEG-IFN Differentiate HIV/HCV Coinfected Sustained Virological Responders From Nonresponders[END_REF].

A first step when dealing with such a complex model is to ensure that model parameters can be uniquely identified, e.g., that there is a one-to-one map between the model outputs and the parameters.

The identifiability study of viral dynamics models has attracted a great deal of attention from mathematicians (6; 7) and has resulted in important theoretical results, such as the identification of parameters that cannot be uniquely determined from data and the minimal number of measurements required for parameter estimation (8; 9). However, these results only provide the minimal amount of information that needs to be collected in order to identify the parameters under perfect conditions. In reality, more data are needed to obtain reliable estimates of model parameters that circumvent measurements error and/or model misspecification. To address this issue, the notion of practical identifiability has been introduced and refers to the ability to provide reliable parameter estimates [START_REF] Han | Design of population studies of hiv dynamics. Deterministic and stochastic models of AIDS epidemics and HIV infections with intervention[END_REF][START_REF] Guedj | Practical Identifiability of HIV Dynamics Models[END_REF][START_REF] Miao | Modeling and Estimation of kinetics parameters and replicative fitness of HIV-1 from flow-cytometry-based growth competition experiments[END_REF]. When the parameters are estimated using maximum likelihood approach, the practical identifiability boils down to the ability to get "satisfactory" standard errors of the parameter estimates. Hence, the practical identifiability is a relative notion, that will depend on the information available, both qualitatively and quantitatively. In nonlinear contexts, such as viral dynamics models, whether a parameter is practically identifiable may highly depend on the data, the schedule of the sampling measurements and on the number of patients, when the parameters are estimated in a population framework, i.e., using nonlinear mixed effect models [START_REF] Lavielle | Maximum likelihood estimation of long-term hiv dynamic models and antiviral response[END_REF].

In order to study and optimize the practical identifiability, it is natural to evaluate the Fisher Information Matrix (FIM). Indeed the Cramer-Rao bound states that maximum likelihood estimated parameters are asymptotically efficient, with an asymptotic variance given by the inverse of the FIM. Nevertheless, in nonlinear mixed effect models, the FIM has no closed form, and its computation involves solving a multi-dimensional integral. The level of complexity and the computational time burden are still increased by the complexity of the nonlinear ODE model used in viral dynamics models. We have proposed an approach that advantageously uses the structure of ODEs to compute "exactly" (more rigorously, as precisely as wanted) the FIM [START_REF] Guedj | Practical Identifiability of HIV Dynamics Models[END_REF]; nevertheless, this approach involves greedy computations that makes its use complicated for design evaluation and impossible for design optimization. Another approach developed by Mentré et al. is to use a linearization of the statistical model around the expectation of the random effects [START_REF] Mentré | Optimal design in random-effects regression models[END_REF]. This approximation provides analytical values for the FIM elements and thus avoids the computational burden. Relevance of this approximation has been shown in models with single and multiple responses but only for models with analytical solutions [START_REF] Bazzoli | Fisher information matrix for nonlinear mixed effects multiple response models: evaluation and comparison by simulation to results from FO, FOCE and SAEM algorithms[END_REF]. An application to multiple responses with nonlinear ODEs was published but without evaluation of the approximation [START_REF] Bazzoli | Design evaluation and optimisation in multiple response nonlinear mixed effect models: PFIM 3.0. computer methods and programs in[END_REF]. Here, we show the relevance of this approximation in the context of HCV viral dynamics models when two responses (Peg-IFN & HCV RNA) modeled by nonlinear ODEs are observed. Using this approximation we compare the practical identifiability that can be expected using five study designs of the literature, and we propose several optimal designs according to the number of sampling measurements allowed in each patient.

Methods

Viral dynamics modeling

The viral dynamics model considers target cells, T, productively infected cells, I and viral particles, V. Target cells are produced at a rate s and die at a rate d. Cells become infected with de novo infection rate b. After infection, these cells are lost with rate δ. In the absence of treatment, virus is produced by infected cells at a rate p and cleared at a rate c. The time origin, t = 0, corresponds to the time of treatment initiation. We assume that treatment is initiated sufficiently long after the acute infection so that the system is in equilibrium, i.e., that viral production is balanced by viral clearance, and dV dt = dT dt = dI dt = 0. Peg-IFN acts mainly by reducing the rate of virion production per infected cell, with effectiveness ǫ, therefore the viral production under treatment is p(1ǫ), where 0 < ǫ < 1. Consequently, treatment initiation leads to an initial rapid decline of viral load with rate c and magnitude log 10 (1ǫ). After this initial drop, the viral load enters a new slower mode of decline that reflects the progressive loss rate of infected cells, δ. However, with weekly administration of Peg-IFN, drug levels can drop rapidly and become too low to maintain a high degree of suppression of viral load. Consequently, a new round of cell infection may occur and give rise to increases in viral load towards the end of the dosing interval.

These rebounds can be captured by assuming an E-max relationship between drug effectiveness, ǫ(t), and the drug concentration, C(t), such that ǫ(t) =

C(t) n C(t) n +EC n 50
, where EC 50 is the drug concentration in blood at which the drug is 50% effective, and n is the Hill coefficient, a parameter that determines how steeply the treatment effectiveness changes with variations in drug concentration [START_REF] Talal | Pharmacodynamics of PEG-IFN Differentiate HIV/HCV Coinfected Sustained Virological Responders From Nonresponders[END_REF]. The changes in Peg-IFN concentrations, C(t), can be described by a first-order absorption and elimination model [START_REF] Talal | Pharmacodynamics of PEG-IFN Differentiate HIV/HCV Coinfected Sustained Virological Responders From Nonresponders[END_REF].

In the long-run, weekly viral load levels in good responders will decline linearly with rate δ (Figure 1, A andB).

Although ribavirin showed only little antiviral efficacy in monotherapy [START_REF] Pawlotsky | Antiviral action of ribavirin in chronic Hepatitis C[END_REF], the combination of IFN or Peg-IFN with ribavirin dramatically improves SVR rates, particularly in patients that respond poorly to interferon-based treatment [START_REF] Layden-Almer | Viral dynamics and response differences in HCV-infected African American and white patients treated with IFN and ribavirin[END_REF]. The mechanisms of action of ribavirin are poorly understood, although one predominant mechanism is thought to be mutagenesis, whereby the proportion of non-infectious virus released is increased [START_REF] Dixit | Modelling how ribavirin improves interferon response rates in hepatitis C virus infection[END_REF]. This can be modeled by assuming that ribavirin reduces the de novo infection rate with effectiveness η (e.g., infection under treatment becomes (1η)b, where 0 < η < 1) [START_REF] Herrmann | Hepatitis C virus kinetics[END_REF]. Since ribavirin is administrated on a daily basis, its effectiveness can be considered constant, although a slow accumulation of ribavirin has been observed in plasma [START_REF] Dixit | Modelling how ribavirin improves interferon response rates in hepatitis C virus infection[END_REF]. Furthermore, the HCV viral dynamics model with Peg-IFN and ribavirin treatment is based on a set of five nonlinear ordinary differential equations (more details on the modeling aspects can be found in [START_REF] Guedj | A perspective on modeling hepatitis C virus infection[END_REF]):

                           dX dt = -k a X + e(t) dA dt = k a X -k e A dT dt = s -(1 -η)bV T -dT dI dt = (1 -η)bV T -δI dV dt = 1 - C(t) n C(t) n +EC n 50 pI -cV (1)
for t = [0, 7, 14, 21, ...], e(t) = F D where D is the drug dose and F is the bioavailability (0 ≤ F ≤ 1). X(t) and A(t) are the amount of drug at the injection site and in the blood, respectively. The concentration of drug in the blood is given by C(t) = A(t)/V d where V d is the drug's volume of distribution.

Since only the viral load and the Peg-IFN levels are measured, the parameters p, s, d, V d , F cannot be identified (8; 9). Therefore, without loss of generality, we fix these parameters to the following values: p = 100, s = 20000, d = 0.001 and F = 1. Although η and b are theoretically identifiable, the lack of information on the dynamics of T and I make them poorly identifiable [START_REF] Wu | Parameter Identifiability and Estimation of HIV/AIDS Dynamic Models[END_REF]. Hence, unless it is specified otherwise, we will consider b = 10 -7 and η = 0.

Statistical Model

To ensure that parameters take on positive values, let us define [log(EC 50 ), log(n), log(δ), log(c), log(k a ), log(k e ), log( the vector of dimension p = 7 of the log-transformed parameters to be estimated. The between-patients variability is taken into account by considering that each individual parameter value, β i , is composed of a fixed part (the fixed effects parameter β) and a random part b i , assumed to be Gaussian with a diagonal variance-covariance matrix Ω = diag(ω 2 1 , ..., ω 2 p ):

     β i = β + b i b i ∼ N (0, Ω) (2)
The observations are of concentrations of viral load (HCV RNA) and Peg-IFN in the blood plasma.

Both markers are observed at the same time points. We denote by

d i = {t ij } j=1,...,
     Y ij1 = log 10 (V ((t ij , β i )) + ǫ ij1 , j = 1, ..., n i Y ij2 = C(t ij , β i ) + ǫ ij2 , j = 1, ..., n i (3)
As usual, for all N patients, ǫ i |b i , i = 1, .., N are assumed to be independent from one individual to the other and from one response to the other. Lastly, for each individual, ǫ i and b i are also independent.

Fisher information matrix

Let ψ be the vector of parameters to be estimated, such as ψ T = β T , ω 2 1 , ..., ω 2 p , σ T , and let λ be the vector of variance terms, such as

λ T = ω 2 1 , ..., ω 2 p , σ T , such that ψ T = β T , λ T . Let l i (ψ; Y i ) be the observed likelihood for individual i: l i (ψ; Y i ) = R p j=1,...,ni 1 σ 1 √ 2π exp - 1 2 Y ij1 -log 10 (V (t ij , β i )) σ 1 2 × 1 σ 2 √ 2π exp - 1 2 Y ij2 -C(t ij , β i ) σ 2 2 dp(b i ). (4)
The log-likelihood is for patient i is defined as

L i (ψ; Y i ) = log(l i (ψ; Y i )),
M F (ψ; d i ) = E ψ - ∂ 2 L i (ψ; Y i ) ∂ψ∂ψ T , (5)
and the FIM for the whole sample is defined by

M F (ψ; D) = E ψ - ∂ 2 L(ψ; Y ) ∂ψ∂ψ T . (6)
By independence between patients, the FIM for the whole sample is simply

M F (ψ; D) = i=1,...,N M F (ψ; d i ).
If all patients have the same design d, then the FIM for the whole sample is

M F (ψ; D) = N × M F (ψ; d).

The PFIM software

Approximation of the FIM.

By making a first-order approximation around the expectation of the random effects in the likelihood expression (Eq. ( 4)) and using the assumptions made in section 2.3, Mentré et al. showed that an analytical approximation of the FIM can be obtained [START_REF] Mentré | Optimal design in random-effects regression models[END_REF]. The approximated FIM is block diagonal with a block for the fixed effects parameters and a block for the random effects parameters. A major advantage to this approximation is avoiding the multiple integrals involved in the exact computation of the FIM (Eq. ( 6)). This approximation of the FIM has been proposed in some free software, including PFIM, PopED and WinPopt. A comparison of the different codes can be found at http://www.pagemeeting.org/?abstract=1179. Here we focus on PFIM, which has the advantages of being written in R, and of using the Federov-Wynn algorithm for design optimization. This software has recently been expanded to deal with multiple response models (PFIM 3.0), including those defined as ODE models (15; 16). The software is available at www.pfim.biostat.fr and the code used in this article is available from the authors upon request.

Comparisons of study designs.

The comparison of the designs is based on their efficiencies. The efficiency criterion, φ(D), is defined by the determinant of M -1 F (ψ; D) normalized by the number of parameters to estimate, namely

φ(D) = |M F (ψ; D)| 1 dim(ψ)
. Then the efficiency of a design D 1 may be compared to that of D 2 by computing the ratio of the efficiency criteria, Eff(D

1 , D 2 ), i.e., Eff(D 1 , D 2 ) = φ(D 1 )/φ(D 2 ). If Eff(D 1 , D 2 ) > 1, then D 1 is more efficient than D 2 and Eff(D 1 , D 2 ) is the mean gain in precision using D 1 instead of D 2 .

Design Optimization.

For a given value ψ of the parameters, a population design, D, is D-optimal, denoted by D D , if it minimizes the determinant of the inverse of the FIM:

D D = Argmin D 1 |M F (ψ; D)| (7)
When the total number of sampling measurements is fixed and there is only a finite set of admissible times, the Fedorov-Wynn algorithm (22; 23) implemented in PFIM can be used to optimize designs based on the criterion φ(D). This algorithm provides an optimal design described by a number of elementary designs with the corresponding proportion of subjects in each elementary design.

Design comparisons

Designs studied and parameter set-up

We study the practical identifiability that can be expected with five designs used in seminal clinical trials aimed at describing viral dynamics with Peg-IFN (5; 24-27) (Table 1). These designs assume that both viral load and Peg-IFN concentrations are measured at each sampling time. Each design uses a different strategy. For example, in D 1 , there are only 8 measurements per patient, spread over four weeks, with sample measurements one day after each injection, close to viral load troughs (and drug peaks). On the other hand, D 2 -D 5 focus their samplings on the first week of treatment and on injection times that correspond to viral load peaks and drug troughs.

Moreover, in order to study how parameter precisions depend on the values of the fixed effects, we studied the SE obtained for four possible parameter set-ups A-D for the fixed effects (Table 2). These four set-ups are meant to cover the different median values that could occur in a clinical study of Genotype 1 . The dose is the same for all the patients (180 µg). All the random effects have a standard deviation ω = 0.5 and the error measurements are such that σ 1 = σ 2 = 0.2.

Patterns of viral kinetics

EC 50 n δ c k a k e V d (µgL -1 ) (d -1 ) (d -1 ) (d -1 ) (d -1 ) (L)
A. Good responders with low inter-dose variation 0.12 2 0. Since drug concentration per se was not found to be associated with long-term response (5; 29), good and poor responders differ primarily in their pharmacodynamics (EC 50 , n) and in their loss rate of infected cells, δ.

Comparison between the approximated standard errors given by PFIM and the empirical standard errors

We first investigate whether the approximate standard errors SE P F IM are a good approximation of the empirical standard errors SE EM P . To address this issue, we simulated, for each design D The fact that the linear approximation seems to work less well for EC 50 , especially for the poorest design D 1 , was already observed in another context [START_REF] Bazzoli | Fisher information matrix for nonlinear mixed effects multiple response models: evaluation and comparison by simulation to results from FO, FOCE and SAEM algorithms[END_REF]. Indeed, EC 50 bridges the pharmacokinetics and the pharmacodynamics and, hence, is associated with a high degree of nonlinearity [START_REF] Bates | Relative curvature measures of nonlinearity[END_REF], which has been shown to hamper the quality of the first-order approximation (32; 33). Of note, ∆(SE) > 0.1 for c in D 1 results from its high degree of correlation with EC 50 in this design, where SE for these parameters are the largest.

Design evaluation

We compared the SE for the designs D 1 -D 5 using parameter set-ups A-D. Using the δ-method, the SEs of the estimators of the log-transformed fixed effects β are asymptotically equal to the relative error RE that can be expected for the parameters on their natural scale. Hence the results given in Table 3 can be read in the following way (first line): with design D 1 , a patient sample size N = 30, and parameter set-up A, the expected relative error RE for δ is RE(δ) = SE(δ)/δ ≃ SE(log(δ)) = 0.10.

Hence, the confidence interval (at the level of 95%) that can be expected for δ is CI 95% (δ) = [0.16; 0.24].

Not surprisingly, the RE tend to improve as the number of sampling points increases. However since the model is nonlinear, the relevance of a sampling point depends on the sensitivity of the biological and statistical model at this particular timepoint and for this particular parameter set-up. For instance, the most intensive design, D 5 , improves the RE of {EC 50 , n, c} compared to D 1 -D 4 , but degrades the RE of δ for good responders (A and B). Indeed, D 5 concentrates most of the sampling points in the first week of treatment with no sampling point at week 3 or week 4. Given the mathematical properties of the biological model, the large number of early sampling points allows for a good estimation of {EC 50 , n, c}.

However, since the effects of the loss rate of infected cells, δ, are mainly seen in the long-term viral decline, D 5 is not the most appropriate for estimating δ. This result is less pronounced for the parameter set-ups C and D, for which patients respond poorly to therapy and have a large viral rebound from days 2-3. Indeed, the kinetics of the rebound depend on the infected cell turnover and, hence, substantial information on δ can be gained in this case within the first week of treatment.

In general, a poor response (design B and D, with less than 1 log 10 reduction at the end of the first week) leads to a poor estimation of the viral parameters. As expected, the study of the FIM reveals that the correlation between the estimators of EC 50 and n is important and increases with EC 50 . For design D 1 , good responders with low (resp. high) inter-dose variation exhibited a correlation equal to 0.25 (resp. 0.58) compared to 0.70 (resp. 0.77) in poor responders with low (resp.) inter-dose variation. This is also true for the correlation between c and δ that increases from -0.12 (resp. -0.08) to -0.71 (resp.

-0.49). The study of the eigenvectors of M F (ψ; D) reveals that the difficulty in practical identifiability stems from the ratio δ/c. In practice, parameter set-up associated with poor responders (i.e. a low δ)

would make the estimation of both parameters difficult. Since δ is known to have greater importance and to be more variable than c (3), it should be recommended to fix c. For instance, when fixing c and using D 1 with parameter set-up C, we found RE(δ) ≃ 0.12 as compared to 0.20 when estimating both parameters. This leads to a confidence interval CI 95% (δ) ≃ 0.10 ± 0.024 as compared 0.10 ± 0.04 when both parameters are included in the estimation procedure.

For the pharmacokinetic parameters, all the profiles except the parameter set-up D show comparable levels of precision. Indeed, this set of parameters is characterized by a rapid degradation of Peg-IFN, k e , and a low second order sensitivity n, and only limited information can be gained after day 3. Analyzing the variances of the random effects is less intuitive but follows the same trend already described for the fixed effects parameters. The SE of the error measurements are roughly proportional to the number of sampling points [START_REF] Retout | Fisher information matrix for non-linear mixed-effects models: evaluation and application for optimal design of enoxaparin population pharmacokinetics[END_REF], but our results indicate that the sampling times and the parameter set-up choices can have substantial effects.

It must be noted that the interpretation of the viral dynamics is seriously limited by the lack of practical identifiability for the parameters related to the host dynamics. Indeed, if the de novo infectivity, b, is included in the set of parameters to estimate, we found RE(β) > 20 in all cases (not shown).

Moreover, including b seriously impaired the estimation of the other parameters, particularly the loss of infected cells, δ, and resulted in a correlation larger than 0.8 in all cases, with RE(δ) > 0.18. This result reflects that the loss rate of infected cells and the infectivity rate share similar effects on the viral load.

Yet, mathematical simulations of the model show that the effects of these two parameters on the host dynamics are different, and thus observation of T and/or I would help to solve this indetermination.

Similarly, we studied the case in which ribavirin is effective by reducing the ability of virus to infect new cells, with 90% effectiveness (η = 0.90). We found that RE(η) > 3 and CI 95% (η) ≃ [0.30; 1.0] for all cases: thus, the ribavirin effectiveness is not practically identifiable. The non-identifiability of the ribavirin effect has already been discussed in empirical sensitivity studies (3; 19). However studying the FIM enables a rigorous assessment.

Fixed effects Random effects

Error measurements

Φ(D) M log(EC50) log(n) log(δ) log(c) log(ka) log(ke) log(V d ) ω 2 log(EC50) ω 2 log(n) ω 2 log(δ) ω 2 log(c) ω 2 log(ka) ω 2 log(ke) ω 2 log(V d ) σ 1 σ 2 A. Good

Design optimization

Computational limitations prevent extensively searching for an optimal study design over continuous sampling ranges in the first month of treatment. Nevertheless, study design can be optimized over a smaller grid consisting of all the potential sampling points described in D 1 -D 5 . Two constraints were imposed. First, both markers were measured at the same sampling time points, including t = 0. Second, the total number of samples for each marker was N * M = 240, where N is the number of patients and M is the number of sampling times for each marker; this value corresponds to the one used in the most parsimonious design, D 1 , described above (N = 30, M = 8).

Table 4 displays the optimal designs according to the maximal number of sampling times M for each marker and each patient. For instance, if M = 3, then the number of patients that can be enrolled in the study is N = 80. To take advantage of the population framework, the optimal design can be composed of different elementary designs where all the patients do not necessarily have the same number of measurements or the same schedule for the sample measurements. Predicting the best balance between the number of patients, N , and the number of samples per patient, M , is not straightforward. In other words, what should be preferred: rich information within each individual and poor information on between-patient variability, or vice-versa ?

Interestingly, the design efficiency, as measured by Φ(D), varies substantially according to this choice.

We found the best efficiency when M = 4 (and hence N = 60 enrolled patients). This design outperforms any other designs considered above. Indeed, while D 5 used twice as much information (N * M = 30 * 16 = 480 samples per marker), its efficiency is equal to 213.1 (Table 3) compared to 230.2 for our optimized design, i.e., a mean gain in precision of 1.05 while the raw amount of information to be collected is reduced by 50%. By using the population framework and the possibility of subgroups with different elementary designs, the fixed effects parameters are substantially improved. For instance, this optimal design gives RE(δ) = 0.057, versus RE(δ) ≥ 0.095 in all other contexts in Table 3. Hence, using this design, δ can be estimated with very good precision (CI 95% (δ) = [0.18; 0.22]). As expected, estimating the variances of the random effects is also strongly improved since there are N = 48 patients compared to N = 30 in A-D. However, since the number of points per patient is only equal to 4, the SE for the error terms is elevated, compared to previous cases. For instance, SE(σ 1 ) = 0.024, compared to SE(σ 1 ) = 0.0075 obtained with the same parameter set-up A and the design D 5 .

Although an intuitive approach could lead us to focus on early time points to estimate the initial viral kinetics, these results show that designs focused on the first week of treatment can generate redundant information and are not optimal. Furthermore, it would be preferable to enroll more patients and to complete long-term follow-up, rather than to conduct intensive sampling with fewer patients. Although this can imply a slight disadvantage in the estimation of treatment effectiveness, it is compensated by a more accurate estimation of the loss rate of infected cells δ, which is the most critical parameter for predicting treatment outcome [START_REF] Neumann | The 2nd phase slope of HCV decline is highly predictive of sustained virologic response following consensus IFN treatment for chronic hepatitis C and is determined by genotype but not dose[END_REF]. This is even more evident if the pharmacokinetic parameters are not constant over time, but could change over time as found in (5; 29).

Conclusion

We have proposed a methodology for rationalizing the measurement sampling times in models of HCV dynamics. Although some limitations for EC 50 in sparse designs have been pointed out, this methodology has been shown to provide good and rapid approximations of the standard errors in this challenging context.

We have shown that an appropriate choice of the sampling times and the elementary designs allows for a substantial reduction of the information needed to estimate the parameters. In particular, we showed that, as compared to intensive designs used in previous studies, the number of sampling times could be reduced by up to 50%. Furthermore, parameters could be obtained with excellent precision with only four measurements of each marker. Interestingly, these data should not be focused on only the very first days that follow treatment initiation, but rather should be balanced throughout the entire follow-up period.

Indeed, long-term sampling points (third and fourth week) provide substantial information, particularly for investigating the viral parameters that could be responsible for a poor response to treatment.

Although population approaches advantageously use the between-patient variability to optimize inference, we showed that the measurement of only the viral load and the drug level is not sufficient for practically estimating all the parameters of the current model, in particular those related to the dynamics of the host and to the antiviral effectiveness of ribavirin [START_REF] Guedj | Practical Identifiability of HIV Dynamics Models[END_REF]. This feature highlights more generally that in the absence of data on the host dynamics, the understanding of HCV viral kinetics provided by modeling approaches remains speculative.

There is a paradox in the approach based on the study of the FIM. It is assumed that the values of the parameters are known, whereas the aim of the study is ultimately to estimate these parameters.

Hence, it could be argued that the figures found here are valid only for the specific parameter set-up that was chosen. While it is true that the precise values of the parameters are unknown, their order of magnitude, in general, is guaranteed. Indeed, these parameters have a biological and physiological meaning and, hence, they must take on values within certain ranges. We have proposed to partially circumvent this issue by showing that the main conclusions were consistent when changing the values of the fixed effects parameters. Another approach is to use a priori distributions for the parameters (36; 37), but a major drawback of Bayesian approaches is the computational cost, since the expectation of the posterior variance-covariance matrix of the parameters must be (numerically) evaluated.

This work constitutes a first step towards improved rationalization of data sampling designs for HCV viral kinetics studies, but several extensions of this work are necessary. First, the information brought by data under the limit of detection (lod), whose proportion is increased with the development of new, highly potent, direct-acting antiviral agents, needs to be rigorously taken into account (38; 39). Another natural extension of this work is to take into account the constraint of the financial cost of patient [START_REF] Retout | Design optimization in nonlinear mixed effects models using cost functions: application to a joint model of infliximab and methotrexate pharmacokinetics[END_REF]. Lastly, predicting treatment outcome is done at the individual level, and the precision at both the population and the individual levels are important. New approaches need to be developed to study the appropriate trade-off between these two levels of prediction. Simulation of the median viral dynamics (black) and Peg-IFN profile (dashed) according to the four parameter set-ups given in Table 2 A) Good responders with low inter-dose variations. B) Good responders with high inter-dose variations. C) Poor responders with low inter-dose variations. D) Poor responders with high inter-dose variations. For each case, we indicate the data points where the viral load is sampled according to the five designs given in Table 1:

is D 1 , • is D 2 , is D 3 is D 4 , is D 5 .
For the sake of readibility, the symbols are drawn above and under the curve and are not indicated for the Peg-IFN profiles since the measurements are done at the same time point.

FIGURE 2

Difference between the approximated and the empirical standard errors according to the study design. The difference is denoted by ∆(SE) = SE P F IM -SE EM P and is displayed separately for the fixed effects parameters (left) and the variance parameters (right). The vertical lines define the area where -0.01 < ∆(SE) < 0.01 and the symbol "•" indicates that ∆(SE) = 0 (up to three significant digits); the code color for the study designs is similar to Figure 1.

  ni the elementary design of the measurement times j for patient i, and D = {d i } i=1,...,N the design of the study (for all patients). Because of measurement errors, random fluctuations, and model limitations, the observations are not the exact output predicted by the mathematical model. We ensure a Gaussian framework and homoskedasticity of the error measurements by working on the log-transformed values of the viral load data. We summarize all sources of deviations by introducing for each observation an independent and additive white noise process, ǫ k , with variances σ 2 k . We set σ = (σ 1 , σ 2 ) T , where k = 1 indicates HCV RNA and k = 2 indicates Peg-IFN. Furthermore, the model for the observations is:

  and the global observed log-likelihood is L(ψ; Y ) = i=1,...,N L i (ψ; Y i ) by independence between patients. The elementary Fisher Information Matrix (FIM)

  responders with low or large inter-dose variations (A, B), and poor responders with low or large inter-dose variation (C, D) (Figure1). Large inter-dose variations (B, D) are typically observed with Peg-IFN-α2b treatment (5; 28), whereas low variations (A, C) are more likely to be observed with Peg-IFN-α2a treatment (characterized by a slower absorption rate but also a slower clearance rate)(28; 29).

1 -D 5 ,

 15 500 data set replications of N = 30 patients with parameter set-up A. For each data set, parameters were estimated using MONOLIX version 2.4 (http://software.monolix.org), a software program based on a Stochastic Approximation Expectation-Approximation (SAEM) algorithm[START_REF] Kuhn | Maximum likelihood estimation in nonlinear mixed effects models[END_REF]. Parameters initial values used with the SAEM algorithm were the true values of the parameters. The empirical SE were obtained by computing the standard deviation of the 500 estimates. Figure2displays, for each design and each parameter, the difference ∆(SE) = SE P F IM -SE EM P , which is a measure of the error made in the estimation of the confidence interval amplitude. We found that PFIM provides, in this challenging context, a very satisfactory approximation of the empirical SE, with ∆(SE) < 0.01 for all fixed effects and all designs, except in three cases: EC 50 with D 1 and D 3 , and c with D 1 . Moreover, the gain in computation speed is considerable: for each design, PFIM takes only about 1 minute, as compared to > 5 days of computation needed to get the empirical standard error (on an Intel Core 2 Duo 2.5 GHz).

  FIGURE 1

Table 1 :

 1 Design used in five studies of viral dynamics. M is the number of sampling points for each marker and each individual

	Design name Authors	Measurement times (days)	M
	D 1 D 2 D 3 D 4 D 5	Zeuzem et al. (2005) (24) Sherman et al. (2005)(25) Herrmann et al. (2003)(26) {0, 0.25, 0.5, 1, 2, 3, 4, 7, 10, 14, 21, 28} {0, 1, 4, 7, 8, 15, 22, 29} {0, 0.25, 0.5, 1, 2, 3, 7, 10, 14, 28} Zeuzem et al. (2001)(27) {0, 0.040, 0.080, 0.12, 0.20, 0.33, 1, 2, 3, 4, 7, 14, 21, 28} 14 8 10 12 Talal et al. (2006)(5) {0, 0.25, 0.5, 1, 2, 3, 5, 6, 7, 7.25, 7.50, 8, 9, 14, 15, 16} 16

Table 2 :

 2 Values of the parameters used in the simulations (on their natural scale), according to four different patterns of viral decline in the population (Fig 2)

Table 3 :

 3 Expected SE for a sample size of 30 patients, according to five popular designs and four different sets of parameter values, corresponding to four patterns of viral kinetics during Peg-IFN based therapy.

Table 4 :

 4 Optimal designs according to the maximum number of samples M allowed per patient for a fixed number of N * M = 240 sample measurements for each marker, similarly to design D 1 (Table1). The parameter values correspond to the set-up A (Table2). (sampling times),number of patients per group } log(EC 50 ) log(n) log(δ) log(c) log(k a ) log(k e ) log(V d ) Although increasing the number of patients improves the identifiability of the parameters, it often increases the cost of the study. Recently, a new development of PFIM integrates the financial cost of study designs into the design optimization procedure but this has not been yet implemented in the available version of PFIM

	Optimal design

Acknowledgments

We thank L. Murillo, A. Smith and A.S. Perelson (LANL) for helpful comments; we are also grateful to both anonymous reviewers and the AE for very insightful remarks. This work was supported by the French Ministry of Foreign Affairs, the Israel Science Foundation (ISF 939/2008), and the US Department of Energy under contract RR06555 and P20-RR18754.