Sébastien Chabod
email: sebastien.chabod@lpsc.in2p3.fr

Julien Giraud

Marine Hervé

Daniel Santos

Nadine Sauzet

Optimal heavy water neutron moderators for an AB-BNCT treatment unit

Keywords: boron neutron capture therapy, neutron moderator, topology optimization

In this study, we use a topology optimization algorithm, developed at the CNRS LPSC, to design a heavy-water neutron moderator for a BNCT treatment unit. The solutions generated by this algorithm are compact yet succeed in limiting the exposure of patient's healthy tissues to levels below recommended limits. They present subtle, original geometries inaccessible to standard design techniques. The versatility of this novel approach makes it possible to automate the moderator design, and fit it to the configuration of the BNCT unit considered, e.g. the neutron source and materials it uses or the room it occupies, as well as to the biological parameters of its patients, e.g. the volume and depth of the tumor to be treated or the characteristics of the targeted organs.

Introduction

Context. Boron Neutron Capture Therapy (BNCT) is a radiotherapy technique that uses a flux of epithermal neutrons to treat deep and/or diffuse tumors [START_REF]Current status of neutron capture therapy[END_REF]. In the early days, these fluxes were generated by nuclear reactor cores, limiting the use of this approach to a small number of sites, the number of which is also declining with the decommissioning of a growing number of research reactors. The development of compact accelerator-based neutron sources that can be deployed in a hospital environment is removing this limitation, leading to renewed interest in this approach, renamed Accelerator-Based BNCT (AB-BNCT). Concretely, a patient absorbs before treatment a boron-10 delivery agent, which targets preferentially cancer cells [START_REF]Current status of neutron capture therapy[END_REF]. Boron-10 captures neutrons with a high cross-section (~3800 barns at 0.025 eV) and splits into 2 nuclei, 4 He and 7 Li, with a probability of 94% both carrying 2.3 MeV of kinetic energy and 6% carrying 2.78 MeV. These heavy charged particles stop after a short track (~5-8 µm) through tissue, depositing most of their energies in the cancerous cells, destroying them. Some amount of boron-10, however, remains in the patient's blood, leading to unwanted exposure of healthy tissue subjected to the neutron flux. In addition, the primary and secondary particles that propagate in the treatment unit and in the patient induce potentially damaging exposure, which should be kept to a minimum [START_REF] Coderre | [END_REF]. For these reasons, the AB-BNCT neutron sources, which emit mostly epithermal to fast neutrons, must be assisted by moderators, which lower the energies of neutrons to levels less risky for healthy tissue, <~ 10 keV [START_REF]Current status of neutron capture therapy[END_REF]. The scientific and technological challenges that must be addressed to make AB-BNCT available, efficient and safe, therefore include (1) the development of intense compact neutron sources, (2) the development of targets that can withstand a high beam power (~30-75 kW), necessary to reduce the treatment time to the order of one hour; (3) the synthesis of boron delivery agents allowing a higher contrast between the concentration in 10 B in the tumor and the blood; (4) improving knowledge on neutron and photon transport, and on the calculation of deposited doses and their biological effects; [START_REF] Hervé | Fast and epithermal neutron fields for Accelerator Based Neutron Capture Therapies[END_REF] high-accuracy monitoring of the doses delivered to patients; [START_REF] Capoulat | [END_REF] the development of compact neutron moderators capable of adjusting their energies and directions to target a tumor while best preserving healthy tissue. The present study addresses this last issue [START_REF] Capoulat | [END_REF]. Aim of the study. The design of BNCT moderators currently relies on parametric approaches. First, a set of potential configurations is assembled, by varying the composition and dimensions of their components. The choice of these configurations is based in part on the analysis of neutron-matter interaction cross-sections, in part on human intuition and experience. For each possible configuration, a simulation of the treatment unit is then performed, using Monte-Carlo codes nowadays. Moderator configurations that offer interesting compromises between chosen therapeutic and economic objectives are thus identified, and possibly refined. These parametric studies yield interesting, multi-material design solutions, which deliver high tumor doses while keeping the doses deposited in healthy tissues under control [3][START_REF] Capoulat | The 9 Be(d,n) 10 B reaction as a neutron source for Boron Neutron Capture Therapy[END_REF][START_REF] Hervé | Fast and epithermal neutron fields for Accelerator Based Neutron Capture Therapies[END_REF]. However, being constrained by human creativity, they can only explore a small fraction of the solution space, and may thus miss out on more efficient but unintuitive designs. For intermediate to deep-seated tumors e.g., parametric based designs cannot yet prevent the deposition of large doses in upstream tissues, in particular in the scalp during the treatment of glioblastoma. Therefore, an approach that would allow the space of possible configurations to be explored with little to zero human bias could provide the BNCT community with innovative design solutions. For this purpose, we will leave in this study the determination of the design of a heavy-water neutron moderator in the hands of a topology optimization (TOPOPT) algorithm, the principle of which will be described in section 3.1. Examples of successful applications of this TOPOPT approach will be presented in sections 3.2 and 3.3. Finally, we will conduct in section 4 a set of sensitivity studies required to evaluate the robustness of the TOPOPT computations, before concluding in section 5. Computation means. The simulations carried out for this study were run on servers with 24 modern CPUs (+24 virtual), dating from the end of the 2010s, at the rate of one simulation per server. The computation times mentioned in the rest of this study correspond without exception to the durations of the computations launched on these machines.

Modeling of the treatment unit

