Julien Bordes
email: julien.j.bordes@gmail.com

Sébastien Incerti
email: incerti@cenbg.in2p3.fr

Erick Mora-Ramirez
email: erick.mora@ucr.ac.cr

Jonathan Tranel
email: j.tranel@hotmail.fr

Cédric Rossi
email: cedric.rossi@chu-dijon.fr

Chris- Tine Bezombes
email: christine.bezombes@inserm.fr

Julie Bordenave
email: julie.bordenave@inserm.frmanuel.bardies@inserm.fr

Manuel Bardiès

Richard Brown
email: richard.brown@ucl.ac.uk

Marie-Claude Bordage
email: marie-claude.bordage@inserm.fr

Monte Carlo dosimetry of a realistic multicellular model of follicular lymphoma in a context of radioimmunotherapy

Keywords: Small-scale dosimetry, realistic biological data, rituximab, Auger electrons, β -particles

Purpose. Tumors geometry and radiopharmaceutical biodistribution impact the energy deposition in targeted radionuclide therapy. However, small-scale dosimetry studies generally consider an artificial environment where the tumors are spherical and the radionuclides are homogeneously biodistributed. To bring realism, we developed a dosimetric methodology based on a three-dimensional in vitro model of follicular lymphoma incubated with rituximab, an anti-CD20 monoclonal antibody used in the treatment of non-Hodgkin lymphomas, which might be combined with a radionuclide. The effects of the realistic geometry and biodistribution on the absorbed dose were highlighted by comparison with literature data. Additionally, to illustrate the possibilities of this methodology, the effect of different radionuclides on the absorbed dose distribution delivered to the in vitro tumor were compared.

Methods. The starting point was a model named Multicellular Aggregates of Lymphoma Cells (MALC). Three MALCs of different dimensions and their rituximab biodistribution were considered. Geometry, antibody location and concentration were extracted from selective plane illumination microscopy. Assuming antibody radiolabeling with Auger electron ( 125 I and 111 In) and β -particle emitters ( 177 Lu, 131 I and 90 Y), we simulated energy deposition in MALCs using two Monte Carlo codes: Geant4-DNA with "CPA100" physics models for Auger electrons emitters and Geant4 with "Livermore" physics models for β -particles emitters.

Results. MALCs had ellipsoid-like shapes with major radii, r, of ~ 0.25, ~ 0.5 and ~ 1.3 mm. Rituximab was concentrated in the periphery of the MALCs. The absorbed doses delivered by 177 Lu, 131 I and 90 Y in MALCs were compared with literature data for spheres with two types of homogeneous biodistributions (on the surface or throughout the volume). Compared to the MALCs, the mean absorbed doses delivered in spheres with surface biodistributions were between 18% and 38% lower, while with volume biodistribution they were between 15% and 29% higher. Regarding the radionuclides comparison, the relationship between MALC dimensions, rituximab biodistribution and energy released per decay impacted the absorbed doses. Despite releasing less energy, 125 I delivered a greater absorbed dose per decay than 111 In in the r ~ 0.25 mm MALC (6.78 • 10 -2 vs. 6.26 • 10 -2 µGy • Bq -1 • s -1 ). Similarly, the absorbed doses per decay in the r ~ 0.5 mm MALC for 177 Lu (2.41 • 10 -2 µGy • Bq -1 • s -1 ) and 131 I (2.46 • 10 -2 µGy • Bq -1 • s -1 ) are higher than for 90 Y (1.98 • 10 -2 µGy • Bq -1 • s -1 ). Furthermore, radionuclides releasing more energy per decay delivered absorbed dose more uniformly thorough the MALCs. Finally, when considering the radiopharmaceutical effective half-life, due to the biological half-life of rituximab being best matched by the physical half-life of 131 I and 177 Lu compared to 90 Y, the first two radionuclides delivered higher absorbed dose.

Conclusion.

In the MALCs considered, β -emitters delivered higher and more uniform absorbed dose than Auger electron emitters. When considering radiopharmaceutical half-lives, 131 I and 177 Lu delivered absorbed doses higher than 90 Y. In view of real irradiation of MALCs, such a work may be useful to select suited radionuclides and to help explain the biological effects.

Introduction

In targeted radionuclide therapy (TRT), a tumor-seeking agent is combined with an emitter of short-range radiation, such as Auger electrons (typically inferior to 1 µm in water), α particles (inferior to 100 µm) and β -particles (in the order of millimeters). The objective of this combination is to deliver in patients a lethal absorbed dose to tumor cells while minimizing the exposure of healthy tissues. Thus, the development and the evaluation of TRT requires small-scale dosimetry studies [START_REF] Strand | Radioimmunotherapy dosimetry-a review[END_REF][START_REF] Roeske | Small-scale dosimetry: challenges and future directions[END_REF] . However, previous studies have been generally based on simplifying assumptions (e.g. [START_REF] Bardiès | Absorbed doses for internal radiotherapy from 22 beta-emitting radionuclides: beta dosimetry of small spheres[END_REF][START_REF] Falzone | Monte Carlo evaluation of Auger electronemitting theragnostic radionuclides[END_REF][START_REF] Champion | Comparison between three promising β-emitting radionuclides, 67 Cu, 47 Sc and 161 Tb, with emphasis on doses delivered to minimal residual disease[END_REF][START_REF] Hindié | Dose deposits from 90 Y, 177 Lu, 111 In, and 161 Tb in micrometastases of various sizes: implications for radiopharmaceutical therapy[END_REF][START_REF] Amato | An analytical model to calculate absorbed fractions for internal dosimetry with alpha, beta and gamma emitters[END_REF][START_REF] Amato | Absorbed fractions for electrons in ellipsoidal volumes[END_REF] ) owing to the lack of data at the submillimeter scale regarding the tumor geometry and the radionuclide biodistribution.

Firstly, tumors were generally supposed spherical, or, alternatively, ellipsoidal. Absorbed dose in spheres was found to be higher than in the ellipsoids of the same volume [START_REF] Rahman | Influence of thyroid volume reduction on absorbed dose in 131 I therapy studied by using Geant4 Monte Carlo simulation[END_REF] , emphasizing the influence of the geometry. Secondly, small-scale dosimetry calculations have regularly assumed a homogeneous radionuclide biodistribution, either at the surface or within the volume of the tumor. In practice, the targeting of a tumor is not homogeneous, partly due to the heterogeneous distribution of the receptors (i.e. the antigens). Thus, mathematical models of heterogeneous radionuclide biodistributions within the tumors have been studied [START_REF] Bao | Theoretical study of the influence of a heterogeneous activity distribution on intratumoral absorbed dose distribution[END_REF][START_REF] Spaic | Generation of dose-volume histograms using Monte Carlo simulations on a multicellular model in radionuclide therapy[END_REF][START_REF] Uusijärvi | Tumour control probability (TCP) for non-uniform activity distribution in radionuclide therapy[END_REF][START_REF] Howell | Monte Carlo simulation of irradiation and killing in threedimensional cell-populations with lognormal cellular uptake of radioactivity[END_REF] . As could be expected, these models led to absorbed doses significantly less uniform within the tumor compared to homogeneous biodistributions. Therefore, given the impact of both tumor geometry and radionuclide biodistribution, there are clear calls for the development of small-scale dosimetry based on realistic biological data.