In this section are presented the models, data and codes used for computing the structure of a BNCT neutron moderator and its performances, for a specific configuration described below. Neutron source. For this study, we consider one of the main neutron sources envisaged for an AB-BNCT treatment unit, which makes use of the 9 Be(d(1.45 MeV), n) reaction. This reaction produces neutrons with good efficiency, 3.3×10 11 neutrons/mC, with a not too hard energy spectrum, and its manufacture and operation at high power are considered less challenging than 7 Li(p,n) sources [START_REF] Capoulat | The 9 Be(d,n) 10 B reaction as a neutron source for Boron Neutron Capture Therapy[END_REF]. The deuteron beam, here parallel to the Y axis, is assumed to be uniformly distributed over the surface of the 9 Be target, a 10 cm 2 disc visible fig. 1 at x = -6.7 cm. Such a surface is considered necessary to dissipate the high power deposited by the beam, expected to be around 40 kW. The 9 Be deposit is assumed to be homogeneous. Therefore, in our simulations, the neutron emission coordinates are sampled uniformly across the surface of the target disk. The energy and angle distribution of the generated neutrons was measured by Capoulat et al. [START_REF] Capoulat | [END_REF]. It is rigorously implemented in the transport code used for this study, MCNP 6, by using the functionalities of its SDEF card [7]. The target casing model, simplified for study purposes, is an aluminum and graphite structure, in yellow fig. 1. Between the casing and the moderator is inserted a layer of lead, in blue fig. 1, to reduce the fluence of primary gamma rays. Neutron moderator. A 9 Be(d(1.45 MeV), n) source generates neutrons of ~2 MeV in average, which must be slowed down to the desired energies for a BNCT treatment, <~ 10 keV, in order to reduce the neutron dose in the healthy tissues. This moderation must be carried out in a device with a low footprint, compatible with use in a hospital environment. Such a moderator must therefore be manufactured mainly using materials having light nuclei, e.g. 1 H (water, plastics), 2 H (heavy water), Be, C, F (inelastic diffusion), etc. For this study, we will optimize the structure of a moderator made of heavy water, 99% atomic pure. Heavy water is an interesting material for several reasons: (i) it has excellent neutron moderating power, 2 H having a low nucleus mass, close to that of the neutron; (ii) 2 H has a low neutron capture cross-section, which increases the output neutron fluence and reduces the production of dosing capture gamma rays. This material has already been used as a moderator in BNCT installations, at the JRR-4 reactor for example [START_REF] Torii | The BNCT Irradiation Facility of JRR-4[END_REF]. Its use was recently studied by Hervé et al., who showed during a parametric study that a hemispherical moderator made of heavy water presented the best tumor dose over brain dose ratios for most of the parameters they explored (composition and radius of the moderator, depth of the tumor) [START_REF] Hervé | Fast and epithermal neutron fields for Accelerator Based Neutron Capture Therapies[END_REF]. This moderator, with a maximum radius of 20 cm and a length of ~30 cm (variable), is inserted into a heavy concrete wall, in green fig. 1, which serves as biological protection. In order to limit the computation time, this moderator has an axis of symmetry, coinciding with the axis of the beam. Unless otherwise stated, the moderator will be segmented into N = 375 voxels Θj = 1…N (25 steps in length, 15 in radius). These voxels, drawn in fig. 1, are cylindrical crowns with axis that of the beam, each containing a density ρj of heavy water varying between 0 and 1.11 g/cm 3 , the natural density of heavy water. The objective of the study is to find the densities ρj that optimize the BNCT treatment quality, using a topology optimization algorithm. Patient's modeling. For this study, we will consider the case of the treatment of glioblastoma. We will only model the patient's head, simplifying its composition. Only the brain, skull and scalp are modeled, using spheres with radii 7.05 cm, 8.21 cm and 8.71 cm, visible in fig. 1. These average radii, consistent with anatomical data for a standard adult human [START_REF] Bushby | Centiles for adult head circumference[END_REF]10], give volumes of brain, skull and scalp equal to those of Snyder's model [11]. The Clinical Target Volume (CTV) containing the tumor is modeled by a sphere of volume 20 cm 3 , represented in red fig. 1. This sphere is positioned at a depth Ptum, defined as being the smallest distance between the center of the CTV and the surface of the patient's head, cf. fig. 1. Its center is on the axis of the beam. The center of the patient's head is also on the axis of the beam, at a distance chosen so that the minimum distance between the surface of the head and the right side of the moderator is equal to 1 mm, cf. fig. 1 (quasi contact to minimize neutron losses). The patient's head is subdivided into M = 89 voxels ∆i = 0…88: 1 voxel ∆0 for the tumor, 88 others for healthy tissues (64 for the brain, 16 for the skull, 8 for the scalp). These 88 voxels are delimited by: (i) cones whose vertices are the center of the CTV and whose opening angles vary between 0° and 180°; (ii) and ellipsoids of revolution, whose axes of symmetry coincide with that of the beam, and whose ends are chosen to regularly pave the patient's head. This paving is presented in fig. 1 (right). The choice of this head model is motivated by several constraints: (i) the computation time: the time required for the optimization algorithm to converge increases with the number M of voxels used to pave the head, in ~ 4M for M small. For M ≫1, this time is expected to increase exponentially with M. Thus, for 89 voxels ∆i for the head and 375 voxels Θj for the moderator, it takes 2 months of computation on a modern 24 CPUs server to determine the optimal structure of the moderator, a long time but which remains humanly compatible. Without axial symmetry, it would be necessary to add a discretization in angle ϕ around the axis of symmetry, by cutting the structure using planes containing the axis of the beam. By taking for example 10 angles ϕ, the number M of voxels would increase to 881, and the computation time to ~20 months; (ii) uncertainty on the volume of tissue concerned by a peak dose: axial symmetry eliminates the uncertainty on the volume of tissue affected by a peak dose, most often not provided in the literature. Indeed, by cutting the model of the patient's head using planes containing the axis of symmetry, again, we can subdivide each voxel ∆i into an arbitrarily large number of sub-voxels, of arbitrarily small volumes but nevertheless all exposed to the same dose; (iii) sensitivity to morphological variability: analytical models, e.g. of the Snyder type, are quite rigid, and their representativeness for all patients may raise debate. In order to study the impact of the morphological variability of patients on the quality of their treatments, the head model considered in this study can be of interest. This sensitivity study will be conducted in section 4.3. Biological dose. To evaluate the biological doses deposited in each voxel of the head, we will use the ICRU 46 compositions and densities of the tissues, recalled in Table 1 [11]. The composition of the tumor will differ from that of the brain only in its higher concentration of 10 B. The 10 B delivery agent considered is BPA, the corresponding 10 B concentrations CB are given in Table 2 for each tissue (factor 1 brain/blood, 3.5 tumor/blood, 1.5 scalp/blood) [START_REF] Capoulat | The 9 Be(d,n) 10 B reaction as a neutron source for Boron Neutron Capture Therapy[END_REF]. In this study, we will use a standard definition of the total biological dose, D = wγDγ + wnDn + wpDp + wBDB, where Dγ,n,p,B are the gamma, neutron, proton and boron doses, and wγ,n,p,B the corresponding biological weighting factors, given in Table 2 for each tissue [START_REF] Coderre | [END_REF]. To evaluate each dose component Dγ,n,p,B, we will use kerma factors for each tissue. The supplementary material of Ref. [11] contains files which allow to reconstruct the neutron kerma, Kn(E), using the elementary neutron kerma given Ref. [12,[START_REF]Nuclear data for neutron and proton radiotherapy and for radiation protection[END_REF], for each tissue. These kerma factors notably take into account the contributions of the 14 N(n,p) or nat Cl(n,p) reactions. The Kn curves obtained are presented in fig. 2. The Kγ photon kerma are reconstituted for each tissue using Ref. [START_REF][END_REF] and log-log fits for the missing energy intervals. We checked fig. 3 (left) that this approach returns the photon kerma for the brain provided in Ref. [START_REF] Kerr | Photon and neutron fluence-to-kerma conversion factors for ICRP-1975 reference man using improved elemental compositions for bone and marrow of the skeleton[END_REF]. The photon kerma obtained for each tissue are shown fig. 3 (right). Finally, the KB boron kerma comes from [11], and is given in fig. 2 per ppm of 10 B in tissue. The neutron kerma Kn(E) including the contributions of (n, p) reactions, and the Coderre coefficients wn and wp being equal, the contributions Dn and Dp will be combined in a single dose, Dn. The equation for the total biological dose Di in each voxel ∆i of the patient's head is then written:

( ) ( ) ( ) ( ) ( ) , , , , , , , , , , , , , , i
B i B i n i n i i i B i B i n i n i n i i i i E D w D w D w D w K E w K E E w K E dE γ γ γ γ γ φ φ = + +   = + +   ∫ (1)
where φn,i(E) and φγ,i(E) are the neutron and photon fluences per unit of energy in the voxel ∆i. Gy-Eq (brain) and 16.7 Gy-Eq (scalp), also without justification nor indication of the volume concerned [START_REF] Capoulat | The 9 Be(d,n) 10 B reaction as a neutron source for Boron Neutron Capture Therapy[END_REF]. An average peak dose for the whole brain, 7 Gy-Eq, is however mentioned there. Hervé et al. used L = 14.1 ± 1.8 Gy-Eq per cm 3 of brain [START_REF] Hervé | Fast and epithermal neutron fields for Accelerator Based Neutron Capture Therapies[END_REF]. This value comes from a study by Coderre et al., which showed that exposure of 14.1 Gy-Eq in 1 cm 3 of brain induces in 50% of treated patients a somnolence syndrome [START_REF] Coderre | [END_REF]. This proportion is reduced to 5% for L = 11 Gy-Eq (1 cm 3 of brain) [START_REF] Coderre | [END_REF], a more acceptable proportion if BNCT were to develop, which is the value used by Capoulat et al. , which shows that doses to the skin between 16.5 and 24 Gy-Eq have manageable toxicities [START_REF] Menéndez | [END_REF]. The study by Menéndez et al. however concerns a small number of patients, whose area treated with BNCT was not the scalp but various leg skins. The transposition of the results of [START_REF] Menéndez | [END_REF] to the determination of L (scalp) is therefore not immediate.

In conclusion, as a precaution, we will use in this study the most conservative L values proposed in the literature, i.e. 11 Gy-Eq (brain), 30 Gy-Eq (skull) and 16.7 Gy-Eq (scalp). These conservative values pose an interesting challenge to the design of the moderator, which should ensure a high contrast between the tumor dose and the doses deposited in healthy tissues.

Transport code and nuclear data. The simulations carried out for this study were performed with the transport code MCNP 6 [7], simultaneously propagating neutrons and photons. With the exception of the results presented in section 4.2, the neutron and photon transport data used come from the ENDF/B-VII.0 database (/B-VI.6 for 1 H and 138 Ba). All cross-sections, including S(α,β) data, are taken at room temperature. With the exception of calculations performed with N = 150, all MCNP simulations performed in this study were run using 5×10 9 source neutrons.

Topology optimization of the moderator

In this section, we will let a topology optimization algorithm determine the structure of an AB-BNCT heavy water neutron moderator, with a level of human intervention as low as possible.

Still, choices will have to be made, which like any human choice are debatable, and amendable 1 if necessary. We will choose to optimize the structure of a D2O moderator, favoring, over all 2 other considerations, the cleanliness of the treatment, i.e. the ability of the AB-BNCT unit: (i) 3

to generate a high dose in the tumor; (ii) while ensuring that the doses deposited in the patient's 4 head do not exceed the local peak doses, or exceeds them as little as possible if this goal proves 5 untenable. This choice, which gives priority to the therapeutic objective over the economic 6 objective, will lead to moderator configurations that will not necessarily maximize the dose rate 7 deposited in the tumor, and therefore will not minimize the treatment time. This choice is 8 present in the watermark in many studies of moderator design, cf. e.g. [3][START_REF] Capoulat | The 9 Be(d,n) 10 B reaction as a neutron source for Boron Neutron Capture Therapy[END_REF][START_REF] Hervé | Fast and epithermal neutron fields for Accelerator Based Neutron Capture Therapies[END_REF]. In part because 9 treatment time is function of the maximum beam power, which itself depends on manufacturing 10 processes and technological innovations not all off the shelf. Partly also because there is no 11 rigid guideline on this treatment time. Some recommend it to be less than 60 min, but what 12 about the possibility of multi-fractionating the treatment for example. In the absence of certainty 13 on this subject, the choice made in section 3 seems reasonable. It is however important to note 14 that an optimization of the moderator design with a different objective, e.g. an objective of 15 minimizing the treatment time, can be carried out using the same methodology than that 16 presented in section 3.1, by simply replacing the FPEN function in Eq. ( 2) by the treatment time. 17 We want to determine the optimal composition of a D2O moderator, i.e. the densities ρj of heavy [START_REF] Icesl Slicer | Progressive infils[END_REF] water contained in each voxel Θj, allowing a Dobj dose to be deposited in the tumor while 22 minimizing the exposure of healthy tissues. To achieve these goals, we must ensure that each 23 dose Di, given Eq. ( 1), deposited in each voxel ∆i of the patient's head does not exceed the 24 corresponding peak dose Li. As a reminder, Li = 11 Gy-Eq ∀∆i ∊ brain, 30 Gy-Eq ∀∆i ∊ skull 25 and 16.7 Gy-Eq ∀∆i ∊ scalp. This problem can be formalized as follows: 26 27 tumor must be equal to Dobj, equal to 30 Gy-Eq for a single-fraction glioblastoma treatment, cf. 39 pp 201 of [START_REF]Current status of neutron capture therapy[END_REF]; (C4) the D2O densities ρ can vary between 0 and ρmax = 1.11 g/cm 3 , the natural 40 D2O density. The FPEN penalization function in problem (2) quantifies the cleanliness of the treatment: the lower its value, the more the treatment respects the peak doses in healthy tissues. Concretely, for each voxel ∆i > 0 of the patient's head, FPEN is incremented by a quantity (Di/Li) n weighted by the volume Vi of the voxel, Vhead being the volume of the head. By taking a number n≫1, FPEN is therefore incremented by a small quantity if Di < Li, since x n →0 if x < 1, n→+∞, and on the contrary by a quantity high if Di > Li, since x n →+∞ if x > 1, n→+∞. We will take n = 30 in this study. An analysis of the impact of this choice is carried out in section 4.1. Problem ( 2) is a constrained optimization problem of which we can see the extreme complexity. The FPEN function is a non-linear functional of the angular fluences of neutrons and photons, themselves non-linear, highly complex functionals of the densities ρ of the moderator. The space of solutions to be explored is also gigantic: for 375 voxels, each containing a density varying between 0 and ρmax, say by small step δρ = ρmax/50, the number of possible moderator configurations is worth 50 375 = 1.3×10 637 , that is to say a numerical infinity. Despite this extreme complexity, it is now possible to solve problem (2) in a humanly compatible time using the topology optimization (TOPOPT) algorithm presented in Ref. [18,19]. This algorithm solves the Lagrangian equations associated with the optimization problem, given Eq. ( 3), using an iterative procedure. Starting from a uniform configuration of densities, ρj = ρ0 ∀j, it gradually modifies, small step ±δρ by small step, the densities ρj of the voxels Θj until reaching the solution of the optimization problem at convergence. The doses Dn,i, Dγ,i and DB,i which make up each dose Di, cf. Eq. ( 1), can be calculated with MCNP using its DE and DF cards [7]. The calculation of the derivatives ∂FPEN/∂ρj involves M×N derivatives ∂Di/∂ρj, which can all be calculated in a single MCNP simulation using its PERT card [7] and the procedure described in section 1 of Ref. [18]. In this study, we will take δρ = ρmax/50 and ρ0 = 0.2 g/cm 3 . . 0, 0

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 1 0 4 1 2 0 3 max 4 min , , , , , 1 , , , , , s.t.
(C ), (C ) (C ) (C ) i M n i PEN i i i i head i B B n n n i i E r n n n n obj j V F D L V w r K r E r E w r K r E D d rdEd V w r K r E r E Β Q Β Q D D j ρ π γ γ γ γ γ γ ρ α ρ α ϕ ρ ρ ϕ ρ ρ ϕ ρ ϕ ϕ ρ ρ ρ = +∞ = Ω∈ ∈∆ = =       Ω     +   = Ω     + Ω     = = = ≤ ∀ ∑ ∫ ∫ ∫ (
PEN L F C L L λ ρ λ = - ∂ ∂ = = ∂ ∂ (3) 3.2

. Example of application of the TOPOPT approach

To illustrate the operation of this TOPOPT algorithm, we propose to determine in this section the optimal composition of a D2O moderator for a 20 cm 3 CTV at a depth Ptum = 4.5 cm, cf. section 2. We will take for this calculation N = 375 voxels for the moderator and M = 89 voxels for the patient's head. The thickness H of the moderator and the heavy concrete wall is set at 30 cm, cf. fig. 1. For this configuration, the algorithm converges in 95 iterations, which required 2 months of calculation. Fig. 4 (left), we present in the XY plane some moderator configurations obtained before and after convergence for iteration numbers NUMITER ranging from 0 (starting configuration, ρj = 0.2 g/cm 3 ∀j) to 95 (convergence). The 3D structure of the moderator is generated by rotating these graphs around the Y = 0 axis of symmetry. The gray scales in fig. 4 (left) give the densities ρj in D2O, in g/cm 3 , in the voxels Θj of the moderator. Fig. 4 (center), we present the maps of the total doses Di associated with these moderators, deposited in each voxel ∆i of the patient's head. The color scales in fig. 4 (center) give the values of the Di doses in Gy-Eq. Finally, fig. 4 (right), we present the values of the dose excesses, Ei = (Di -Li)/Li, in each voxel of the head, where Li is as a reminder the value of the peak dose in the voxel ∆i. The color scale in fig. 4 (right) gives the values of Ei in percent. For Ei <0, the objective of the moderator optimization, the voxel is colored white. So, in summary, the more white and blue in these maps, the cleaner the treatment. These figures show how the TOPOPT algorithm improves iteration after iteration the structure of the moderator to improve the treatment cleanliness. To complete this information, we show fig. 5 the evolution of FPEN (left) and of the total dose rate Dtum deposited in the tumor, expressed in fGy-Eq/source neutron (s.n.) (right), as a function of the iteration number NUMITER. At convergence, Dtum = 1.05 fGy-Eq/s.n.. The efficiency of a 9 Be(d(1.45 MeV), n) source being 3.3×10 11 neutrons/mC [START_REF] Capoulat | The 9 Be(d,n) 10 B reaction as a neutron source for Boron Neutron Capture Therapy[END_REF], for a beam intensity of 30 mA [START_REF] Capoulat | The 9 Be(d,n) 10 B reaction as a neutron source for Boron Neutron Capture Therapy[END_REF], the treatment time necessary to deposit 30 Gy-Eq in the tumor will be 48 min. This time could be further reduced by adding to the simplified model of the unit treatment shown fig. 4 a reflector around the source at x < 0 cm. Indeed, a fraction of the source neutrons are emitted at angles greater than 90°, and are therefore artificially lost. [3,[START_REF] Capoulat | The 9 Be(d,n) 10 B reaction as a neutron source for Boron Neutron Capture Therapy[END_REF] described in section 2. For a dose to the tumor of 30 Gy-Eq and a CTV of intermediate depth, Ptum = 4.5 cm, we find that the doses deposited in healthy tissues remain mostly within the imposed limits. They exceed the peak values used by Capoulat et al. [START_REF] Capoulat | The 9 Be(d,n) 10 B reaction as a neutron source for Boron Neutron Capture Therapy[END_REF], 11 Gy-Eq (brain) and 16.7 Gy-Eq (scalp), only in a small volume of the brain (5.4%) and scalp (11.8%), and this by little in both cases. They nowhere exceed the peak values used by Torres-Sánchez et al., 12.5 Gy-Eq (brain) and 24 Gy-Eq (scalp), with a comfortable margin for the scalp. As for the skull, the maximum deposited dose is 9.9 Gy-Eq, well below its peak value, 30 Gy-Eq. In the end, the doses deposited exceed the peak doses of Capoulat et al. only in 4.8% of the head volume. To achieve these conclusive results with a compact moderation volume, 0.04 m 3 (20 cm radius, 30 cm in length) despite the high energy, ~2 MeV in average, of the 9 Be(d,n) source neutrons, the TOPOPT algorithm generates a subtle device, visible fig. 4 iteration 95 or fig. 7 in more detail, composed of: (i) a neutron guide which drives the neutrons towards the tumor while slowing them down. This guide contains a heavy water ring, which scatters the neutron trajectories. At its exit, it recreates in an original way the equivalent of a multi-field treatment, which consists in simultaneously exposing a patient to several moderated neutron beams with different orientations [3]; (ii) a central needle, which acts as a small guide of moderated neutrons, and adds another directional exposure; (iii) and a main moderation body and a reflector. The details of this sophisticated structure are inaccessible to a parametric study, let alone to human intuition, thus illustrating the potential of the TOPOPT approach. For the moderator obtained fig. 7, we give fig. 8, for each voxel ∆i of the patient's head, the weights wB,iDB,i/Di, wn,iDn,i/Di (all neutron contributions except 14 N(n,p)), wp,iDp,i/Di (only 14 N(n,p)) and wγ,iDγ,i/Di, of each of the 4 components of the total dose Di deposited, presented in section 2. The color scales of these figures give these weights in percent. We observe that: (i) the boron dose represents 76% of the total dose deposited in the tumor. Outside the tumor, however, it is only really important in areas of the scalp on the beam axis, where it matches the neutron dose contribution; (ii) the negligible nature of the 14 N(n,p) dose, except in the regions of the skull facing the moderator; (iii) the dominant character of the neutron dose, which accounts for ~50% of the total dose in all voxels except tumor, and which logically dominates in the areas of the scalp close to the exit of the neutron guide. The cleanness of the treatment thus strongly depends on the way in which the moderator slows down and focuses the neutrons on the head; (iv) the low relative contribution of gamma rays in the areas of the brain facing the moderator, which becomes important in its rear areas. This observation confirms the necessity to propagate neutrons and photons simultaneously in the simulations.

n. on the interval [E, E+∆E],

where the ∆E are the widths of the bins of the histograms fig. 9. We note that the moderator does not generate a ~10 keV mono-energetic, spatially homogeneous spectrum, contrary to what is sometimes recommended in the literature [START_REF]Current status of neutron capture therapy[END_REF]. Its neutron guide (position 3) brings a contribution between 1 keV and 200 keV, visible fig. 9. The spectra generated by the central needle (position 1) and the moderation body (position 2) are epithermal from 1 eV to 100 keV, and contain a smaller proportion of neutrons of ~1 MeV. The spectrum emitted by the reflector (position 4) is a little harder, which makes it possible to compensate for the greater quantity of tissues to cross so that the neutrons arrive at the tumor with an energy favoring the captures on the 10 B. It is thus observed that the shape of the spectra at the exit of the moderator is coupled with the details of its structure. The TOPOPT algorithm cleverly adjusts the energies and directions of the neutrons exiting the moderator to minimize exposure of healthy tissue.

Exploration of several parameters

In this section, we will show how the shape of the moderator and the quality of treatment evolve with the maximum available thickness H or with the depth Ptum of the tumor. These studies will illustrate one of the strengths of the TOPOPT approach, which allows it to fit the design of the moderator to the various possible configurations of a BNCT treatment unit, as well as to the various biological parameters of its patients, including the nature of the neutron source used, the volume and the depth of the tumors to be treated, the structures of patients' heads, or the time and possible fractionation of the treatment, among many other parameters. 3.3.1. Evolution of the structure and efficiency of a TOPOPT D2O moderator with its thickness The calculations presented in section 3.2 were carried out for a thickness, H = 30 cm, of the moderator and concrete wall. This thickness has not been set arbitrarily: its order of magnitude comes from the work of Hervé et al. [START_REF] Hervé | Fast and epithermal neutron fields for Accelerator Based Neutron Capture Therapies[END_REF], who showed that a heavy water hemisphere of ~30 cm in radius moderates neutrons enough to deposit a high dose in the tumor while limiting the exposure of healthy brain tissue. However, one can wonder how the shape and performances of an optimal D2O moderator evolve with thickness H. Fig. 10 (left), we present the TOPOPT configurations obtained for N = 375 voxels and Ptum = 4.5 cm, by varying H between 20 and 40 cm. The corresponding maps of the total doses Di (resp. dose excesses Ei in each voxel of the head are presented in fig. 10 (center) (resp. right). We observe fig. 10 a foreseeable evolution of the optimal shape of the moderator with H. For low thicknesses, H ≤ 20 cm, the too small available volume does not allow the algorithm to keep the central needle visible fig. 7, which is replaced in favor of a compact moderation body. Conversely, for high H values, there is more volume than necessary to slow down the neutrons. The algorithm has more room, and can create more complex and more efficient structures. In all cases, whatever the value of H, we observe that the algorithm systematically recreates a neutron guide similar to that shown fig. 7, which mimics a multi-field exposure. Fig. 11, we present the dose-volume histograms obtained in healthy brain (left) and scalp (right) as a function of H, for N = 375 voxels and Ptum = 4.5 cm. We note the cleanliness of the treatments obtained for H ≥ 30 cm. For H = 40 cm e.g., the doses exceed the most conservative peak doses, 11 Gy-Eq (brain), 30 Gy-Eq (skull), 16.7 Gy-Eq (scalp), only in 0.41% of the volume of the patient's head, and this by little, +1-2% max. We also give fig. 12 the evolution of FPEN (left) and Dtum (right) with H. The corresponding treatment times T, necessary to reach 30 Gy-Eq in the tumor with 30 mA of beam intensity, are also indicated in fig. 12 (right). As one would expect, the more the thickness of the moderator increases, the better the neutrons are moderated and guided, the cleaner the treatment. But at the same time, the more H increases, the lower the neutron flux at the exit of the moderator, the longer the treatment. The choice of the thickness of the moderator will therefore require a compromise between the therapeutic objective (minimization of exposure to healthy tissues) and the economic objective (minimization of treatment time).

Evolution of the structure and efficiency of a TOPOPT D2O moderator with tumor depth

In this section, we will study how the optimal shape of the D2O moderator and the associated treatment cleanliness evolve with the depth Ptum of the tumor. For this study, which is costly in computing power, we reduce the number N of voxels of the moderator to 150. As the volume of moderator voxels is subsequently increased, the MCNP simulations can be run with a lower number of source neutrons, here 2x10 9 . This reduction in the spatial resolution makes it possible to reduce the computation time per Ptum value to ~15 days, instead of 2 months beforehand. The considered Ptum depths vary from 3.5 cm (tumor in contact with the skull) to 8.5 cm (tumor at the center of the head). The thickness H of the moderator is set at 30 cm. The TOPOPT configurations obtained are presented in fig. 13-14 (left). For each configuration, the total doses Di and the dose excesses Ei in the voxels of the head are presented fig. 13-14 (center and right). For Ptum ≤ 5.5 cm, we observe fig. 13 that the algorithm recreates a moderation body, a reflector, a neutron guide and a central needle, as discussed in section 3.2. At 5.5 cm, however, a topology transition occurs; the moderator shape evolves to be closer to that of a moderated gun, whose moderation volume decreases as Ptum increases. Such a development is consistent with the fact that the deeper the tumor, the harder the neutron energy spectrum must be to reach it. For Ptum ≤ 5.5 cm, the treatment maintains a good degree of cleanliness, with doses exceeding peak values only in a small volume of the patient's head, and by little. Beyond that, the quality of treatment deteriorates. This development is made clearly visible in fig. 15, where we present the dose-volume histograms obtained for healthy brain (left) and scalp (right) for each depth. For tumors deeper than 5.5 cm, it will probably be necessary to modify the thickness H and the radius (so far set at 20 cm) of the moderator, as suggested by the perhaps too small thickness of the reflectors visible in fig. 14. Materials other than heavy water will also have to be tested, which may perform better at very high depths. This specific study will be conducted as part of a project currently being submitted. We give fig. 16 for information the evolution of the Dtum dose rate with Ptum. The corresponding treatment times T, necessary to reach 30 Gy-Eq in the tumor with 30 mA of beam intensity [START_REF] Capoulat | The 9 Be(d,n) 10 B reaction as a neutron source for Boron Neutron Capture Therapy[END_REF], are indicated.

Manufacture of a TOPOPT D2O moderator

The densities ρj of heavy water obtained by topology optimization in the voxels Θj of the moderator vary between 0 and ρmax, the density of heavy water. Strictly compliant manufacture of a TOPOPT moderator will therefore require the use of porous materials, e.g. beads of varying sizes, or alveolar media manufactured by 3D printing. Such technologies exist, cf. e.g. [20,[START_REF] Icesl Slicer | Progressive infils[END_REF]. However, a simplification of the design, of the binary type, where we force each moderator voxel to contain either 0% or 100% heavy water, would, at least initially, be easier to machine. Fig. 17 (left) presents the design obtained for H = 30 cm and Ptum = 4.5 cm by taking ρj = ρmax if ρj > ρmax/2, ρj = 0 otherwise. The FPEN value obtained for this simplified configuration, called TOPOPT Black & White (BW) in the rest of this study, is equal to 0.1998 ± 0.0062, against 0.1892 ± 0.0058 for the original TOPOPT configuration. These two values are statistically compatible, the simplification of the design does not lead to a significant deterioration in the quality of treatment. This result is consistent, since apart from a few border voxels, the TOPOPT design shown in fig. 7 is already predominantly binary. This simplification of the design facilitates the manufacturing of the TOPOPT moderator. One possible solution is to machine, e.g. from aluminum or steel, a cylindrical tank 20 cm in radius and 30 cm in length. The neutron guide and the central needle of the TOPOPT BW moderator can then be 3D printed and positioned in the tank, which is then filled with heavy water (32 liters for the design in fig. 17). We give fig. 17 (right) a photo of a section of the neutron guide, printed in PLA using a RAISE3D machine from the CNRS LPSC printing platform. There is a small hole, made for the heavy water to fill the ring in the neutron guide visible in fig. 7

Robustness of the design of a moderator

In this section, we present a set of sensitivity studies, useful for evaluating the robustness of the TOPOPT calculations. We have explored the sensitivity of the quality of treatment to the details of the structure of a TOPOPT moderator, to the voxelization of the patient's head, to the parameterization of the FPEN function, to the uncertainties on transport data and to the morphological diversity of the patients.

Robustness of the topological optimization

Global or local optimum. The moderator shapes computed in section 3 are complex but proved efficient. One can nevertheless legitimately wonder if each of these shapes really constitutes a global optimum of the optimization problem (2), and not a local optimum. Could there be even more efficient shapes? The TOPOPT algorithm is a complex iterative algorithm, and predicting the outcome of an iterative algorithm, even a very simple one, is in the vast majority of cases impossible, cf. Langton's ant or Syracuse conjecture. Nonetheless, there are arguments in favor of the global optimum. In Ref. [18], we verified, for several optimization problems constructed to have a solution obvious to a human or an analytical solution that the TOPOPT algorithm does indeed reach the global optimum. In addition, it is possible to test the robustness of the solutions found, e.g. the solution proposed in fig. 7, by modifying some geometry details, for example by removing the D2O ring or the central needle described in section 3.2, then by looking at the impact of these modifications on the treatment quality. The geometries of these modified moderators are presented in fig. 18. The FPEN values obtained are 0.2030 ± 0.0063 (no needle) and 0.2202 ± 0.0067 (no ring), against 0.1892 ± 0.0058 for the TOPOPT moderator. We see a degradation in the cleanliness of the treatment, indicating that the intricacies of the structure of a TOPOPT moderator (at least this one) are not superfluous. Existence of dose hot spots. The head model used in this study, presented in section 2 and fig. 1, is composed of M = 89 voxels: 1 for the tumor, 88 for the healthy tissues. Even if its axial symmetry eliminates the need for a discretization in angles around the axis, one can legitimately wonder if this small number of voxels can mask, by averaging it on the volume of a voxel, a (or several) local deposit of very high dose, of the hot spot type, in particular in tissues located near the exits of the moderators. To study this possibility, we took the TOPOPT moderator shown fig. 7, and recalculated with it the doses deposited in a head model comprising M = 881 voxels, 1 for the tumor, 880 for the healthy tissues, drawn in fig. 19 (left). The resulting map of dose excesses is given fig. 19 (left). By comparing it with that obtained in fig. 4 for M = 89, we see no hot spot, only a slight increase in the dose delivered in the scalp in front of the central needle, cf. fig. 7. To complete the analysis, we plot fig. 19 (right) the dose-volume histograms in healthy brain and scalp obtained for M = 89 and 881. We again note a slightly higher dose deposition in an area encompassing 0.04% of the volume of the scalp, located in front of the central needle. We therefore retained M = 89 in this study, a value which ensures a compromise between the spatial resolution of the dose deposition and the large computation time (as a reminder, 2 months for a TOPOPT computation with N = 375 and M = 89, against ~20 months probably with M = 881). Sensitivity to optimization parameters. The TOPOPT moderators presented in section 3 were obtained by minimizing the FPEN function formulated in Eq. ( 2). This formulation of FPEN is a human choice, reasoned, but which constitutes nonetheless a weak point. One can legitimately wonder if the shapes and performance of moderators could be significantly improved with a different formulation of FPEN. Here again the question of the global character of the calculated optima arises. Ideally, to answer this question, it would be necessary to propose a large number (~ 100) of different yet credible FPEN formulas, to compute for each of them the shape and the properties of the TOPOPT moderator, then to conclude on the robustness of the design to the choice of the formulation of the optimization problem. Such a study would however require a gigantic computing power, which would amount in decades per 24 CPUs server. For the time being, we therefore propose a less ambitious but more economical study, which consists in keeping the formulation (2) of FPEN, but by modifying the only parameter which is arbitrarily chosen there, the number n. Indeed, the only condition on n is that it is ≫ 1. To do this, we recalculated the shape of the optimal moderator for H = 30 cm, N = 150 voxels, Ptum (right) corresponding map of the dose excesses, Ei, in percent. Compared to the design shown fig. 13 for n = 30, we note an evolution in the shape of the moderator, with a reduction in the length of the neutron guide and the central needle, both in favor of a larger moderation volume. However, these design changes do not lead to a significant change in the quality of treatment, as can be seen by comparing the maps of the dose excesses fig. 13 for n = 30 and fig. 20 for n = 200. By recalculating FPEN with n = 30 for the design shown fig. 20, we get FPEN = 0.2367 ± 0.0124, against 0.2245 ± 0.0112 for the design shown fig. 13, two values compatible in their error bars. The treatment quality therefore seems not very sensitive to the choice of the number n, which constitutes a first element of response on the robustness of the optimization to the formulation of FPEN. However, we observe that the Dtum dose rate obtained for the design shown fig. 20 is equal to 0.9468 ± 0.0016 fGy-Eq/s.n., against 1.0495 ± 0.0017 fGy-Eq/s.n. for the design shown fig. 13. The treatment time for the TOPOPT design with n = 200 is therefore increased by 10% compared to the TOPOPT design with n = 30. Taking a value of n that is too large is therefore counterproductive. Indeed, for very large n values, function FPEN varies considerably during convergence, over 70 orders of magnitude for n = 200: FPEN = (2.6 ± 1.4)×10 74 at NUMITER = 0, (9.7 ± 3.9)×10 51 at NUMITER = 10, (4.0 ± 1.8)×10 3 at convergence. Above all, its values are affected by very large statistical errors, since they grow proportionally with n (derivative of FPEN). Such variations and stat. uncertainties alter the convergence of the TOPOPT algorithm, by penalizing in a counterproductive way the transfers of matter between the voxels of the moderator which would make it possible to exceed, even slightly, the peak doses at a place of the patient's head, to significantly reduce exposure elsewhere. We therefore kept n = 30 throughout the study, a reasonable value which satisfies the condition n ≫ 1 while limiting the statistical error on FPEN at a few %. 4.2. Sensitivity to nuclear data We will address here the question of the uncertainty on the calculated quality of treatment induced by the uncertainties on the transport data. These transport data are produced through models whose use or parameters can differ from one evaluator to another. They are affected by statistical and systematic uncertainties, due to imperfections or unknowns in the experimental devices used to measure them, to statistical fluctuations in measurements, to approximations in modeling and analysis, etc. Carrying out such a sensitivity analysis & uncertainty quantification study is a considerable work, which goes beyond the scope of this study. To address this nonetheless important question, we propose here a simplified and common approach, which consists in performing the calculations using several different databases, then in comparing the results obtained. For this study, we used the moderator shown fig. 7, obtained for N = 375 voxels and Ptum = 4.5 cm, using transport data from the ENDF/B-VII.0 database (/B-VI.6 for 1 H and 138 Ba), cf. section 2. We then recalculated the doses deposited in the patient's head, the FPEN function, and the Dtum dose rate in fGy-Eq/s.n., using two other databases, JEFF-3.1 and JENDL-3.3. This study required some adjustments: (i) as JENDL-3.3 lacks the gamma-ray production cross-sections for 2 H, 35 Cl, 37 Cl and 138 Ba (.42c), we used for these isotopes data from ENDF/B-VII.0 & /B-VI.6 (.70c for 2 H, 35 Cl, 37 Cl, .66c for 138 Ba); (ii) same for 138 Ba in JEFF-3.1 (.03c), replaced by ENDF/B-VI.6 (.66c); (iii) The S(α,β) thermal data were all taken from ENDF (hwtr.10t for the moderator, grph.10t and al27.12t for the source casing, lwtr.10t for the head). Fig. 21, we present, as a function of the voxel number i, the relative deviations, Di(x)/Di(ENDF) -1), with x = JENDL-3.3 or JEFF-3.1, between the 10 B, 14 N(n,p), neutron (without 14 N(n,p)), γ and total doses, calculated with ENDF/B-VII.0 & /B-VI.6, JEFF-3.1 and JENDL-3.3 (note: these are doses expressed in Gy-Eq, obtained by setting the total dose Di=0 deposited in the tumor at 30 Gy-Eq, see fig. 21). We observe that, whatever the voxel i > 0 considered, the JEFF-3.1 and JENDL-3.3 total doses are systematically lower than the ENDF total doses, sometimes significantly (note: with the exception of voxels 80 and 88 for JEFF-3.1, where the deviations appear positive but are in fact consistent with 0 in the statistical error bars). It follows that the treatments predicted with JEFF-3.1 and JENDL-3.3 are cleaner than that predicted with ENDF, but also a little longer, as shown by the values of FPEN and Dtum calculated for these 3 databases given table 3. It is probable that JENDL-3.3 underestimates the exposure of healthy tissues, because its data, particularly of gamma-ray production, are less complete. The difference between JENDL-3.3 and ENDF gamma doses is clearly visible in fig. 21 (left). Throughout our study, we adopted a conservative approach, consistently using the most unfavorable database, ENDF/B-VI.6 & VII.0, as it predicts a higher exposure to healthy tissues than JENDL-3.3 or JEFF-3.1.

Sensitivity to morphological parameters

The TOPOPT moderators obtained in the previous sections were computed using a simplified model of the patient's head, presented and discussed section 2. The shape and composition of a real human head are much more complex, so one can legitimately wonder what impact this simplification, common to any analytical or even voxelized model, can have on the quality of treatment. In addition, the morphology of a human head varies considerably from patient to patient. It is therefore interesting to study how (1) the modeling of the head and (2) the great morphological variability of patients can impact the design of the moderator and its associated treatment quality. The TOPOPT procedure presented in this study may provide an element of solution to these problems, thanks to its adaptability. In this section, for the time being, we propose to provide an element of response to the problem (2).

To answer this, we generated 92 different patients' heads, by independently sampling for each patient its brain radius Rb, its skull radius Rs and its head radius Rh, according to Gaussians of centers <Rb> = 7.05 cm, <Rs> = 8.21 cm and <Rh> = 8.71 cm, and 5.35% standard deviation. These values comply with Ref. [START_REF] Bushby | Centiles for adult head circumference[END_REF]10], which give size distributions of adult heads. However, since a variable distributed according to a Gaussian can take any value between -∞ and +∞, we have prohibited values of Rh less than 0.899×<Rh> (1st percentile) and greater than 1.115×<Rh> (99th percentile) [10]. In the absence of data found on the minimum and maximum thicknesses of skull and scalp, Rb and Rs radii are sampled freely, provided that Rb < Rs < Rh. These 92 generated head configurations all have different radii, and different brain, skull or scalp volumes. For each of them, we have adapted the voxelization shown in fig. 1: 1 voxel for the CTV (i = 0), 64 for healthy brain tissue (i = 1…64), 16 for the skull (i = 65…80), 8 for the scalp (i = 81…88), 89 in total. Each patient is positioned at the exit of the TOPOPT BW moderator shown fig. 17. The positions of the centers of the heads are adjusted so that the minimum distance between the surface of the head and the exit of the moderator remains equal to 1 mm, cf. section 2. The volume of the CTV is equal to 20 cm 3 , its depth Ptum equal to 4.5 cm. Fig. 22, we plot the distributions of the dose-volume histograms obtained in healthy brain and scalp for the 92 patients, indicating the 1st, 5th, 25th, 50th, 75th, 95th and 99th percentiles. We find that the doses deposited in the brain exceed the peak dose used by Capoulat et al. [START_REF] Capoulat | The 9 Be(d,n) 10 B reaction as a neutron source for Boron Neutron Capture Therapy[END_REF], 11 Gy-Eq, only in a small fraction of the brain for all patients. Doses exceed the peak dose used by Torres-Sánchez et al. [3], 12.5 Gy-Eq, in 25 out of 92 patients. However, for 23 of these 25 patients, the affected brain volume is negligible, less (or even much less) than 0.2%. For 2 out of 92 patients, who thus stand out exceptionally, the volume of brain exposed to a dose higher than 12.5 Gy-Eq reaches 4.7% for one, 4% for the other. These 2 patients have unusual morphological parameters: (i) for the first, a small head, Rh = 7.99 cm (the 2nd smallest in the sample), and thin thicknesses of skull and scalp, resp. 4.1 mm and 0.5 mm; (ii) for the second, an abnormal skull thickness, 1.3 mm, and radiuses of brain, skull and head, 7.60 cm, 7.73 cm and 8.30 cm, all quite far from the average values. The findings are similar for scalp. The doses deposited in the scalp never exceed the peak dose used by Torres-Sánchez et al. [3], 24 Gy-Eq. They remain below 18 Gy-Eq in the bulk of the scalp volume for all patients. These results are important, and rarely presented in the literature. They show that, for a tumor of intermediate depth, here 4.5 cm, reusing a moderator designed for a generic model patient, whose head has the standard proportions, will most often have no major harmful impact on the healthy tissues of a random patient, whatever its morphology of his head. However, for a fraction of patients, of the order of a few percent, side effects are to be feared, either in the brain (somnolence syndrome) or in the scalp.

Conclusion and perspectives

In this study, we applied a topology optimization algorithm developed at the CNRS LPSC to compute the structure of a heavy water neutron moderator for an AB-BNCT unit, for various tumor depths or moderator thicknesses. Calculations performed for deep glioblastomas using a 9 Be(d(1.45 MeV), n) 10 B source give convincing results. Despite the high energy, ~2 MeV in average, of the source neutrons, the TOPOPT moderators manage, in a compact volume, less than 0.04 m 3 , to deliver targeted doses, which reach 30 Gy-Eq in the tumor in a reasonable time while sparing healthy tissues. The local doses deposited in these tissues remain below the recommended limits, in almost the entire head volume, for tumor depths of up to 6 cm at this time. TOPOPT moderators have sophisticated structures, inaccessible in their detail to intuition or to previous parametric methods. They contain unexpected components, e.g. neutron guides mimicking multi-field treatments, the efficiency of which could inspire other designs. The progress made in manufacturing processes, 3D printing for example, now makes the machining of such components accessible to a small structure, laboratory or hospital unit. The versatility of the TOPOPT approach makes it possible to automatically fit the design of the moderator to the configuration of the BNCT treatment unit considered and to the biological parameters of the patients. However, despite these first promising results, progress remains to be made, on the computation time, on the head model, or on the maximum tumor depth that can be properly reached. As such, we plan to study, as part of a submitted project, materials other than heavy water that could be of therapeutic and economic interest. Heavy water is indeed an efficient material for slowing down neutrons in a compact volume without capturing them, but it remains expensive and is a special nuclear material.

Acknowledgments

We would like to thank Véronique Ghetta, CNRS LPSC, for her advice and the driving role she played in initiating this work.

Figure 1 .

 1 Figure 1. (Left) simplified models of the neutron source (yellow), the biological protection (blue), the concrete wall (green), the moderator and its voxelization (gray), the patient's head and the CTV (red) ; (right) simplified model of the patient's head, with its voxelization: CTV (red), healthy brain (orange), skull (green), scalp (blue).

Figure 2 .

 2 Figure 2. Neutron kerma Kn(E) in Gy.cm 2 /neutron(n) in the brain, the skull and the scalp; boron kerma, in Gy.cm 2 /neutron(n), per ppm of 10 B in a tissue.

Figure 3 .

 3 Figure 3. (Left) comparison of a reconstituted photon kerma in healthy brain with data of Ref. [15]; (right) photon kerma, Kγ(E), in Gy.cm 2 /gamma(g), for brain, skull and scalp.

  [START_REF] Capoulat | The 9 Be(d,n) 10 B reaction as a neutron source for Boron Neutron Capture Therapy[END_REF]. The values L = 16.7 and 24 Gy-Eq respectively used by Capoulat et al. and Torres-Sánchez et al. for the scalp are consistent with the findings of a study by Menéndez et al.

18 3 . 1 .

 31 Formalization and resolution of the optimization problem 19 20

  2) 28 29 The problem (2) consists in finding the densities ρ = (ρ1, …, ρN) in D2O of the moderator voxels 30 which minimize a penalization function FPEN, while respecting 4 constraints: (C1) the neutron 31 transport must obey the Boltzmann equation, Bn(ρ)ϕn = Qn, where Bn is the Boltzmann operator, 32 a function of the physicochemical characteristics of the materials crossed, therefore of the 33 densities ρ. Function ϕn(r,E,Ω,ρ) is the angular fluence of neutrons at position r = (x, y, z), of 34 energies E and directions Ω. The Qn term is the distribution in positions, energies and directions 35 of the neutrons generated by the source, as a reminder 9 Be(d(1.45 MeV), n) for this study; (C2) 36 photon transport must also obey a Boltzmann equation, with a source term Qγ depending on the 37 neutron fluence (production of gammas by neutron captures); (C3) the dose D0 deposited in the 38

Figure 4 .

 4 Figure 4. (Left) D2O moderator configurations obtained for N = 375 voxels, Ptum = 4.5 cm, H = 30 cm and NUMITER = 0, 20, 50 and 95; (center) corresponding maps of the total doses, Di, in Gy-Eq deposited in the patient's head; (right) corresponding maps of the dose excesses, Ei, in percent.

Figure 5 .

 5 Figure 5. Evolutions of FPEN (left) and Dtum (right) with NUMITER, for Ptum = 4.5 cm, N = 375 and H = 30 cm.

Figure 6 .

 6 Figure 6. Dose-volume histograms in healthy brain (left) and scalp (right), obtained for the TOPOPT moderator calculated for Ptum = 4.5 cm, N = 375 and H = 30 cm.

Figure 7 .

 7 Figure 7. Details of the TOPOPT moderator calculated for N = 375, Ptum = 4.5 cm and H = 30 cm.

Figure 8 .

 8 Figure 8. Relative contributions, in %, to the total dose Di of each of the 4 dose components presented section 2, wBDB, wnDn (all reactions except 14 N(n,p)), wpDp (only 14 N(n,p)) and wγDγ, calculated for the TOPOPT moderator obtained for N = 375 voxels, Ptum = 4.5 cm and H = 30 cm (convergence at iteration 95) and shown fig. 7.

Figure 9 .

 9 Figure 9. Energy spectra at the exit of the TOPOPT moderator obtained for N = 375, Ptum = 4.5 cm and H = 30 cm, at positions 1 to 4 indicated on the left.

Figure 10 .

 10 Figure 10. (Left) configurations of the TOPTOP D2O moderators obtained for N = 375 voxels, Ptum = 4.5 cm and H = 20, 25, 35 and 40 cm; (center) corresponding maps of the total doses, Di, in Gy-Eq deposited in the patient's head; (right) corresponding maps of the dose excesses, Ei, in percent.

Figure 11 .

 11 Figure 11. Dose-volume histograms in healthy brain (left) and scalp (right), obtained for the TOPOPT moderator calculated for Ptum = 4.5 cm, N = 375 voxels and H = 20, 25, 30, 35, 40 cm.

Figure 12 .

 12 Figure 12. Evolution with H of FPEN (left), Dtum (right) and T (right) at convergence, for Ptum = 4.5 cm and N = 375 voxels.

Figure 13 .

 13 Figure 13. (Left) configurations of the TOPOPT D2O moderators obtained for N = 150 voxels, H = 30 cm and Ptum = 3.5, 4.5 and 5.5 cm; (center) corresponding maps of the total doses, Di, in Gy-Eq deposited in the patient's head; (right) corresponding maps of the dose excesses, Ei, in percent.

Figure 14 .

 14 Figure 14. (Left) configurations of the TOPOPT D2O moderators obtained for N = 150 voxels, H = 30 cm and Ptum = 6.5, 7.5 and 8.5 cm; (center) corresponding maps of the total doses, Di, in Gy-Eq deposited in the patient's head; (right) corresponding maps of the dose excesses, Ei, in percent.

Figure 15 .

 15 Figure 15. Dose-volume histograms in healthy brain (left) and scalp (right), obtainded for N = 150 voxels, H = 30 cm and Ptum ranging from 3.5 to 8.5 cm.

Figure 16 .

 16 Figure 16. Evolution of Dtum and T with Ptum, at convergence, for H = 30 cm and N = 150 voxels.

  or 17.

Figure 17 .

 17 Figure 17. (Left) BW design of the TOPOPT moderator calculated for N = 375 voxels, Ptum = 4.5 cm and H = 30 cm; (right) 3D printing in PLA of a section of the neutron guide of this moderator.

Figure 18 .

 18 Figure 18. (Left) design obtained by removing the central needle in the TOPOPT moderator shown fig. 7; (right) design obtained by removing the D2O ring in the neutron guide of the TOPOPT moderator shown fig. 7.

Figure 19 .

 19 Figure 19. (Left) map of the dose excesses Ei, in %, obtained for the head model with M = 881 voxels and the TOPOPT moderator shown fig. 7; (right) dose-volume histograms in healthy brain and scalp, obtained for M = 89 and 881 voxels.

  = 4.5 cm, this time taking n = 200 instead of 30 so far. The configuration of the TOPOPT moderator thus obtained is shown fig. 20 (left); the corresponding maps of doses and dose excesses in the head voxels are presented in fig. 20 (center and right).

Figure 20 .

 20 Figure 20. (Left) configuration of the TOPOPT D2O moderator obtained for N = 150 voxels, H = 30 cm, Ptum = 4.5 cm and n = 200; (center) corresponding map of the total doses, Di, in Gy-Eq deposited in the patient's head;

Figure 21 .Table 3 .

 213 Figure 21. Deviations in % between the 10 B, 14 N(n,p), neutron (w/o 14 N(n,p)), γ and total doses, calculated for the TOPOPT moderator of fig. 7 with the ENDF/B-VI.6 & VII.0 and JENDL-3.3 databases (left), and with the ENDF/B-VI.6 & VII.0 and JEFF-3.1 databases (right). Note: the tumor voxel is the voxel i = 0, brain is i = 1…64, skull is i = 65…80, scalp is i = 81…88.

Figure 22 .

 22 Figure 22. Distributions of the dose-volume histograms in healthy brain (left) and scalp (right) obtained for the 92 patients, using the TOPOPT BW moderator shown fig. 17 for Ptum = 4.5 cm and H = 30 cm.

Table 1 .

 1 Chemical compositions (in mg/g) and densities (in g/cm3 ) in head tissues[11].

		Brain Skull Scalp
	H	107	50	100
	C	145	212	204
	N	22	40	42
	O	712	435	645
	Na	2	1	2
	Mg	0	2	0
	P	4	81	1
	S	2	3	2
	Cl	3	0	3
	K	3	0	1
	Ca	0	176	0
	Density (g/cm 3 )	1.04	1.61 1.09
	Tissue wB wn wp wγ CB (µg/g)
	tumor 3.8 3.2 3.2 1	52.5
	brain 1.3 3.2 3.2 1	15.0
	skull 1.3 3.2 3.2 1	15.0
	scalp 2.5 3.2 3.2 1	22.5

Table 2 .

 2 Biological weighting factors w and 10 B concentrations CB (in µg/g) in head tissues.


Gildas Morvan

Daniel Jolly

Multi-level agent-based modeling with the Influence Reaction principle

Keywords: multi-level simulation, influence reaction model, micro-macro link, cellular automata

This paper deals with the specification and the implementation of multi-level agent-based models, using a formal model, IRM4MLS (an Influence Reaction Model for Multi-Level Simulation), based on the Influence Reaction principle. Proposed examples illustrate forms of top-down control in (multi-level) multi-agent based-simulations.

Introduction

Four main theoretical issues emerge in the literature on multi-level1 agent-based modeling: the conception of a meta-model allowing a non ambiguous characterization of a multi-level model at the conceptual level [START_REF] Maus | Hierarchical modeling for computational biology[END_REF][START_REF] Morvan | IRM4MLS: the influence reaction model for multi-level simulation[END_REF][START_REF] Picault | An interaction-oriented model for multi-scale simulation[END_REF][START_REF] Uhrmacher | Combining micro and macro-modeling in devs for computational biology[END_REF], the introduction of a dynamic level of detail [START_REF] Navarro | Dynamic level of detail for large scale agent-based urban simulations[END_REF][START_REF] Soyez | Multi-agent multi-level modeling -a methodology to simulate complex systems[END_REF], the detection and reification of emergent phenomena [START_REF] Caillou | Simanalyzer : Automated description of groups dynamics in agent-based simulations[END_REF][START_REF] Chen | Identifying multi-level emergent behaviors in agent-directed simulations using complex event type specifications[END_REF][START_REF] David | See emergence as a metaknowledge. a way to reify emergent phenomena in multiagent simulations?[END_REF][START_REF] Moncion | Automatic characterization of emergent phenomena in complex systems[END_REF][START_REF] Prévost | Detection and reification of emerging dynamical ecosystems from interaction networks[END_REF][START_REF] Vo | A modelling language to represent and specify emerging structures in agent-based model[END_REF] and the representation of aggregated entities [START_REF] Parunak | Between agents and mean fields[END_REF]. This paper focuses on the first one, with respect to the Influence Reaction (IR) principle, shortly, action as a two step process: (1) agents produce "influences", i.e., individual decisions, according to their internal state and perceptions (2) the system "reacts", i.e., computes the consequences of influences, according to the state of the world [START_REF] Ferber | Influences and reaction: a model of situated multiagent systems[END_REF]. This model has been extended in several ways, and notably for multi-agent based-simulation (MABS), by adding an explicit representation of time [START_REF] Michel | The IRM4S model: the influence/reaction principle for multiagent based simulation[END_REF]. An IR-based meta-model, IRM4MLS (an Influence Reaction Model for Multi-Level Simulation), and its Java implementation are introduced in the section 2. Using, and then extending a simple example: the Conway's game of life, two toy-models of increasing complexity are presented in the section 3, illustrating forms of top-down control in (multi-level) MABS. Results are discussed in the section 4.

IRM4MLS: an Influence Reaction Model for Multi-Level Simulation

In this section, IRM4MLS, an Influence Reaction Model for Multi-Level Simulation, is introduced [START_REF] Morvan | IRM4MLS: the influence reaction model for multi-level simulation[END_REF]. It extends IRM4S (an Influence Reaction Model for Simulation) in order to deal with multi-level models [START_REF] Michel | The IRM4S model: the influence/reaction principle for multiagent based simulation[END_REF]. From a technical perspective, levels can be viewed as interacting IRM4S simulations 2 .

A multi-level model is defined by a set of levels L and a specification of the relations between levels. Two types of relations are specified: an influence relation (agents in a level l are able to produce influences in a level l = l) and a perception relation (agents in a level l are able to perceive the dynamic state of a level l = l), represented by directed graphs denoted respectively < L, EI > and < L, EP >, where EI and EP are two sets of edges, i.e., ordered pairs of elements of L. Influence and perception relations in a level are systematic and thus not specified in EI and EP (cf. eq. 1 and 2).

The in and out neighborhood in < L, EI > (respectively < L, EP >) are denoted N - I and N + I (resp. N - P and N + P ) and are defined as follows:

∀l ∈ L, N - I (l) (resp. N - P (l)) = {l} ∪ {l ∈ L : l l ∈ EI (resp. EP )}, ( 1
)
∀l ∈ L, N + I (l) (resp. N + P (l)) = {l} ∪ {l ∈ L : ll ∈ EI (resp. EP )}, (2)
The set of agents in the system at time t is denoted A(t). ∀l ∈ L, the set of agents belonging to l at t is denoted A l (t) ⊆ A(t). An agent belongs to a level iff a subset of its physical state φa belongs to the state of the level. Thus, an agent belongs to zero, one, or more levels. An environment models the natural dynamics of level properties and can be shared by different levels (fig. 1).

The dynamic state of a level l ∈ L at time t, denoted δ l (t) ∈ ∆ l , is a tuple < σ l (t), γ l (t) >, where σ l (t) ∈ Σ l and γ l (t) ∈ Γ l are the sets of environmental properties and influences of l. The behavior of an agent a ∈ A l is defined as Decision l a • M emorizationa • P erception l a , with P erception l a :

l P ∈N + P (l)
∆ l P → l P ∈N + P (l)

P l P a , (3)
M emorizationa :

l∈L|a∈A l l P ∈N + P (l)

P l P a × Sa → Sa, (4)
Decision l a : Sa →

l I ∈N + I (l) Γ l I . (5)
There is no memorization function specific to a level to preserve the coherence of the internal state of the agents. The environment ω of a level l produces influences through a function:

N atural l ω : ∆ l → l I ∈N + I (l) Γ l I . (6)
The reaction function computes next level state and time advance:

Reaction l : Σ l × Γ l → ∆ l × T l . (7)
The time representation is inspired by DEVS (Discrete EVent System specification) [START_REF] Zeigler | Theory of Modeling and Simulation[END_REF]. T = l∈L {T l } denotes the time vector of the simulation, such as ∀l ∈ L, T l =< t l , dt l >, where t l represents when the current event (or step, depending on the simulation model) time and dt l its lifespan. The final simulation time is denoted t f . The algorithm 1 ensures the scheduling of these different functions with respect to temporal constraints of perception and memorization, influence production and reaction [START_REF] Morvan | IRM4MLS: the influence reaction model for multi-level simulation[END_REF].

The implementation of IRM4MLS is based on the idea of micro kernel, taken from MadKit3 [START_REF] Gutknecht | The madkit agent platform architecture[END_REF]. Thus, in this approach, a technical agent, e.g., an observer or a message broker, would be scheduled with respect to the IR principle (cf. algo. 1), the concept of level ensuring a clear separation between system and simulation agents.

The API is minimal (seven high-level abstractions) and specifies only the methods needed to schedule a model (fig. 2). Most methods are generic and then are implemented at an abstract level. Basically, to implement a model one only has to override perception, memorization, influence production, reaction and initialization functions.

Agent and behaviors (such as environment and natural dynamics) are represented by different entities to clearly distinguish the core side of an agent, its state and memorization function, and the level sides (perception and influence production functions) that can change according to simulations. Algorithm 1: simulation model of IRM4MLS

Input: < L, E I , E P >, A, δ, T, t f Output: δ(t f ) 1 while ∃t l ≤ t f do 2 foreach a ∈ A(t) do 3 foreach l ∈ L : a ∈ A l (t) ∧ ∀l P ∈ N + P (l), t l ≥ t l P do 4 p a (t l ) = P erception l a (< δ l P (t l P ) : l P ∈ N + P (l) >); 5 end 6 s a (t l + dt l ) = M emorization a (p a ); 7 end 8 foreach l ∈ L : ∀l I ∈ N + I (l), t l ≤ t l I ∨ t l + dt l < t l I + dt l I do 9 foreach l I ∈ N + I (l) : t l ≤ t l I ∧ t l + dt l > t l I do 10 γ l I ω (t l I ) = N atural l ω (δ l (t l )) ; 11 foreach a ∈ A l (t) do 12 γ l I a (t l I ) = Decision l a (s a (t l + dt l )); 13 end 14 end 15 end 16 foreach l ∈ L : t l + dt l ∈ min(t + dt) do 17 γ l (t) = {γ l (t) l I ∈N - I (l) γ l I ω (t) a∈A l I γ l I a (t)}; 18 <δ l (t l + dt l ), T l >= Reaction l (σ l (t l ), γ l (t l )); 19 end 20 end
3 IRM4MLS in practice: multi-level games of life 3.1 Introducing a macroscopic parameter (top-down control)

In this section, a toy IRM4S (or 1-level IRM4MLS) model, is presented: a modified agent-based version of the Conway's game of life (or simply Life). This simple example illustrates where a macroscopic parameter, should (but should not) be introduced in an agent-based model that relies on the influence reaction principle: in the reaction function (but not in the behavioral functions of agents). Therefore, this parameter has a non-ambiguous semantics that does not depend on the updating scheme of the simulations, even in the case of strong interaction 4 . Each agent represents a cell that can be dead or alive and that has eight neighbors

Cells evolve in parallel: the reaction function can then simply be defined as "applying agent influences" (algo. 2). If a cell is dead and has three living neighbors or is alive and has two or three living neighbors, it will be alive at the next step; in other cases, it will be dead. Let specify the behavior of the agents:

1. they perceive the number of living cells in their neighborhood (algo. 3), 2. memorize their internal state, i.e., their next state (algo. 4), 3. and then, decide whether or not they will be alive at the next step according to their internal state (algo. 5).

The environment is static: there is no natural dynamics and thus, N aturalω returns ∅.

One cruel and ironic aspect of Life is that a cell has generally little chance to remain alive in the long run (fig. 3 is a board composed of small still lifes and 1-period oscillators. This behavior is predictable knowing the λ parameter of the game of life 5 . λ is a complexity measure of cellular automata introduced by [START_REF] Langton | Computation at the edge of chaos: Phase transitions and emergent computation[END_REF] that depends on the number of cell states K, the neighborhood N and the number n of transitions to a quiescent state Sq in the transition function such as

λ = 1 - n K N , (9)
with K = 2, N = 9 and n = 8 2 + 2 • i=0,1,4-8 8 i for Life. To improve λ lif e to a value λ + , one needs to change the rules. In this example, λ + is regarded as a macroscopic parameter, explicitly introduced in the model and independent from cell behaviors 6 : influences of dying cells are not taken into account by the reaction function with a probability p such as

n = 8 2 + (2 -p) • i=0,1,4-8 8 i . (10)
Thus, their is a simple linear relation between p and λ + : p = (λ + -λ lif e )/0.3359. For λ + ∈ [0.48, 0.6], the number of dying cells tends to decrease in time and large structures of vertical or horizontal rows eventually emerge, shaped by moving groups of switching state cells that seem to work at their boundaries, and eventually vanish when the board becomes a dense still life of density ≈ 0.5 (fig. 4). 5 λ lif e = 0.2734375. For λ ≈ 0.25, "structures of period 1 appear. Thus, there are now three different possible outcomes for the ultimate dynamics of the system, depending on the initial state. The dynamics may reach a homogeneous fixed point consisting entirely of state Sq, or it may reach a heterogeneous fixed point consisting mostly of cells in state Sq with a sprinkling of cells stuck in one of the other states, or it may settle down to periodic behavior" [15, p. 17]. 6 An other macroscopic parameter, the asynchrony, has been previously introduced in Life is such way [START_REF] Fatès | Does "life" resist asynchrony?[END_REF].

Top-down feedback control

In the previous example, while the macroscopic parameter λ + has an influence on agents, it is not related to the state of the system and therefore, there is no need to observe it. The reaction function can then be viewed as an open-loop controller. In the next example, a top-down feedback control is introduced.

The goal of the multi-scale model presented in this section is to keep Life boards at the desired density (ρ + ), by controlling the proportion of dying cells at the mesoscopic level, to account for the natural variability of density between regions of the grid (fig. 3(b)). Moreover, the control should affect as less as possible simulations at the microscopic level and, to keep it simple, should be tuned by a single linear parameter such as λ + in the previous model.

Two levels are considered: the cell (or microscopic) level, lm and the cell region (or mesoscopic) level, lM . At the mesoscopic level, the model properties are the expected density and the cells in each region:

∀t, σ l M (t) = {ρ + } a l M ∈A l M {cells(a l M )}. ( 11
)
The cells behave according to the game of life rules (algo. 3, 4 and 5). However, Reaction lm depends on mesoscopic influences (algorithm 6). EP and EI are equal to {lM lm}. Mesoscopic agents have a proportional control behavior. They 1. perceive the density of living cells in a region, 2. memorize their internal state, i.e., the difference between expected and actual densities, 3. and then decide the influence sent to agents of lm: The kP parameter has to be carefully tuned to run realistic simulations: too small simulations do not achieve the desired solution (ρ = ρ + ), too big the board density tends to oscillate around ρ + and the number of micro influences not taken into account by Reaction lm becomes too important. However, for appropriate kP values, this simple linear controller achieves good results and allows to find a good compromise between conflicting micro and meso knowledge (fig. 5).

∀a l M ∈ A l M ∀a lm ∈ cells(a l M ), command(a lm ) = kP • . (12

Conclusion 4.1 Discussion

An important issue of multi-level agent-based modeling, only briefly discussed here, is to define the adequate methodology. Indeed, the traditional MABS methodology is purely bottom-up: microscopic knowledge is used to construct models while macro- scopic knowledge is used to validate models [START_REF] Drogoul | Multi-agent based simulation: Where are the agents?[END_REF]. Thus, it seems irrelevant in a multilevel context. Three general conclusions can be drawn from the previous examples:

• a parameter should be introduced at its observation level in the model; therefore, each observed level should be explicitly represented in the model,

• an intuitive way to model an external control on a level l is to modify the reaction function of l, i.e., to modify the way influences of agents of l are taken into account (but not agent behavioral functions), according to external influences; an external feedback control implies both observation and influence relations:

N - I (l) = N - P (l),
• such a controller can be viewed as a technical tool that aims to find a compromise between conflicting knowledge from the different studied levels and achieve realistic simulations from different studied points of view [START_REF] Li | Multi-scale compromise and multi-level correlation in complex systems[END_REF].

Pattern oriented modeling (POM) consists in "the multi-criteria design, selection and calibration of models of complex systems" [START_REF] Grimm | Pattern-oriented modelling: a 'multi-scope' for predictive systems ecology[END_REF]. Many aspects of this methodology, developed in the context of ecology, seem particularly relevant for multi-level agentbased models since "patterns" are generally observed at different levels of organization in complex systems. However, the problem is far from being solved. For instance, the introduction of a dynamic level of detail raises several questions regarding, e.g., the validation of the model or the representation of composite agents [START_REF] Navarro | Dynamic level of detail for large scale agent-based urban simulations[END_REF][START_REF] Parunak | Between agents and mean fields[END_REF]. A case study of three real world multi-level agent-based models reveals other interesting methodological issues [START_REF] Gil-Quijano | From biological to urban cells: Lessons from three multilevel agent-based models[END_REF].

Related works

At least three works could be related to this one:

• ML-DEVS is an extension of DEVS that allows the simulation of multi-level models (and not only coupled models in which the behavior of a model is determined by the behaviors of its sub-models) [START_REF] Uhrmacher | Combining micro and macro-modeling in devs for computational biology[END_REF]. Two types of relation between levels are defined: information propagation and event activation which are quite similar to those defined in IRM4MLS. However, ML-DEVS supports only pure hierarchies of models, i.e., interaction graphs are viewed as trees [START_REF] Maus | Hierarchical modeling for computational biology[END_REF]. DEVS, as a generic event-based simulation framework, has also been extended to support agent-based models [START_REF] Müller | Towards a formal semantics of event-based multi-agent simulations[END_REF]. A major design difference between IR and DEVS based approaches is the technical orientation of the latter leading to an important gap between conceptual and computational models.

• PADAWAN (Pattern for Accurate Design of Agent Worlds in Agent Nests) is a multi-scale agent-based meta-model based on a compact matricial representation of interactions: IODA (Interaction-Oriented Design of Agent simulations) [START_REF] Kubera | Interaction-oriented agent simulations: From theory to implementation[END_REF][START_REF] Picault | An interaction-oriented model for multi-scale simulation[END_REF]. Moreover, authors analyze the structure of what is called here interaction graphs in multi-scale models (a relevant issue for IRM4MLS as well), and conclude they should be viewed as upper semilattices and not simply trees as suggested elsewhere. A major design difference between IR and IODA based approaches is that the latter constraints the definition of interactions, leading to a simple but restrictive simulation framework.

• GAMA 7 is a MABS platform with a dedicated modeling language, GAML, that offers multi-level capabilities [START_REF] Taillandier | GAMA: A simulation platform that integrates geographical information data, agent-based modeling and multi-scale control[END_REF]. Moreover, it includes a framework (a set of predefined GAML commands) to agentify emerging structures [START_REF] Vo | A modelling language to represent and specify emerging structures in agent-based model[END_REF]. It is certainly the most advanced platform, from an end-user point of view, that integrates a multi-level approach.

Perspectives

The main perspectives of this work concern the implementation of existing works with IRM4MLS:

• the concept of PolyAgent [START_REF] Parunak | Pheromones, probabilities and multiple futures[END_REF][START_REF] Parunak | Concurrent modeling of alternative worlds with polyagents[END_REF],

• multi-level organizational models widely used in engineering sciences such as holonic multi-agent systems (cf. footnote 2), system of systems and heterarchical control [START_REF] Cossentino | Aspecs: an agent-oriented software process for engineering complex systems[END_REF][START_REF] Maier | Architecting principles for system of systems[END_REF][START_REF] Morvan | Engineering hierarchical complex systems: an agent-based approach -the case of flexible manufacturing systems[END_REF],

• multi-scale tools: generic scaling operators and emergence detection and reification algorithms [START_REF] Chen | Identifying multi-level emergent behaviors in agent-directed simulations using complex event type specifications[END_REF][START_REF] Navarro | Dynamic level of detail for large scale agent-based urban simulations[END_REF].

Moreover, the first model presented in this paper could be used to explore the relations between computational capabilities of a cellular automaton (λ + ), noise (a function of λ + -λ) and entropy. Moreover, finding the conditions under which cells arrange themselves in a steady state could be an interesting way to solve heuristically large instances of the maximum density still life problem [START_REF] Elkies | The still-life density problem and its generalizations. Voronoi's Impact on Modern Science, Book I[END_REF].

The second model illustrates a form of simple proportional top-down feedback control. Such approach could be generalized to model more complex of cross-level feedback control, using integrations and derivates of observed variables.

Figure 1 :

 1 Figure 1: Main concepts of IRM4MLS (as a simplified class diagram)

Figure 2 :

 2 Figure 2: Java API of IRM4MLS

Figure 3 :

 3 Figure 3: Dynamics of the density of living cells, ρ, starting from a random grid with 100 replications (a) Macroscopic dynamics: expected density and variability of living cells in the whole grid. (b) Mesoscopic dynamics: expected density and variability of living cells in 10 × 10 cell clusters.

Figure 4 :Algorithm 6 : 5 else 6 alive

 4656 Figure 4: (a) Mean number of steps needed to converge to a steady state starting from a random grid (100 replications), simulations are stopped after 2 • 10 4 steps (simulations converge in the dark area). (b) Example of still life found for λ + = 0.5.

Figure 5 :

 5 Figure 5: Simulation examples with k P = 10•ρ + and initial density ρ(t 0 ) = 2•ρ + . r represents the rate of microscopic influences not taken into account (in %) by Reaction lm .

  (a)). Moreover, what you get most of the time,

	Algorithm 2: Reaction	Algorithm 4: M emorization a
		Input: σ(t), γ (t)		Input: p a (t), s a (t)
		Output: δ(t + 1)		Output: s a (t + dt)
	1 foreach a ∈ A(t) do		1 if alive(a) ∧ p a (t) ∈ {2, 3}∨
	2	alive(a) = γ a (t) ;		2 ¬alive(a) ∧ p a (t) ∈ {3} then
	3 end		3	s a (t + dt) = 1 ;
				4 else
				5	s a (t + dt) = 0 ;
	Algorithm 3: P erception a	6 end
		Input: δ(t)		Algorithm 5: Decision
		Output: p a (t)
	1 p a (t) =	alive(n) ;
		n∈neighbors(a)

a Input: s a (t + dt) Output: γ a (t) 1 γ a (t) = s a (t + dt) ;

The term multi-scale may also be found. Intuitively, a level and scale are similar concepts that both mean viewpoint. However, this notion should be clarified. In the following we assume that two agents are not at the same scale iff they represent processes that have different spatial and/or temporal extents. Two agents are not in the same level iff they cannot interact directly, i.e., with a single interaction function. It should follow from the previous definitions that it exists multi-scale models that are mono-level and conversely.

Therefore, each level has a microscopic side: the agent behaviors, and a macroscopic side: the reaction function. This aspect can also be found in holonic multi-agent systems[START_REF] Cossentino | Aspecs: an agent-oriented software process for engineering complex systems[END_REF].

http://www.madkit.org

Strong interaction implies that agents agree on the outcome of the interaction[START_REF] Michel | Weak interaction and strong interaction in agent based simulations[END_REF]. Thus, such model should not be simulated with a STRIPS-like meta-model, i.e., that views an action as a change of the state of the system. It would lead to problems of result replication[START_REF] Bigbee | Replication of sugarscape using mason[END_REF], but also of parameter and result interpretation[START_REF] Fatès | How important are updating schemes in multi-agent systems? an illustration on a multi-turmite model[END_REF][START_REF] Michel | Weak interaction and strong interaction in agent based simulations[END_REF].