In the context of radioimmunotherapy (RIT), a single in vitro cultured carcinoma cell was considered in order to reproduce the effect of its realistic geometry and monoclonal antibody biodistribution on the dosimetry [START_REF] Arnaud | Complex cell geometry and sources distribution model for Monte Carlo single cell dosimetry with iodine 125 radioimmunotherapy[END_REF] . These biological data were extracted from microscope images. Assuming antibody radiolabeling with iodine-125, the energy deposition in the cell was calculated with Monte Carlo simulations. However, the relevance of this work was limited since tumors targeted in RIT such as non-Hodgkin lymphoma (NHL) grow as multicellular three-dimensional structures. Thus, a single cell dosimetric model neglects the absorbed dose delivered by the radiopharmaceutical bound to other cells of the tumor (crossfire effect). Furthermore, studies regarding unlabeled monoclonal antibody treatments highlighted that the biodistribution might be impacted by the tumor size and the extracellular matrix (i.e. the assemblage of non-cellular molecules providing structural support for the surrounding cells), limiting the effectiveness of the treatment. Increasing tumor size is regarded as possibly reducing the monoclonal antibody concentration and decreasing survival of a mouse model [START_REF] Daydé | Tumor burden influences exposure and response to rituximab: pharmacokinetic-pharmacodynamic modeling using a syngeneic bioluminescent murine model expressing human CD20[END_REF] and patients [START_REF] Meignan | Baseline metabolic tumor volume predicts outcome in high-tumor-burden follicular lymphoma: a pooled analysis of three multicenter studies[END_REF] . Furthermore, the extracellular matrix might hamper the monoclonal antibody penetration [START_REF] Netti | Role of extracellular matrix assembly in interstitial transport in solid tumors[END_REF] . To address these questions, a three-dimensional in vitro culture system of follicular lymphoma (a NHL subtype) expressing extracellular matrix was designed in a separate project at the Cancer Research Center of Toulouse 18- 21 . The model was termed multicellular aggregate of lymphoma cells (MALC). The MALCs were incubated with rituximab. This chimeric monoclonal antibody targeting the membrane-associated CD20 antigen has been combined with radionuclides for RIT treatments (see clinical review of Read et al. [START_REF] Read | The status of radioimmunotherapy in CD20+ non-Hodgkin's lymphoma[END_REF] ).

To improve the realism and the relevance of the single cell dosimetry, this paper presents a dosimetric model based on selective plane illumination microscopy (SPIM) images of three MALCs, of different size, treated by rituximab. Data were extracted to define the tumor geometry and the radiopharmaceutical biodistribution. Rituximab was virtually radiolabeled with two types of electron emitters.

The Auger electron emitters iodine-125 ( 125 I) and indium-111 ( 111 In) were studied because they proved their efficiency in destroying lymphoma cells in vitro in RIT-related research [START_REF] Griffiths | Cytotoxicity with Auger electron-emitting radionuclides delivered by antibodies[END_REF]24 . β -particle emitters lutetium-177 ( 177 Lu), iodine-131 ( 131 I) and yttrium-90 ( 90 Y) were also considered. For the RIT treatment of patients with relapsed / refractory NHL, 177 Lu-lilotomab satetraxetan is currently on trial [START_REF] Blakkisrud | Tumor-absorbed dose for non-Hodgkin lymphoma patients treated with the anti-CD37 antibody radionuclide conjugate 177 Lu-lilotomab satetraxetan[END_REF] , while both 131 I-tositumomab (Bexxar) and 90 Y-ibritumomab tiuxetan (Zevalin) have been used [START_REF] Goldsmith | Radioimmunotherapy of lymphoma: Bexxar and Zevalin[END_REF] . In RIT, the dosimetry is highly impacted by the radiopharmaceutical effective half-life, a parameter that would depend on the rituximab radiolabeling. We studied the energy deposition patterns in MALCs with Monte Carlo modeling of radiation transport. Additionally, the effective half-life was accounted for by using the simulated data and a simple pharmacokinetic model.

As a first step, the relevance of implementing realistic biological data was highlighted by comparing average absorbed doses to literature data. For the following step, to illustrate the possibilities offered by this new dosimetric model, we compared the radionuclides' absorbed fraction, average absorbed dose and absorbed dose volume distribution. As a last step, the influence of the radiopharmaceuticals half-life on the absorbed dose was studied.

Material and methods

Three-dimensional model of follicular lymphoma

The development of the MALCs as well as the microscopy acquisition were detailed in prior studies [START_REF] Gravelle | Innate predisposition to immune escape in follicular lymphoma cells[END_REF][START_REF] Gravelle | Cell growth in aggregates determines gene expression, proliferation, survival, chemoresistance, and sensitivity to immune effectors in follicular lymphoma[END_REF][START_REF] Decaup | Anti-tumor activity of obinutuzumab and rituximab in a follicular lymphoma 3D model[END_REF][START_REF] Rossi | Boosting γδ T cell-mediated antibody-dependent cellular cytotoxicity by PD-1 blockade in follicular lymphoma[END_REF] . This section focuses on information relevant to the present work.

MALC growth and labeling

To compare the antibody biodistribution in different tumor volumes, the MALCs, created with follicular lymphoma cells expressing green fluorescent protein (GFP), were cultured during 2, 5 or 10 days (hereafter referred to as D2, D5 and D10). Then, they were treated by a saturating concentration of 10 µg/ml of CF555 labeled rituximab [START_REF] Bezombes | Rituximab antiproliferative effect in B-lymphoma cells is associated with acid-sphingomyelinase activation in raft micodomains[END_REF] for a duration of 24 h.

SPIM imaging

Fluorescence microscopy technologies, such as SPIM detect the light emitted by fluorescent markers [START_REF] Huisken | Selective plane illumination microscopy techniques in developmental biology[END_REF] . SPIM was particularly adapted to this study since it allows the observation of whole organisms or three-dimensional multicellular models over a few millimeters while keeping a low phototoxicity [START_REF] Ichikawa | Live imaging and quantitative analysis of gastrulation in mouse embryos using light-sheet microscopy and 3D tracking tools[END_REF][START_REF] Weber | Multilayer mounting for long-term light sheet microscopy of zebrafish[END_REF] . 24 h after rituximab administration, z-stack of images were acquired. Each image corresponded to a slice of the MALC in the x-y plane. For example, Figure 1

Image processing

To extract the data required for the dosimetry of the MALCs, SPIM images were processed with Fiji [START_REF] Schindelin | Fiji: an open source platform for biological image analysis[END_REF] ("Fiji Is Just ImageJ", downloadable from fiji.sc). This program has specific features for the processing of microscope images of three-dimensional biological samples [START_REF] Eliceiri | Tools for visualizing multidimensional images from living specimens[END_REF] and can read TIFF files. First, the stacks of images were transferred to Fiji for 3D reconstructions (e.g. MALC D5 in Figure 2). Then, two macros were developed and implemented in Fiji:

-The first one treated the green signal in order to extract the MALCs geometry. MALC voxels were separated from the background by applying an intensity threshold. A clear delineation of the MALCs boundary allowed a straightforward determination of the threshold value. Since the black hollow sections of the MALC D5 and D10 were in the center, they were considered as a part of the volume of interest (VOI). Finally, in view of the definition of the MALC geometry during Monte Carlo simulations, the macro recorded different information into a ".dat" file; namely, the number of voxels, the voxel dimensions along the three axes, and each voxel coordinate (x, y, z).

-The second macro treated the red signal from the rituximab CF555 label. It removed the background noise, then recorded the voxel coordinates and intensity into a ".dat" file. During simulations, the source voxels were selected with a probability directly derived from their red signal intensity. It was supposed that the CF555 biodistribution would be identical to the radiopharmaceutical. This assumes that the CF555 labeling would neither detach from the antibody nor impact its biodistribution. To the best of our knowledge, no publication addresses these questions. This also assumes that the radiolabeling would not affect the rituximab biodistribution. This last point is supported by the requirement that, in nuclear medicine, the radiolabeling strategies preserve the biological properties of the proteins 33 such as rituximab.

2.3

Voxel-based dosimetry

General settings

Energy depositions were calculated with the Geant4 Monte Carlo open-source platform [START_REF] Agostinelli | Geant4-a simulation toolkit[END_REF][START_REF] Allison | Geant4 developments and applications[END_REF][START_REF] Allison | Recent developments in Geant4[END_REF] since it enables the implementation of voxelized geometries and its physics models are well adapted to energies encountered in TRT. Furthermore, Geant4, being a well-established object-oriented code base, has a high degree of reuse meaning that little code must be written for a new application. The Geant4 "microbeam" example served as a basis. It was developed at the Centre d'Études Nucléaires de Bordeaux-Gradignan (CENBG) to simulate the irradiation of individual cells with a microbeam of protons or α particles [START_REF] Incerti | Monte Carlo dosimetry for targeted irradiation of individual cells using a microbeam facility[END_REF][START_REF] Barberet | Monte-Carlo dosimetry on a realistic cell monolayer geometry exposed to alpha particles[END_REF] . In order to suit this example for TRT, we added the possibility of simulating the irradiation of a VOI (in this work the MALC) by internalized radiopharmaceuticals. To reconstruct the MALCs, the class G4PVParameterised placed the SPIM images voxels (3.7 • 3.7 • 3.1 µm 3 ) according to their coordinates stored in the output file of the first macro (see 2.2). The MALCs were positioned into a sphere of liquid water. This enabled electron backscattering into the MALCs with the consequence of an additional energy deposition. The simulation of the irradiation of each MALC by each radionuclide considered 10 6 decays. The activity was assumed homogeneously biodistributed within a voxel source.

Accordingly, the coordinates of a decay site inside such a voxel were sampled from a uniform distribution. The particles were then emitted isotropically. Depending on the type of electron emitter, energy depositions in the MALCs were simulated with different methods that will be presented in the following sections (2.3.2 and 2.3.3).

Auger electron emitters

125 I and 111 In mean emission spectra of monoenergetic electrons (Auger and internal conversion electrons, Figure 3) and photons (X-rays and γ-rays) were extracted from the publication of the American Association of Physicists in Medicine 39 (AAPM), as it is still a reference for Auger electron emitters (e.g. [START_REF] Arnaud | Complex cell geometry and sources distribution model for Monte Carlo single cell dosimetry with iodine 125 radioimmunotherapy[END_REF][START_REF] Taborda | Dosimetry at the sub-cellular scale of Auger-electron emitter 99m Tc in a mouse single thyroid follicle[END_REF][START_REF] Piroozfar | The effect of 111 In radionuclide distance and Auger electron energy on direct induction of DNA double-strand breaks: a Monte Carlo study using Geant4 toolkit[END_REF] ). Table 1 shows important properties of the decay of 125 I and 111 In. Half-lives were extracted from the Medical Internal Radiation Dose (MIRD) book [START_REF] Eckerman | MIRD: Radionuclide data and decay schemes[END_REF] . The mean penetration corresponds to the length of the vector connecting the electron emission and stopping points. Hence, this distance is lower than the path length because the latter considers the direction changes of the electrons. Although nuclear medicine dosimetric studies usually report the path length, in the present work the penetration was more relevant for analyzing the results. Table 1. Some physical properties of the Auger electron emitters. Half-lives come from MIRD data [START_REF] Eckerman | MIRD: Radionuclide data and decay schemes[END_REF] , spectra data were extracted from the AAPM publication [START_REF] Howell | Radiation spectra for Auger electrons emitting radionuclides: Report No. 2 of AAPM Nuclear Medicine Task Group No 6[END_REF] and mean penetrations were determined with Geant4-DNA. Auger electrons induce energy deposition clusters localized at the nanometer scale [START_REF] Kassis | The amazing world of Auger electrons[END_REF] . Such phenomena can be studied with step-by-step Monte Carlo codes since they give the track structure of the particles down to the nanometer scale. Accordingly, 125 I and 111 In emissions were simulated with Geant4-DNA [START_REF] Bernal | Track structure modeling in liquid water: A review of the Geant4-DNA very low energy extension of the Geant4 Monte Carlo simulation toolkit[END_REF][START_REF] Incerti | Comparison of GEANT4 very low energy cross section models with experimental data in water[END_REF][START_REF] Incerti | The Geant4-DNA project[END_REF][START_REF] Incerti | Geant4-DNA example applications for track structure simulations in liquid water: a report from the Geant4-DNA Project[END_REF] (geant4-dna.org). "CPA100" physics models ("option 6") [START_REF] Bordage | Implementation of new physics models for low energy electrons in liquid water in Geant4-DNA[END_REF][START_REF] Bordes | Low-energy electron transport with alternative physics models within Geant4-DNA code and radioimmunotherapy applications[END_REF] were used. They describe ionization, electronic excitation and elastic scattering processes and are applicable from 256 keV down to 11 eV in liquid water.

Auger electron emitter

For the simulations of photon interactions (Compton and Rayleigh scatterings as well as photoelectric absorption), Geant4-DNA does not provide any physics models to-date. Thus, Geant4 "Livermore" physics models for photons were selected [START_REF] Cullen | EPDL97: the evaluated photon data library `97 version[END_REF][START_REF] Ivanchenko | Combination of electromagnetic physics processes for microdosimetry in liquid water with the Geant4 Monte Carlo simulation toolkit[END_REF] . During simulations, electrons and photons could produce vacancies in the atomic shells so atomic deexcitation emission was tracked (Auger electron emission, Auger cascade and fluorescence).

β -particle emitters

AAPM's publication [START_REF] Howell | Radiation spectra for Auger electrons emitting radionuclides: Report No. 2 of AAPM Nuclear Medicine Task Group No 6[END_REF] only provides data for Auger electron emitters. Therefore, 177 Lu, 131 I and 90 Y mean emission spectra of β -particles (Figure 4), monoenergetic electrons and photons were extracted from the MIRD book [START_REF] Eckerman | MIRD: Radionuclide data and decay schemes[END_REF] . Table 2 presents the main properties of the decay of 177 Lu, 131 I and 90 Y. Table 2. Some physical properties of the β -particle emitters. Half-lives and spectra data come from MIRD book [START_REF] Eckerman | MIRD: Radionuclide data and decay schemes[END_REF] and mean penetrations were determined with Geant4. Part of the emission of 177 Lu, 131 I and 90 Y has an energy too high to be modeled with Geant4-DNA physics models. Thus, simulations of the irradiation with these emitters were accomplished with Geant4 in a condensed history mode, as this provides physics models applicable up to GeV energies. "Livermore" [START_REF] Perkins | Tables and graphs of electron-interaction cross sections from 10 eV to 100 GeV derived from the LLNL Evaluated Electron Data Library[END_REF][START_REF] Perkins | Tables and graphs of atomic subshell and relaxation data derived from the LLNL Evaluated Atomic Data Library (EADL)[END_REF] physics models were used to describe particles interactions.

The ionization process was simulated along with bremsstrahlung. This latter slightly increased the stopping power of the electrons / β -particles at these energies, thus marginally reducing their penetration. "Livermore" electron models are defined down to about 10 eV. Nevertheless, in accordance with guidelines, the cut-off was set to 250 eV. Energy production cuts for secondary particles was also fixed at 250 eV. In order to limit the number of electrons crossing a voxel with no energy deposition, step size limit was arbitrary chosen to 0.6 µm, corresponding to less than 1/5 of the voxel dimensions.

Dosimetric quantities

Simulation results

Geant4 simulations generated output files containing the energy deposited in the MALCs voxels.

Firstly, different dosimetric quantities were directly calculated from these data: the absorbed fraction and the absorbed dose.

-The absorbed fraction ϕ represents the fraction of energy released from the source region that is deposited in the target region. Here, the MALCs represent both the source and the target regions. The results focused on electron absorbed fractions ϕ e -(monoenergetic electrons and eventually β -particles), because photon absorbed fractions are insignificant at this scale:

ϕ e -MALC = E e - E e -0 (1)
Where E e -is the deposited energy and E e -0 the electron-released energy. The absorbed fraction is frequently reported in nuclear medicine dosimetry studies [START_REF] Loevinger | MIRD primer for absorbed dose calculations[END_REF] .

-The energy deposited E i (in J) in a VOI i (either a whole MALC or a voxel) characterized by a mass m i (in kg) was converted to absorbed dose D i (in Gy):

D i = E i m i (2)
This quantity was then divided by 10 6 to obtain the mean absorbed dose per decay (or per unit of time-integrated activity, in Bq • s, according to MIRD most recent terminology [START_REF] Bolch | MIRD Pamphlet No. 21: a generalized schema for radiopharmaceutical dosimetry-standardization of nomenclature[END_REF] ). First, this quantity was determined for the whole MALCs. However, the absorbed dose was not uniformly delivered within the volume. In order to specify the volume fractions having received either low or high absorbed dose, cumulative dose-volume histograms (cDVH) were plotted. They give the fraction of the volume (i.e. the fraction of MALC voxels) in which an absorbed dose greater than or equal to a given value was delivered.

These graphs are routinely used in external radiotherapy, but also in some TRT studies [START_REF] Ljungberg | The accuracy of absorbed dose estimates in tumours determined by quantitative SPECT: a Monte Carlo study[END_REF][START_REF] Marcatili | Model-based versus specific dosimetry in diagnostic context: comparison of three dosimetric approaches[END_REF] .

Considering radiopharmaceutical effective half-life

Subsequently, in order to consider the effective half-life of the radiopharmaceutical X-rituximab (X being one of the radionuclides), the MALCs and voxels absorbed doses per unit administered activity (A 0 in Bq) were calculated. In nuclear medicine, this quantity allows the assessment of the absorbed dose delivered to a VOI for a given administered activity. It was calculated as the product of the absorbed doses per decay (see 2.4.1) and the time-integrated activity coefficient a ̃ (also known as the residence time and given in s). The value of a ̃ represents the total number of decays occurring in a MALC per unit administered activity [START_REF] Bolch | MIRD Pamphlet No. 21: a generalized schema for radiopharmaceutical dosimetry-standardization of nomenclature[END_REF] :

a ̃= 1 A 0 ∫ A(X -rituximab, t) • dt ∞ 0 (3)
Where A(X -rituximab, t) is the time-dependent activity of the radiopharmaceutical in a MALC (in Bq). Assuming that this radiopharmaceutical has an instantaneous uptake in a MALC (i.e. instantaneously bound to the MALCs according to the biodistribution extracted from the images) and its elimination is monoexponential:

a ̃= 1 A 0 ∫ A 0 • e -(ln(2) T eff (X-rituximab) ⁄ )•t ∞ 0 • dt (4) a ̃= 1.443 • T eff (X -rituximab) (5)
Where T eff (X -rituximab), the effective half-life of the radiopharmaceutical (in s), depends on T phys (X)

the physical half-life of the radionuclide X (in s) and on T bio (rituximab) the biological half-life of rituximab (in s):

T eff (X -rituximab) = T phys (X) • T bio (rituximab) T phys (X) + T bio (rituximab) (6)
The physical half-lives are well documented (see [START_REF] Gustafsson | Uncertainty propagation for SPECT/CT-based renal dosimetry in 177 Lu peptide receptor radionuclide therapy[END_REF][START_REF] Gear | EANM pratical guidance on uncertainty analysis for molecular radiotherapy absorbed dose calculations[END_REF] . According to equation ( 6), the 3.4 days effective half-life corresponds to a rituximab biological half-life of 5.8 days. Using this value and the equation ( 6), the effective half-lives of rituximab radiolabeled with 125 I, 111 In, 177 Lu or 90 Y were deducted.

Results

MALC and rituximab biodistribution characteristics

Table 3 presents the dimensions of the three MALCs along x-, y-, z-axes and their masses in Geant4 simulations which were used to determine the absorbed doses according to equation (2). The masses were evaluated by multiplying the number of voxels in a MALC by the mass of a voxel filled with liquid water (as simulations were performed in this material). Figure 1 shows that MALCs considered in this study have ellipsoid-like shape. For MALC D2, D5 and D10, major radii are ~0.25 mm, ~0.5 mm and ~1.3 mm, respectively. Rituximab is heterogeneously biodistributed inside a shell-like volume at the periphery of the MALCs (Figure 1). Generally, the shell-like volumes are characterized by a thickness and a distance to the surface of the MALCs both of the order of a few hundredths of a millimeter. It was previously demonstrated that the extracellular matrix limits the diffusion of therapeutic agents in tumors [START_REF] Netti | Role of extracellular matrix assembly in interstitial transport in solid tumors[END_REF] . It was also shown that MALCs expressed extracellular matrix proteins [START_REF] Decaup | Anti-tumor activity of obinutuzumab and rituximab in a follicular lymphoma 3D model[END_REF] . Accordingly, in an ongoing MALC-based study, the possibility that the rituximab peripheral location is due to the extracellular matrix is being investigated. were compared to the values for spheres with radius, r, of 0.2 mm, 0.5 mm and 0.9 mm, respectively.

Furthermore, only 177 Lu, 131 I and 90 Y were considered because Bardiès and Chatal focused their study on β -particle emitters. The peripheral biodistribution of rituximab (discussed in section 3.1) can be seen as an intermediate between the surface and volume biodistributions assumed by Bardiès and Chatal.

Hence, this is coherent that the absorbed doses per decay for MALCs lie between the values for these two theoretical situations. It is also noticeable that results with these two biodistributions own a similar degree of agreement with results for the MALCs. This indicates that for the purpose of absorbed dose calculations, in these particular cases, the surface and the volume biodistribution are both equally relevant. Furthermore, non-negligible discrepancies between the theoretical model and realistic data were observed. These results highlight the need to account of realistic biological data.

Comparison of the radionuclides

For clarity, the results for each radionuclide are reported by ascending order of energy released per decay, that is: 125 I, 111 In, 177 Lu, 131 I and 90 Y.

Figure 6 displays, for each radionuclide, the electron absorbed fraction in the MALCs and Table 4 shows the corresponding numerical values. As MALCs get bigger, the electrons / β -particles deposit more energy inside before exiting them, increasing the absorbed fractions. Assuming a homogeneous surface biodistribution of rituximab (which is at first sight the closer approximation of the real biodistribution as seen on Figure 1), for very large MALCs the absorbed fraction would tend toward 0.5 (semiinfinite medium condition). However, for MALC D5 and D10, the values become greater than 0.5 for all radionuclides but 90 Y (presumably, absorbed fraction of 90 Y would become greater than 0.5 in bigger MALCs). This can be explained by the rituximab penetration inside the MALCs, enabling the electrons / β -particles to systematically deposit energy in the VOI, regardless their trajectories. This point is critical since it further demonstrates the importance of using realistic biological data. As radionuclides release more energy per decay, the penetration of the electrons / β -particles increase (see Table 1 andTable 2). Thus, more energy is deposited outside the MALCs, and in turn the absorbed fractions decrease. 125 I deposits nearly its whole electron energy within the three MALCs. The electrons it emits are unlikely to deposit energy outside the MALCs since their penetrations are limited to 1.5 •10 -2 mm while the rituximab is mostly located at more than ~3 •10 -2 mm from the edge. Note that the real values are inferior to 0.998 (rounded to 1.00 for consistency), due to scarce rituximab biodistribution closer to or at the surface. The absorbed fractions for 111 In, lying between 0.55 and 0.74, are notably lower than that of 125 I. The reason is that ~80% of the electron energy released per decay of 111 In is taken away by emissions of at least 145.0 keV, characterized by penetrations above 1.8 •10 -1 mm (right-hand side of Figure 3). Due to this property, the absorbed fractions for 111 In are remarkably close to 177 Lu, as their emissions have similar penetrations. Furthermore, the absorbed fractions for 131 I are higher than for 177 Lu as its penetration are higher (by 1.0 •10 -1 mm at 𝐸 β - 𝑚𝑒𝑎𝑛 ). Finally, the absorbed fractions for 90 Y, ranging between 0.04 and 0.12, are by far the lowest. Such low values are due to the penetrations of its emissions (up to 8.4 mm) being significantly superior to the size of the MALCs. Lower energy released by 177 Lu, is compensated by its higher absorbed fractions. Among β -particle emitters, 90 Y (933.2 keV/decay) emits by far the more energy, yet the absorbed doses delivered are lowest in MALCs D2 and D5 due to extremely low absorbed fractions. In MALC D10, however, it delivered the highest absorbed doses because, compared to MALC D2, absorbed fractions were increased by a factor 3. This benefit, illustrated by Figure 8 (d), (e) and (f), was brought by the higher penetration of their emission. It translates to higher cDVH values on the left-hand side (Figure 9). Among β -particle emitters, 90 Y energy deposition was the most uniform throughout the MALCs. In comparison, 177 Lu and 131 I deposited less energy in the central part but more in the periphery. The improved absorbed dose uniformity of 90 Y is further expressed by steeper gradients in its cDVH. It is linked to the variation of penetration and stopping power of β -particles with their energy. β -particles emitted by 90 Y having higher penetrations than for 177 Lu and 131 I, were more likely to reach the center of the MALCs. This difference is more pronounced as MALCs get bigger. Compared to 177 Lu and 131 I, emissions of 90 Y have lower stopping powers (according to NIST database [START_REF] Berger | ASTAR: computer programs for calculating stopping-power and range tables for electrons[END_REF] , at 𝐸 β - 𝑚𝑒𝑎𝑛 : 0.34, 0.29 and 0.19 keV/µm, respectively), thus they deposited less energy close to the decay site.

Absorbed dose per unit administered activity

The influence of the radiolabeled rituximab effective half-life on the absorbed dose is presented.

Table 6 summarizes the effective half-lives for radiolabeling with 125 I, 111 In, 177 Lu, 131 I or 90 Y as calculated from equation (6). The differences between the physical half-lives of the radionuclides were decreased by the introduction of the biological half-life of the rituximab (5.8 days). For each radionuclide, Figure 10 displays the mean absorbed dose per unit administered activity within the MALCs and Table 7 shows the corresponding numerical values. Regarding Auger electron emitters, the effective half-life of 125 I-rituximab was 2.8 times higher than 111 In-rituximab. Therefore, despite lower absorbed dose delivered per decay in MALC D5 and D10, 125 I delivered more absorbed dose per unit administered activity than 111 In. For β -particle emitters, the effective half-life of 177 Lurituximab and 131 I-rituximab were 1.7 and 1.9 times greater than that of 90 Y-rituximab. As a result, absorbed dose per unit administered activity for 177 Lu and 131 I were more than 3 times higher in MALC D2 and more than 2 times higher in MALC D5 than that of 90 Y. Furthermore, in MALC D10, despite higher absorbed dose per decay, 90 Y delivered less absorbed dose per unit administered activity than 177 Lu and 131 I. As discussed by Scheidhauer et al. [START_REF] Scheidhauer | Biodistribution and kinetics of 131 I-labelled anti-CD20 MAB IDEC-C2B8 (rituximab) in relapsed non-Hodgkin's lymphoma[END_REF] , the relatively long biological half-life of rituximab (5.8 days in NHL) is best matched by 131 I physical half-life (8.0 days) than by 90 Y short physical half-life (2.7 days).

177 Lu physical half-life of 6.6 days is also well adapted to the rituximab biological half-life. When the effective half-life of radiolabeled rituximab increases, the absorbed doses per administered activity increase. As a consequence, the region covered by high absorbed dose is largest for 131 I and smallest for 90 Y as illustrated by the isodoses (Figure 11). Furthermore, the higher the effective half-life, the greater the shift of the cDVH curves toward higher absorbed doses. Accordingly, in the region of high absorbed doses, cDVH values of 125 I become significantly higher than that of 111 In. In this region, cDVH values of 177 Lu and 131 I are noticeably higher than those of 90 Y.

Discussion

Our methodology allowed the observation of the absorbed dose in an in vitro tumor model which is more realistic than the spherical tumor with uniform biodistribution simplification. It is recognized, however, that the realism of this work is limited with respect to clinical context. It is therefore essential to review the degree of relevance of the biological and physical aspects of this work.

The present study focused on three MALCs with dimensions of the order of millimeter. On one hand, such a sample did not reproduce the diversity of tumor geometry and rituximab pharmacokinetics in patients, which are critical for dosimetry. For example, in clinical routine, the size of the diagnosed follicular lymphomas ranges from a few millimeters to several centimeters [START_REF] Wilder | Long-term results with radiotherapy for stage I-II follicular lymphomas[END_REF] . Additionally, multiple studies have highlighted the need for dosimetry in smaller follicular lymphomas (down to ~10 µm) [START_REF] Bousis | Monte Carlo single-cell dosimetry of Augerelectron emitting radionuclide[END_REF][START_REF] Bousis | Monte Carlo single-cell dosimetry of I-131, I-125 and I-123 for targeted radioimmunotherapy of B-cell lymphoma[END_REF][START_REF] Fourie | Microdosimetry of the Auger electron emitting 123 I radionnuclide using Geant4-DNA simulations[END_REF] .

Our study did not extend to the biggest and the smallest of the follicular lymphomas. It is expected that, with the increasing / decreasing size of the tumor model, radionuclides with higher (e.g. 90 Y) / lower (e.g. 125 I) penetration emissions, would become more favored from a dosimetric perspective. On the other hand, previous studies have highlighted the ability of three-dimensional aggregates to mimic the physiological in vivo aspects of tumors, enabling the identification of new targets and drugs [START_REF] Duval | Modeling physiological events in 2D vs. 3D cell culture[END_REF][START_REF] Thoma | 3D cell culture systems modeling tumor growth determinants in cancer target discovery[END_REF][START_REF] Hirschhaeuser | Multicellular tumor spheroids: An underestimated tool is catching up again[END_REF] . In particular, this tool has long been recognized as very effective to characterize the penetration of antibodies [START_REF] Sutherland | Penetration and binding of radiolabeled anti-carcinoembryonic antigen monoclonal antibodies and their antigen binding fragments in human colon multicellular tumor spheroids[END_REF][START_REF] Minchinton | Drug penetration in solid tumors[END_REF] . These previous works give confidence in the clinical relevance of the rituximab peripheral distribution in the MALCs, a key point of the present study.

The accurate Monte Carlo tracking of the radionuclide emissions in the MALCs relied largely on the choice of the physics models. Geant4-DNA provides two sets of physics models able to entirely simulate the electron emissions of 125 I and 111 In: the default ("option 2") and the CPA100 ("option 6") sets.

The latter set was selected following two main arguments. Differential cross sections for ionization and cross section for electronic excitation of "option 6" are in better agreement than that of "option 2" with experimental data. Furthermore, for low energy electrons, dose-point kernels generated by "option 2" are considered unrealistic because they are very diffusive [START_REF] Kyriakou | Technical Note: Improvements in GEANT4 energy-loss model and the effect on low-energy electron transport in liquid water[END_REF] , while those calculated by "option 6" are less so [START_REF] Bordes | Low-energy electron dose-point kernel simulations using new physics models implemented in Geant4-DNA[END_REF] . Concerning β -particle emitters, "Livermore" were considered more relevant than "PENELOPE-2008" [START_REF] Salvat | PENELOPE-2008: a code system for Monte Carlo simulation of electron and photon transport[END_REF] electron physics models based on two arguments. Firstly, to be in accordance with Geant4-DNA "option 6", the simulations were realized in liquid water. For this material, "Livermore" ionization cross sections are in better agreement with experimental data [START_REF] Allison | Recent developments in Geant4[END_REF] . Secondly, regarding deexcitation processes, "Livermore" simulates every atomic shell, while "PENELOPE-2008" only takes into account K, L and M shells [START_REF] Salvat | PENELOPE-2008: a code system for Monte Carlo simulation of electron and photon transport[END_REF] .

Monte Carlo simulations were informed with data from one SPIM acquisition, showing biodistribution in MALCs after 24 h of incubation. To approximate the time-dependent activity in the MALCs, we used the straightforward assumptions of a radiopharmaceutical instantaneous uptake and monoexponential elimination. In clinical studies, a more detailed description of the uptake and elimination is provided. Typically, 3 to 4 SPECT or PET acquisitions are performed, the first a few hours and the last between a few tenth of hours and nearly 200 h after administration [START_REF] Blakkisrud | Tumor-absorbed dose for non-Hodgkin lymphoma patients treated with the anti-CD37 antibody radionuclide conjugate 177 Lu-lilotomab satetraxetan[END_REF][START_REF] Ljungberg | The accuracy of absorbed dose estimates in tumours determined by quantitative SPECT: a Monte Carlo study[END_REF][START_REF] Avila-Rodriguez | Biodistribution and radiation dosimetry of [ 64 Cu]copper dichloride: first-in-human study in healthy volunteers[END_REF][START_REF] Santoro | Implementation of patient dosimetry in the clinical practice after targeted radiotherapy using [ 177 Lu-[DOTA0, Tyr3]-octreotate[END_REF] . Generally, these data are fitted by a mono/biexponential function. The instantaneous uptake does not account for the decays occurring outside the tumor during the uptake in a clinical context. Hence, the absorbed dose to the tumor delivered by radionuclides with short half-lives, such as 111 In and 90 Y, is overestimated compared to those with longer half-live. We note, however, that SPECT acquisitions of mice and patients showed early uptake of rituximab in NHL [START_REF] Camacho | Technetium-99m-or Cy7-labeled rituximab as an imaging agent for non-Hodgkin lymphoma[END_REF][START_REF] Stopar | 99m Tc-labelled rituximab, a new non-Hodgkin's lymphoma imaging agent: First clinical experience[END_REF] (a few hours after administration while 111 In and 90 Y physical half-lives are ~2 days). Regarding the elimination of the rituximab in NHL, no data were found in the literature.

Implementation of more sophisticated pharmacokinetic models in our methodology would be an important next step for a more relevant dosimetry. Ideally, these models should consider rituximab data in patients' NHL, upon availability.

Conclusion

We have developed a Monte Carlo dosimetric model based on realistic NHL geometry and rituximab biodistribution. It was employed to describe the absorbed dose that would be delivered at the microscopic scale by Auger electron emitters ( 125 I and 111 In) and β -particle emitters ( 177 Lu, 131 I and 90 Y) in the context of RIT.

A comparison between the mean absorbed dose per decay for realistic geometry and biodistribution with literature data for a sphere with either a surface or a volume biodistribution demonstrated significant differences. Then, for the five radionuclides, the absorbed fractions, the mean absorbed dose and its associated distribution were compared. This comparison further highlighted the importance of both the realistic MALC geometry and rituximab biodistribution as they often impacted the outcomes.

Aside from these two parameters, the decay properties -the electron / β -particle energy released per decay and penetration -had an important impact. The uniformity of the absorbed dose delivered throughout the MALCs increased for radionuclides releasing more energy per decay (in ascending order: 125 I, 111 In, 177 Lu, 131 I and 90 Y). The mean absorbed dose of β -particle emitters was higher than that of Auger electron emitters. 177 Lu and 131 I mean absorbed doses and were higher than 90 Y in the two smallest MALCs. Furthermore, when considering the radiolabeled rituximab effective half-life, the absorbed dose for 177 Lu and 131 I remained close to each other and considerably greater than that of 90 Y in the three MALCs.

The methodology developed may be applied to various situations of TRT. Voxelized images of others tumoral / healthy tissue models and tumor-seeking agents can be considered. It is possible to simulate the emission of other radionuclides of interest such as potential theranostic Auger electron emitters (e.g. copper-64 [START_REF] Avila-Rodriguez | Biodistribution and radiation dosimetry of [ 64 Cu]copper dichloride: first-in-human study in healthy volunteers[END_REF] , useful for PET imaging) and β -particle emitters (e.g. terbium-161 [START_REF] Müller | Therapeutic radiometals beyond 177 Lu and 90 Y: production and application of promising α-particle, β --particle, and Auger electron emitters[END_REF][START_REF] Müller | Scandium and terbium radionuclides for radiotheranostics: current state of development towards clinical application[END_REF] , useful for SPECT imaging) or therapeutic α particle emitters (e.g. actinium-225 [START_REF] Kratochwil | Targeted α-therapy of metastatic castrationresistant prostate cancer with 225 Ac-PSMA-617: dosimetry estimate and empiric dose finding[END_REF][START_REF] Robertson | Development of 225 Ac radiopharmaceuticals: TRIUMF perspectives and experiences[END_REF] ). To further improve the realism of the simulations, future work may integrate the temporal variation of the vector biodistribution and the growth of MALC. Such work could be useful in the case of real irradiation experiment to select relevant radionuclides and assess biological effects.

Figure 1 .

 1 displays a central slice of MALCs (a) D2, (b) D5 and (c) D10. Green signal corresponds to GFP 488 nm emission and reveals MALC geometry. The black hollow at the center of MALCs D5 and D10 are due to a lack of follicular lymphoma cells. Red signal from the 561 nm emission of CF555 represents the rituximab biodistribution. To optimize the contrast between follicular lymphoma cells and rituximab, the colors on the images differ from the real colors produced by the fluorescent markers. SPIM images were composed of voxels with (x, y and z) dimensions of 3.7 • 3.7 • 3.1 µm 3 and were stored into TIFF files. SPIM images of the central slices of the MALC: (a) D2, (b) D5 and (c) D10. Follicular lymphoma cell nuclei and rituximab appear in green and red, respectively.

Figure 2 .

 2 Figure 2. 3D reconstruction of MALC D5 from SPIM images. The axes are given in millimeters.

Figure 3 .

 3 Figure 3. Mean electron emission spectra (Auger and internal conversion electrons) of 125 I and 111 In from the AAPM publication 39 . For each energy decade, the top x-axis provides the corresponding mean penetration according to Geant4-DNA.

Figure 4 .

 4 Figure 4. β -particle emission spectra of 177 Lu, 131 I and 90 Y from MIRD book 42 . The top x-axis provides the corresponding Geant4 mean penetration every 500 keV.

Figure 5

 5 Figure 5 displays the relative difference (%) of absorbed dose per decay between the reference data of Bardiès and Chatal 3 for spheres with either (a) a surface or (b) a volume biodistribution and MALCs calculations (relative difference = 100 • (Dsphere -DMALC) / DMALC). In this previous study, calculations were carried out by an analytical method based on scaled dose-point kernels in spheres with radius ranging from 0.01 to 22 mm. For a valid comparison, the results for a MALC were compared to the results of the sphere having the closest mass. Accordingly, absorbed doses for MALCs D2, D5 and D10

Figure 5 .

 5 Figure 5. Relative differences (%) between absorbed dose per decay calculated in spheres of various radii with either a (a) surface or (b) a volume biodistribution 3 and MALCs (taken as reference). Comparison spheres vs. MALCs are as follow: r = 0.2 mm vs. MALC D2 (left), r = 0.5 mm vs. MALC D5 (middle) and r = 0.9 mm vs. MALC D10 (right).

Figure 6 .

 6 Figure 6. Absorbed fractions in MALCs D2, D5 and D10 for 125 I, 111 In, 177 Lu, 131 I and 90 Y.

Figure 7

 7 Figure 7 displays the mean absorbed dose per decay delivered to the MALCs by each radionuclide and Table 5 shows the corresponding numerical values. These results are impacted by the energy released by the radionuclides and the absorbed fractions.

Figure 7 .

 7 Figure 7. Mean absorbed dose per decay in MALC D2, D5 and D10 for 125 I, 111 In, 177 Lu, 131 I and 90 Y.

Figure 8

 8 Figure 8 illustrates, as an example the distributions of the absorbed dose per decay (in µGy•Bq - 1 •s -1 ) inside the voxels of the slice of MALC D5 shown in Figure 1 (b). Note that 0.01 µGy•Bq -1 •s -1 corresponds to a mean energy deposited per decay inside a voxel of 2.6 • 10 -3 eV•Bq -1 •s -1 . Figure 9 displays the cDVH (in µGy•Bq -1 •s -1 ) of MALCs (a) D2, (b) D5 and (c) D10 for 125 I, 111 In, 177 Lu, 131 I and 90 Y. Globally, when radionuclides release more energy, the absorbed dose is delivered more uniformly within the MALCs because of the increasing penetration of the emission.

Figure 8 (

 8 Figure 8 (b) highlights the peripheral irradiation of 125 I, with a pattern very similar to the rituximab biodistribution. However, this radionuclide delivered high absorbed dose as revealed by the lesssteep gradient of its cDVH. These observations are consistent with 125 I emitting low energy electrons, which deposit most of the energy in a small distance from the decay site. The highly localized absorbed dose delivered by 125 I was further observed at the cellular level by Arnaud et al.[START_REF] Arnaud | Complex cell geometry and sources distribution model for Monte Carlo single cell dosimetry with iodine 125 radioimmunotherapy[END_REF] . Unlike this radionuclide,111 In sparsely delivered absorbed dose in the central part of the MALCs (Figure8(c)) due to the higher penetrations of its emission. For both Auger electron emitters, the low absorbed dose delivered in the central part of the MALCs (or the lack thereof) is linked to the rituximab peripheral biodistribution. Nevertheless, the limited penetration of rituximab inside the MALCs can be counterbalanced by a radiolabeling with one of the β -particle emitters.

Figure 8 .Figure 9 .

 89 Figure 8. Absorbed dose per decay (µGy•Bq -1 •s -1 ) delivered to the voxel of the slice of MALC D5 shown in Figure 1

Figure 10 .

 10 Figure 10. Mean absorbed dose per unit administered activity to MALCs D2, D5 and D10 for 125 I, 111 In, 177 Lu, 131 I and 90 Y.

Figure 11

 11 Figure 11 illustrates, as an example the isodoses (in cGy•Bq -1 ) inside the voxels of the slice of MALC D5 shown by Figure 1 (b). Only the β -particle emitters were represented since the absorbed dose distribution of Auger electron emitters were very heterogeneous and did not allow for a clear representation of isodoses. Figure 12 shows the cDVH of MALC (a) D2, (b) D5 and (c) D10 per unit administered activity for 125 I, 111 In, 177 Lu, 131 I and 90 Y.

Figure 11 .Figure 12 .

 1112 Figure 11. Isodose distributions in the slice of MALC D5 shown in Figure 1 (b), per unit administered activity of (b) 177 Lu, (c) 131 I and (d) 90 Y (cGy•Bq -1 ). The corresponding color scale is displayed in (a).

Table 1

 1 We emphasize, however, that its precision is limited. It is associated to a standard deviation of 1.1 day. Also, Scheidhauer et al. reported a high inter-patient variability, with effective halflives ranging from 1.7 to 5.5 days. The gamma camera imaging introduces further uncertainties on the dosimetry

and Table

2

), but, rituximab biological half-life in NHL is less-well established. In a clinical study, Scheidhauer et al. evaluated the effective halflife of the radiopharmaceutical 131 I-rituximab in NHL

[START_REF] Scheidhauer | Biodistribution and kinetics of 131 I-labelled anti-CD20 MAB IDEC-C2B8 (rituximab) in relapsed non-Hodgkin's lymphoma[END_REF]

. Biodistribution in 14 patients was determined by gamma camera imaging. Our work focused on the 3.4 days mean effective half-life determined by

Scheidhauer et al.

Table 3 .

 3 Sizes along x-, y-and z-axes of the MALCs D2, D5 and D10.

		MALC	Dimensions x • y • z (mm 3 )	Mass (mg)
		D2	0.34 • 0.53 • 0.50	4.6 • 10 -2
		D5	0.96 • 1.09 • 1.07	5.5 • 10 -1
		D10	2.54 • 1.23 • 2.26	2.8
	3.2	Dosimetric quantities per decay
	3.2.

1 Comparison with literature data obtained in a simplified model

Table 4 .

 4 Absorbed fraction in the MALCs for the five radionuclides.

	Absorbed fraction

Table 5 .

 5 Mean absorbed dose per decay in the MALCs for the five radionuclides.Despite 111 In (34.7 keV/decay) releasing more energy than 125 I (19.5 keV/decay), the former delivered a lower absorbed dose in MALC D2 due to its lower absorbed fraction. Nevertheless, for MALC D5 and D10, the absorbed fraction for111 In increases, leading to greater absorbed doses than for 125 I.

			Mean absorbed dose per decay (µGy • Bq -1 • s -1 )
	Radionuclide	MALC D2	MALC D5	MALC D10
	125 I	6.78 • 10 -2	5.69 • 10 -3	1.10 • 10 -3
	111 In	6.26 • 10 -2	6.41 • 10 -3	1.37 • 10 -3
	177 Lu	2.27 • 10 -1	2.41 • 10 -2	5.34 • 10 -3
	131 I	2.01 • 10 -1	2.46 • 10 -2	5.84 • 10 -3
	90 Y	1.21 • 10 -1	1.98 • 10 -2	6.15 • 10 -3

In each MALC, absorbed doses for β -particle emitters are greater than for Auger electron emitters, because they release considerably more energy. For example, absorbed doses for 177 Lu are greater than that of

111

In by a factor ~4. Both have similar absorbed fractions, but 177 Lu (147.8 keV/decay) releases ~4 times more energy. For β -particle emitters, 177 Lu and 131 I (191.7 keV/decay) results are close.

Table 6 .

 6 Effective half-lives of rituximab radiolabeled with one of the five radionuclides in NHL, determined with equation (6).

	Radionuclide X	Effective half-life of X-rituximab (day)
	125 I	5.3
	111 In	1.9
	177 Lu	3.1
	131 I	3.4
	90 Y	1.8

Table 7 .

 7 Mean absorbed dose per unit administered activity to the MALCs for the five radionuclides.

	Mean

absorbed dose per unit administered activity (cGy • Bq -1 )

	Radionuclide	MALC D2	MALC D5	MALC D10
	125 I	4.50	3.78 • 10 -1	7.27 • 10 -2
	111 In	1.48	1.51 • 10 -1	3.23 • 10 -2
	177 Lu	8.76	9.33 • 10 -1	2.06 • 10 -1
	131 I	8.45	1.04	2.46 • 10 -1
	90 Y	2.78	4.56 • 10 -1	1.42 • 10 -1


F Richardeau
email: frederic.richardeau@laplace.univ-tlse.fr

F Boige

A Castellazzi
email: alberto.castellazzi@kuas.ac.jp

V Chazal

A Fayyaz

A Borghese

G Guibaud

SiC MOSFETs soft and hard failure modes: functional analysis and structural characterization

Keywords: Silicon carbide, SiC MOSFETs, short-circuit robustness, fail-to-open. I

This paper discusses first the discriminating phenomena yielding either a fail-to-short or fail-to-open circuit signature in 1200 V SiC power MOSFETs undergoing shortcircuit electro-thermal stress. Since fail-to-open behavior is of particular relevance to the application, the paper goes on to propose a benchmarking of a number of commercial devices, identifying a single product which offers consistent fail-to-open characteristics with bias voltages up to at least 50% of nominal rating. For that particular device, a thorough functional and structural characterization is presented. In particular, it is shown that: the gate current is an effective monitor of ensuing degradation under short-circuit stress and can be used to assess damage accumulation, as well as the reversible or permanent nature of device degradation; fail-to-open signature is associated with degradation of the gate-structure, with the creation of shortcircuits between the gate and source terminals in regions relatively far away from the active cells and not involving the field oxide. The findings are relevant to application of both discrete devices and multi-chip power modules, including multiple parallel connected dies.

INTRODUCTION

Silicon carbide (SiC) power MOSFETs are rapidly gaining increasing relevance in large market applications, including harsh operation and equally reliability sensitive environments, such as automotive or avionics. A key robustness requirement in such domains, particularly in the case of electric drive applications, is short-circuit (SC) withstand capability. A substantial body of literature has already been produced on the subject, pointing out potential criticalities in this respect, when transitioning from silicon (Si) transistors to SiC ones (see [START_REF] Castellazzi | Short-circuit tests on SiC power MOSFETs[END_REF][START_REF] Castellazzi | Short-circuit robustness of SiC Power MOSFETs: Experimental analysis[END_REF][START_REF] Castellazzi | SiC power MOSFETs performance, robustness and technology maturity[END_REF][START_REF] Romano | Short-circuit failure mechanism of SiC power MOSFETs[END_REF], for instance). In particular, the studies in [START_REF] Romano | A Comprehensive Study of Short-Circuit Ruggedness of Silicon Carbide Power MOSFETs[END_REF] voltage to a value well below the threshold voltage (as a result of damage in the gate structure, the gate-source resistance after damage is lower than the external drive gate resistance, RGS << RG,EXT), which practically results in the device turning off the device; the latter features a shorting of the drain-source terminals, with drain current, ID, thermal runaway and collapse of the drain-source voltage, VDS. The specific failure mode is of great relevance for the application, with a clear preference for FTO type transistors. Indeed, in contrast to FTO, the FTS mode typically implies: a) the fast discharge of energy storage components (e.g., capacitors), with oftentimes explosion-like characteristics and a possible failure propagation typically for a module device where chips are close to each other ; b) loss of functionality of the whole circuit in the case of parallel multichip architectures, either discrete or module based.

A preliminary screening on commercial 1.2 kV-rated MOSFETs delivered the results in Fig. 1 b): today a single manufacturer could be identified, whose device consistently performed as FTO for VDS values up to 50% rated voltage, which is a representative application requirement [START_REF] Boige | Ensure an original and safe "Fail-to-Open" mode in SiC MOSFET devices in extreme short-circuit operation[END_REF]. Because of its application interest, that device was considered further in detail. The gate leakage current, IG,LEAK, is an effective monitor of FTO behavior. Past work reported no damage accumulation and device aging for pulse widths below a minimum value, which could be related to a critical IG,LEAK value [START_REF] Chen | Study of short-circuit robustness of SiC MOSFETs, analysis of the failure modes and comparison with BJTs[END_REF]. So, in this work, first, tests were carried out to identify such critical pulse-width, tPW,CRIT. The results of Fig. 2 show: the SC drain current, ID, a); the gate-source voltage, VGS, b); IG,LEAK, c), for the device subjected to pulses of different duration. Here, VDS=600 V in all tests. During the pulse, at about time=10 s, IG,LEAK shows exponential increase onset, which, according to the results presented in [START_REF] Chen | Study of short-circuit robustness of SiC MOSFETs, analysis of the failure modes and comparison with BJTs[END_REF] and physically analyzed and modeled in [START_REF] Boige | Physical Origin of the Gate Current Surge During Short-Circuit Operation of SiC MOSFET[END_REF] as a Schottky emission current, was retained as the tPW,CRIT for the device under test (DUT) here.

The pulse width was then halved and repetitive stress was applied: Fig. 3 a) shows degradation of the saturation characteristics after thousands of pulses; however, that was not found to be related to any degradation at gate level, as the results of Fig. 3 b) clearly shows. On the other hand, increasing the pulse width yields very different results, Fig. 4: the drain current is still degraded, a), but that is now clearly due to damage of the gate-structure, manifest in the form of substantial IG,LEAK value increase, on an accelerated basis between #3000 and #3200 cycles , Fig. 4 b), and subsequent decrease, as a result of a low value series gate resistance in the transistor driving path, of the bias voltage applied across the gate-source terminals, Fig. 4 c).

III. STRUCTURAL CHARACTERIZATION

The results are interpreted on the basis of two different failure mechanisms intervening, depending on the pulse width: for "short" pulses (<tPW,CRIT) only reduced saturation current is manifest (Type 1 degradation), possibly caused by degradation of the chip metallization and subsequent VGS de-biasing [2; 9]; for "long" pulses, gate structure aging is evident (Type 2 degradation). The latter is caused by cracks around the interlayer dielectric (ILD) triggered by the high difference in thermal expansion between the top Al and the SiO2 and formation of metallic paths in the gate-source ILD, caused by the top Al melting at the edge of the chip during the shortcircuit, as presented in Figs. 5 a) and b) [START_REF] Boige | Ensure an original and safe "Fail-to-Open" mode in SiC MOSFET devices in extreme short-circuit operation[END_REF]. The region of major degradation and thus leading to eventual failure was observed to change on the device surface, but it was always located near the border of the active area, where terminations of various type are present and where geometrical singularity appears between the peripheral gate distribution ring and the interdigital ramifications fingers, Figs. 6 a) andb). However, the failure mechanism remained the same, even in the proximity of active cells, Fig. 6 c), which can be reasonably assumed to be the locations reaching highest temperatures.

IV. CONCLUSION

SiC power MOSFETs exhibit some key differences with Si MOSFETs and IGBTs. In particular, their features-set implies higher electro-thermal stress levels during short-circuit. On top of that, the present level of technology maturity is still associated with a greater spread in the value of some of the main device electro-thermal parameters, such as threshold voltage and on-state resistance. So, rather uneven degradation in parallel multi-chip structures is a likely occurrence. The LEAK,c). In this case, the change in device characteristics is strictly related to degradation of the gate, as manifest by the significant increase in gate current value. ability to ensure a soft (FTO) type of failure ensures the capability to still operate the module, even in the presence of a single die failure, enabling for instance, so-called "homehopping" system level operation mode in critical equipment such as automotive or avionics. Studies have already been undertaken to steer and control the occurrence of the desired FTO event in a single chip, without compromising the controllability of other parallel connected devices [START_REF] Boige | Shortcircuit robustness of parallel SiC MOSFETs andfail-safe mode strategy[END_REF][START_REF] Castellazzi | SiC MOSFET Device Parameter Spread and Ruggedness of Parallel Multichip Structures[END_REF]. A permanent gate leakage current of several mA is a salient precursor, easy to detect and monitor by the gate-driver circuit to reduce ageing and, for instance, trigger a supervisory system to stop operation at a suitable time to replace one or more devices in the system, before any failure occurs.

Fig. 1 :

 1 Fig. 1: short-circuit current ID and VDS, (top) and VGS (bottom), a); summary of SC failure mode (FTO or FTS) for different commercial devices at different VDS, b).

Fig. 2 :Fig. 3 :

 23 Fig. 2: Experimental results for SC tests with varying pulse durations: ID, a); VGS, b); IG,LEAK, c). A drop in VGS is noticeable (circle) from about 10 s, at which moment IG,LEAK features a steeper increase.

Fig. 4 :

 4 Fig. 4: Experimental results for repetitive SC test with "long" pulse widths: ID, a); VGS, b); IG,LEAK, c). In this case, the change in device characteristics is strictly related to degradation of the gate, as manifest by the significant increase in gate current value.

Fig. 5 :Fig. 6 :

 56 Fig. 5: Micrograph of MOSFET (bare-die top view) showing localized IG,LEAK increase, a); Scanning electron microscope (SEB) image of a device focused ion beam (FIB) microsection after destructive test, b).


F Richardeau

F Boige

Circuit-Type modelling of SiC power Mosfet in short-circuit operation including selective fail-to-open and fail-to-short modes competition

SiC Mosfet has very unique properties in extreme operation such as a well-known fail-to-short failure mode but in competition with a lesser known fail-to-open failure mode. These two modes are generally studied and modelled separately, whereas, in practice, they are coupled with the junction temperature of the chip. This paper presents a circuit-type modelling approach of these two modes simultaneously. This modelling allows to simulate the selectivity and competition between these two modes, one is clearly critical and the other is advantageously safe. The proposed model is then compared with short-circuit test of 1.2kV-80m MOSFET SiC.

Introduction and objectives

SiC Mosfet converters have been spectacularly developed to combine both high efficiency and a high level of compactness. Applications in the fields of transport, railway and also grid-connected PV inverters are targeted in the range of 10kW to several 100kW, in a medium range voltage, from 400V to 900V [START_REF] Lim | Design of 100 kVA SiC Power Converter for Aircraft Electric Starter Generator[END_REF]. The robustness of SiC Mosfet device and the management of its protection are also the subject of extensive research in many laboratories because its short-circuit time capability (TSC) is often lower than Silicon devices and varies greatly from one manufacturer to another [START_REF] Boige | Gate leakage-current analysis and modelling of planar and trench power SiC MOSFET devices in extreme short-circuit operation[END_REF]. The very thin gate-oxide (typically 50nm against 100nm for silicone-planar device) under a high gate field (4MV/cm compared to 10MV/cm breakdown field of the SiO2) is highly stressed in short-circuit operation. Indeed, the junction temperature can reach at least 1300K, leading to hot carrier injection in the oxide that could explain the possible weakness of the gate region and its variability in extreme operation [START_REF] Chiu | A Review on Conduction Mechanisms in Dielectric Films[END_REF], [START_REF] Boige | Physical origin of the gate current surge during short-circuit operation of SiC MOSFET[END_REF].

Thus SiC Mosfet has unique properties in short-circuit operation such as a competition between a wellknown drain-source fail-to-short (FTS) failure mode (also named "hard" failure mode, or again "catastrophic" failure mode) and a lesser known drain-source fail-to-open (FTO) failure mode (also named "soft" failure mode) caused by a gatesource short-circuit path. The FTS mode, comes from a thermal runaway and/or the parasitic BJT latch-up leading to the semiconductor local destruction and the creation of a metallic path drain-source [START_REF] Khatir | Experimental and numerical investigations on delayed short-circuit failure mode of single chip IGBT devices[END_REF]. The FTO mode, comes from the Al top metal fusion and diffusion, caused by the extreme heat in short-circuit, between the polysilicon gate and the Al source. The location of this metallic path can be either on top of the chip between the gate fingers and the source pad below the passivation layer, and/or across the ILD (Inter-Layer Dielectric) through metallized cracks [START_REF] Boige | Ensure an original and safe "fail-to-open" mode in planar and trench power SiC MOSFET devices in extreme short-circuit operation[END_REF]- [START_REF] Mbarek | Short-circuit robustness test and in depth microstructural analysis study of SiC MOSFET[END_REF] caused by a very low thermal expansion coefficient of the SiO2 compared to Al and SiC ones. In this last case, gate and source electrodes are shorted leading to a FTO mode.

FTS and FTO modes are generally studied and modelled separately, whereas in practice they are coupled with the temperature of the chip. However, the thermal runaway of the FTS mode is governed by the local temperature at the maximum field along the depletion region, while the FTO mode is dependent of the top metal melting and the ILD cracks which are governed by the oxide thermo-mechanical capability. This paper presents for the first time a circuit-type modelling approach of the SiC Mosfet including the global thermal dynamic of the FTS and FTO modes simultaneously. This modelling allows to simulate the selectivity and competition between these two modes, one is clearly critical and the other is advantageously safe. The proposed model is then compared with short-circuit test of 1.2kV-80m MOSFET SiC.

The paper is organized in four sections. Section 2 presents the circuit-type electro-thermal modelling of a SiC Mosfet based on VGSth, IDsat and IGleak thermaldependency. Section 3, deals with the FTO and FTS modelling separately but using the same circuit-type approach which will be embedded in the same model base. Finally, in section 4, the two failure modes are put in competition situation along two short-circuit test and compared with experimentation. The first one under a reduced drain-source voltage, where a stable FTO mode occurs. The second one, under a nominal drain-source voltage where the FTS mode happens after an instable and incomplete FTO mode.

Circuit-type electro-thermal modelling

The short-circuit operation shows a very fast thermal dynamics over a wide temperature range (from 300K to more than 1500K in a dozen of µs). Such a stress dynamic can reveal specific drift mechanisms or even degradation mechanisms of weak region which are not visible in nominal operation [START_REF] Boige | Gate leakage-current analysis and modelling of planar and trench power SiC MOSFET devices in extreme short-circuit operation[END_REF]. For example, the strong gate-source leakage-current is caused by a hot electron injection through the thin gate oxide [START_REF] Boige | Physical origin of the gate current surge during short-circuit operation of SiC MOSFET[END_REF] and the high drain-source peak current decreases as the channel mobility decreases at high temperature [START_REF] Mnatsakanov | Semiempirical model of carrier mobility in silicon carbide for analyzing its dependence on temperature and doping level[END_REF]. As a matter of fact, in order to model and fit the thermal-dependant physical mechanisms involved during the short circuit operation, an extended electrothermal modelling is needed. The authors previously presented in [START_REF] Boige | SiC power MOSFET in short-circuit operation: Electro-thermal macro-modelling combining physical and numerical approaches with circuit-type implementation[END_REF] a compact and non-segmented electro-thermal model including circuit elements as illustrated in Fig. 1. Note that in the following, all numerical values, unless otherwise specified, are fitted in relation with a 1.2kV-80m DUT, ROHM™ SCH2080 (generation II), based on an active area  12.4mm², in TO247 package.

First thermal modelling

Due to the very short transient operation (<100µs) clearly less than the thermal time constant of the chip, a diffusive 1D thermal model is sufficient. Then an 1D Comsol™ model of the chip's layers has been developed [START_REF] Boige | SiC power MOSFET in short-circuit operation: Electro-thermal macro-modelling combining physical and numerical approaches with circuit-type implementation[END_REF] including the electric field distribution Ex around the P+/N-depletion region. This model is fed by a time-dependant current density Jx(t) coming from experimental curves. At this stage, the gate regions, with polysilicon gate fingers, SiO2 and the ILD are not modelled because of their very small thickness (<1µm) and the possible mesh size problem in a FEM approach. A single temperature value is therefore considered in the thermal modelling, named Tj in the depletion zone at the maximum electric field position. However, the thickness of the top Al metal layer is only taken into account with an increase of the P-well thickness. Note that in SiC structure, the drift region is realized by a very fine epitaxy layer of about ten microns compared to a 300µm substrate. The gate oxide is also very thin (50nm) on which the Al layer is only 5µm thick. Given the small dimensions involved in the active zone, we have shown in [START_REF] Boige | SiC power MOSFET in short-circuit operation: Electro-thermal macro-modelling combining physical and numerical approaches with circuit-type implementation[END_REF] that the thermal gradient between the maximum electric field position and the gate region is very limited. Then authors only consider the temperature at the maximum electric field position of the depletion zone. In order to adapt this model into PLECS™, the 1D model can be either discretized at a constant geometric step x or fitted on a time-domain response to extract an equivalent Cauer's line (Rthi, Cthi) with a reduced number of cells (typically 30) or embedded by means of a vectorization formalism in a PLECS' functional block. Moreover, due to the "snap" effect of the electric field peak localized around the depletion region, a distributed power heat flux along the 1D x-axis is not necessary. Then a localized input instantaneous heat flux [W, t] connected at the input port of the Cauer's line is sufficient.

Electric modelling

It has been shown in [START_REF] Boige | Gate leakage-current analysis and modelling of planar and trench power SiC MOSFET devices in extreme short-circuit operation[END_REF] that the gate leakage-current can be modelled by an ohmic gate-source dipole whose value is in relation to the temperature inside the gate region. Relation [START_REF] Lim | Design of 100 kVA SiC Power Converter for Aircraft Electric Starter Generator[END_REF] gives the behavioral gate leakage-current model in relation to the adiabatic thermal energy accumulated during the short-circuit sequence (En) by means of a polynomial equivalent conductance. It is specifically fitted using the fit Matlab™ optimization function with experimental values. An offset function is also introduced which is independent of the gate polarization. In this relation, pij coefficients are in relation to the drain-source voltage applied to the DUT. For 600V, Tcase = 25°C and VGS  [18V, 21.5V], it was extracted : p01 = 9.51, Relations ( 2) is the well-known physical relations modelling the saturation current (IDsat). Relation (3) model the channel mobility (µn) in function of the temperature (Tj) using the Matheissen's mobility rule [START_REF] Potbhare | A Physical Model of High Temperature 4H-SiC MOSFETs[END_REF]. Relation (4) model the threshold voltage (VGS(th)) as an equivalent linear model [START_REF] Baliga | Modern power devices[END_REF].

I Dsat 𝑇 𝑗 = μ n 𝑇 𝑗 • C ox • Z 2L . 𝑉 𝑔𝑠 -V gs th 𝑇 𝑗 2 2 μ n 𝑇 𝑗 = μ a0 𝑇 𝑗 T 0 a 1 + μ a0 μ b0 𝑇 𝑗 T 0 b-a + μ a0 μ c0 𝑇 𝑗 T 0 c-a (3)
V GS th 𝑇 𝑗 = x 1 • 𝑇 𝑗 + x 2 4
Where: µa0, μb0, μc0 are mobility parameters at ambient temperature (300K) and a, b and c are sensitive thermal coefficients which have been fitted with experimental measures and temperature extracted from the thermal model presented in the previous subsection in order to obtain a wide range temperature model: [in cm².V -1 .s -1 ], µa0 = 2.33, μb0 = 2.95, µc0 = 1.94.

First validation of the electro-thermal model

The previous model is evaluated in a short-circuit simulation with a 12µs pulse for two values of VGS = 18V and 21V, results are presented in Fig. 2. In both cases, the model accurately represents the gate leakage-current with only a few % error. The drainsource current is slightly over-estimated for VGS = 21V while it is slightly under-estimated for VGS = 18V. However, the estimation of the adiabatic energy remains very accurate. Note this one greatly exceeds 1000K after 10µs. Such a temperature level is in agreement with recent publications using the same modelling approach but on different chips' size [START_REF] Pascal | Thermal simulations of SiC MOSFETs under shortcircuit conditions influence of various simulation parameters[END_REF] or considering a different VDS bias [START_REF] Tsibizov | Accurate Temperature Estimation of SiC Power MOSFETs under Extreme Operating Conditions[END_REF].

Fail-to-open and Fail-to-short modelling

In this section the previous electro-thermal model is extended to include FTO and FTS modes. These will be modelled separately because they represent different mechanisms which will be explained in this section.

FTO mode analysis and modelling

Previous thermal simulations shown that the maximum junction temperature near to the top surface can significantly exceed the melting temperature of the Al layer (660°C or 933K) [START_REF] Boige | SiC power MOSFET in short-circuit operation: Electro-thermal macro-modelling combining physical and numerical approaches with circuit-type implementation[END_REF]. Some failure Fig. 3. Experimental FTO and FTS modes in type 1 short-circuit operation for VDS = 600V (cascaded modes) in blue curve and for VDS = 400V in black curve, (@Tjnit = 25°C, Vbuffer(on/off) = 18V/-5V and RG=47Ω). [18], [START_REF] Romano | A Comprehensive Study of Short-Circuit Ruggedness of Silicon Carbide Power MOSFETs[END_REF].

V GS (V) I G (mA) I DS (A) V DS (V) I DS V DS V DS =600V V DS =400V
As shown in Fig. 3, below a given drain-source voltage, the FTO mode is guaranteed. Above this voltage the FTO mode is overtaken by the FTS mode which will be described in §3.2. Indeed, for a reduced drain-source voltage, the heat power density is limited, which therefore reduces the junction thermal dynamics of the chip, allowing a sufficient delay time for the top Al metal to diffuse under the passivation of the chip to form a stable Al metal bridge leading to the device natural turn-off and cooling [START_REF] Romano | A Comprehensive Study of Short-Circuit Ruggedness of Silicon Carbide Power MOSFETs[END_REF]. In addition, thermal shocks at the Al-SiO2/ILD and Al-SiC interfaces yield significant differential strains which causes cracks and metal infiltration between source electrode and gate fingers as it was seen on others SiC devices [START_REF] Boige | Investigation on damaged planar-oxide of 1200V SiC power MOSFETs in non-destructive short-circuit operation[END_REF]- [START_REF] Wang | Short-Circuit Ruggedness and Failure Mechanisms of Si/SiC Hybrid Switch[END_REF]. These two metallurgical and thermo-mechanical mechanisms allow to understand the origin of the FTO mode [START_REF] Boige | Ensure an original and safe "fail-to-open" mode in planar and trench power SiC MOSFET devices in extreme short-circuit operation[END_REF], [START_REF] Othman | Etude d'interrupteurs en carbure de silicium et potentiel d'utilisation dans des applications aéronautiques[END_REF], [START_REF] Romano | A Comprehensive Study of Short-Circuit Ruggedness of Silicon Carbide Power MOSFETs[END_REF].

In Fig. 4, the top metal temperature estimation from the experimentation is presented the colours of Fig. 3 are respected. The temperature is estimated with an enhanced version of the 1D thermal model presented in §2.2. The thermal sensitivity of the thermal conductivity and the specific heat capacity of the chip materials (Al and SiC) are also modelled from [START_REF] Hatch | Aluminum: Properties and Physical Metallurgy[END_REF], [START_REF] Snead | Handbook of SiC properties for fuel performance modeling[END_REF] as well as the Al phase-transition from solid to liquid state around 660°C (933K). Based on [START_REF] Hatch | Aluminum: Properties and Physical Metallurgy[END_REF], the heat capacitance in the solid state and the liquid state is displayed in Fig. 5. Relation [START_REF] Khatir | Experimental and numerical investigations on delayed short-circuit failure mode of single chip IGBT devices[END_REF] gives respectively the relation between the solid state capacitance to the fusion state model to the liquid state model.

C p = C pS T • 1 -α T + C pL • α T + L f • dα dT (5)
With : CpS, solid state heat capacitance which is temperature-dependent [START_REF] Hatch | Aluminum: Properties and Physical Metallurgy[END_REF]; Lf, heat of fusion (Lf = 397 kJ/kg [START_REF] Hatch | Aluminum: Properties and Physical Metallurgy[END_REF]); CpL, liquid state heat capacitance which is typically considered as a constant value (Cpl = 1.175kJ/(kg.K) [START_REF] Hatch | Aluminum: Properties and Physical Metallurgy[END_REF]); 𝛼 𝑇 is a phase transition function chosen as a cumulative function of a normal distribution (6) as also illustrated in Fig. 5. The normal distribution is centered in Tf and with a variance 𝜎 = 5𝐾. The solid and liquid phases of the heat capacitances are combined with the phase transition function to have a smooth transition during fusion. The integration of the normal law gives 1 and multiplying it by the heat of fusion gives the amount of latent heat released during the phase transition. The transition will take place on the interval 6𝜎 = Δ𝑇 𝑠→𝑙 .

For aluminum and also to ensure the convergence of the numerical model, Δ𝑇 𝑠→𝑙 = 30𝐾. During that time, the material is supposed in a mixed state between solid and liquid.

𝛼 𝑇 = ∫ 𝑑𝛼 𝑑𝑇 𝑥 𝑑𝑥 𝑇 -∞ (6)
The proposed modelling is by default implemented in Comsol™. The solid line curves in Fig. 4 show that it allows to simulate by a quasi-plateau corresponding to the quantity of latent energy to be dissipated in order to produce the Al layer phase change. [START_REF] Pascal | Thermal simulations of SiC MOSFETs under shortcircuit conditions influence of various simulation parameters[END_REF] recently proposed the same modelling approach based on our work by including a sensitivity study to the physical parameters. Concerning the circuit-based electro-thermal modelling of the DUT, these curves are fitted by a Cauer's line of 8 elements with constant parameters to easily allow an integration in PLECS™. The responses of the Cauer's lines are superimposed in a dotted line in Fig. 4. This line represents the simulated

FTS area

FTO permanent

Competition area FTO area response in a very average way and therefore above all in a qualitative way.

The conditions to access the FTO area safely is presented by means of a "green safe area" in the transient thermal domain: the junction temperature must be at least 1500K and at least 17µs after the total melting of the Al, called delay time (TD). The reference [START_REF] Boige | Ensure an original and safe "fail-to-open" mode in planar and trench power SiC MOSFET devices in extreme short-circuit operation[END_REF] showed that the quasi same FTO temperature is obtained for a different chip manufacturer but with a smaller size and therefore with a lower Tsc. Only 100K separates the two FTO temperatures of the two chips while the sizes are in a 1.2 ratio and the Tsc in a 2 ratio. These results show that the FTO thermal level is not very variable over a wide range of chip characteristics.

Given the significant temperature difference involved, we believe that the two failure mechanisms are well differentiated. Moreover, experimentally, the authors have verified the existence of these two types of failure-modes on a set of about ten DUTs components. The dispersion is not significant either on the TSC value (17µs) or on the specific drainsource voltage value, VDS  400V, that define the boundary between the two modes.

FTS mode analysis and modelling

As seen in Fig. 3 with the blue curve (VDS = 600V), the FTS mode is triggered before the consolidation of a complete FTO mode. Such a model concerns the representation of the drain-source leakage current by thermal generation [START_REF] Khatir | Experimental and numerical investigations on delayed short-circuit failure mode of single chip IGBT devices[END_REF]. A Shockley-Read-Hall (SRH) carriers generation model is used by means of relation [START_REF] Ouaida | Gate Oxide Degradation of SiC MOSFET in Switching Conditions[END_REF] [START_REF] Khatir | Experimental and numerical investigations on delayed short-circuit failure mode of single chip IGBT devices[END_REF], [START_REF] Wang | Temperature-Dependent Short-Circuit Capability of Silicon Carbide Power MOSFETs[END_REF] because it allowed to represent globally the drain leakage obtained in an experimental way. The thermal runaway dynamic of this leakage leads to the FTS mode. In SRH generation , WD is the space of charges width in short-mode under VDS defined in relation [START_REF] Reigosa | Effect of short-circuit stress on the degradation of the SiO2 dielectric in SiC power MOSFETs[END_REF] and ni the intrinsic density of charge carriers described in relation [START_REF] Wang | Short-Circuit Ruggedness and Failure Mechanisms of Si/SiC Hybrid Switch[END_REF] [START_REF] Baliga | Silicon Carbide Power Devices[END_REF].

𝐼
With: g, is the SRH generation lifetime [1ns, 1s] from [21][24]. We used this parameter g as a fitting variable (fixed at 1ns) so that our model produces the same drain -source current and failure time at 600V than the experiment. This fit was also evaluated at 400V giving a Tsc value in FTO mode very close to the experiment. Our model is therefore robust. The results of this model is presented in the next section in the simulation vs experimentation.

FTO and FTS modes competition: simulation vs experimentation

The complete proposed model is depicted in Fig. 6. In addition to electro-thermal modelling of the gate leakage-current, drain-source current and threshold voltage as described in §2, a dedicated FTO and FTS circuit-parts have been included. FTO mode is modelled by a gate-source parallel failed-resistor Rgdef (0.1) controlled by switch S3 and associated with switch S2 to force the channel turn-off. These two switches are managed by the functional block diagram presented in Fig. 7: if the top metal estimated temperature is higher than the Al fusion temperature (933K) and the FTO threshold temperature TFTO after the delay time Td (17µs), the FTO mode is triggered by a latch gate. As illustrated in Fig. 4, the estimated temperature comes from a simplified fitted Cauer's line on the phase-transition curve of the top metal temperature. This thermal line is named Ttop in Fig. 6. FTS mode is modelled by a parallel leakage drainsource controlled current-source Ith and FIth block that includes the SRH generation model studied in §3.2. This block is fed by the same thermal line Ttop. If the estimated temperature is higher than TFTS, the FTS mode is triggered and a dedicated switch S1 force the turn-on through a failed resistor Rdef (5m).

The evaluation of the model is given in Fig. 8 in comparison with experimental results from a testbench described in [START_REF] Boige | Gate leakage-current analysis and modelling of planar and trench power SiC MOSFET devices in extreme short-circuit operation[END_REF]. At nominal voltage, in Fig. 8a), due to the high thermal dynamic stress applied to the chip, the proposed model is able to highlight the cascaded FTO and FTS modes: the drain-source leakage current is too high to prevent thermal runaway. Note that the dedicated drain-source leakage component is drawn in dashed line. At reduced voltage but with the same pulse duration, in Fig 8b ), due to the low thermal dynamic stress applied to the chip, the FTO mode is the unique and stable mode. The low drain-source leakage current is clearly not enough to trigger the thermal runaway and the FTS mode does not appear.

Conclusion

In this article, the authors proposed a complete electrothermal circuit-type model for SiC Mosfet in short-circuit operation able to reproduce the new FTO mode in conjunction with the well-known FTS failure modes. The modelling of the dynamic drain-source carriers thermal generation combined with the FTO and FTS thermal thresholds allow to predictively simulate and to highlight either a specific cascaded failure modes at nominal voltage or an original safe FTO mode at low voltage as analysed from experimentation. This approach complements and extends in an original way the known electrothermal power device models. Finally, the proposed compact model has the following uses and interest:

-modelling and prediction of the TSC and ESC (robustness study model), -modelling and prediction of failure-modes (failure modes study model), -modelling separately IGS (leakage), VGS drop and IDSS (leakage) (diagnostic study model).

-fully compatible with gate driver protection design or fault propagation analysis in different conversion architectures.

Fig. 1 .

 1 Fig. 1. Electro-thermal PLECS™ Circuit-type modelling-base for SiC Mosfet in short-circuit operation. FIGleak, FVGSth, FIDsat are respectively modelling function for IGleak(Tj), VGSth(Tj) and IDsat(Tj).

  10 3 and: a = -1.48, b = -4.8, c = 2.56. In the same conditions thee threshold voltage parameters, x1 and x2, have been extracted: x1 = -0.025 V/K, x2 = 15.4V.

Fig. 2 .

 2 Fig. 2. Type 1 short-circuit SiC MOSFET model evaluation in PLECS™ software for different buffer voltages bias (@VDS = 600V, Tjinit = 25°C and RG = 47Ω).

Fig. 5 .Fig. 4 .

 54 Fig. 5. Thermal-dynamic heat capacitance modelling of Al layer for FTO modelling. Temperature (K) T f

Fig 6 .Fig 7 .

 67 Fig 6. Final proposed electro-thermal SiC Mosfet model for short-circuit simulation including FTO and FTS modes modelling.

Fig 8 .

 8 Fig 8. Comparison between the proposed model and the short-circuit failure modes experimentation at Vbuffer = 18V : a) cascaded FTO and FTS modes for VDS = 600V, b) unique and stable FTO mode for VDS = 400V.

  𝑡ℎ 𝑇 =

		𝑞𝑊 𝐷 𝑛 𝑖 𝑇 𝜏 𝑔	S	(8)
	𝑊 𝐷 = √ 2𝜖 𝑠 𝑞	•	𝑁 𝑑 +𝑁 𝑎 𝑁 𝑑 𝑁 𝑎	• 𝑉 𝐷𝑆	(9)
	𝑛 𝑖𝑆𝑖𝐶 𝑇 = 1.7 • 10 16 • 𝑇	3 2 • 𝑒 -2.08•10 4 𝑇

Acknowledgements

This research work received financial support from the French National Research Agency (ANR). Project name: HIT-TEMS managed by CNAM Paris and Satie lab.