Albert Saporta

Ane Etxebeste

Théo Kaprelian

Jean Michel Létang

David Sarrut

Modeling families of particle distributions with conditional GAN for Monte Carlo SPECT simulations

Keywords: Conditional GAN, Wasserstein, Monte Carlo, SPECT simulation

Objective. We propose a method to model families of distributions of particles exiting a phantom with a conditional Generative Adversarial Network (condGAN) during Monte Carlo simulation of SPECT imaging devices. Approach. The proposed condGAN is trained on a low statistics dataset containing the energy, the time, the position and the direction of exiting particles. In addition, it also contains a vector of conditions composed of four dimensions: the initial energy and the position of emitted particles within the phantom (a total of 12 dimensions). The information related to the gammas absorbed within the phantom is also added in the dataset. At the end of the training process, one component of the condGAN, the generator (G), is obtained. Main results. Particles with specific energies and positions of emission within the phantom can then be generated with G to replace the tracking of particle within the phantom, allowing reduced computation time compared to conventional Monte Carlo simulation. Significance. The condGAN generator is trained only once for a given phantom but can generate particles from various activity source distributions.

Introduction

Monte Carlo simulation in medical physics is widely used as it is the reference simulation method of particle tracking in matter. It allows accurate modeling of the physical interactions between particles and matter, such as Compton or Rayleigh scattering, photo-electric effect or Bremsstrahlung radiation. This approach is thus heavily used in the design and development of nuclear imaging systems such as Positron Emission Tomography (PET) or Single Photon Emission Computed Tomography (SPECT). For example, the design of new SPECT imaging devices [START_REF] Brown | In-silico optimisation of tileable philips digital SiPM based thin monolithic scintillator detectors for SPECT applications[END_REF][START_REF] Massari | Super Spatial Resolution (SSR) method for small animal SPECT imaging: A Monte Carlo study[END_REF][START_REF] Auer | Preliminary investigation of a Monte Carlo-based system matrix approach for quantitative clinical brain 123 I SPECT imaging[END_REF] or the development of reconstruction algorithms require realistic Monte Carlo simulations in various configurations. Such simulations create a mapping from a given activity source distribution inside a patient or a phantom to a signal captured by the imaging device outside of the patient/phantom by tracking particles one-by-one through the objects present in the simulation. The Monte Carlo simulation of a SPECT imaging device can be decomposed into two steps. During the first step, gammas emitted from an activity distribution of a given radionuclide are tracked in a medium (a phantom or a patient), potentially undergoing Compton scattering, until they are absorbed or exit the medium. During the second step, the interactions of the gammas within the detector are simulated. Decomposing a simulation is useful to avoid redundancy in certain applications, as the computation time to perform such simulations is generally high. For example, when studying different imaging systems, only the second simulation step has to be repeated while the first step is unchanged. For the first step, particles exiting the medium can be stored in a phase-space file and reused later. Depending on the complexity of the simulated configuration, billions of particles should be tracked to reach an acceptable statistical convergence, making the whole process usually very long. Indeed, the time needed for the simulation of the transport and tracking of gammas in a phantom or a patient can be long, due to the various possible interactions. Moreover, phase-space files are generally large (up to several GB) and can be cumbersome to process, use and exchange. Several works [START_REF] Brualla | PENELOPE/PRIMO-calculated photon and electron spectra from clinical accelerators[END_REF][START_REF] Chabert | Development and implementation in the Monte Carlo code PENELOPE of a new virtual source model for radiotherapy photon beams and portal image calculation[END_REF][START_REF] Grevillot | Simulation of a 6 MV Elekta Precise Linac photon beam using GATE/GEANT4[END_REF] provided methods to model radiotherapy Linac phase-space distributions analytically, but they have never been investigated for SPECT simulations. For SPECT simulations, several variance reduction techniques have been developed to improve the simulation efficiency, such as the angular response function (ARF) [START_REF] Song | Fast modelling of the collimator-detector response in Monte Carlo simulation of SPECT imaging using the angular response function[END_REF][START_REF] Descourt | Implementation of angular response function modeling in SPECT simulations with GATE[END_REF][START_REF] Sarrut | Learning SPECT detector angular response function with neural network for accelerating Monte-Carlo simulations[END_REF], forced detection (FD) [START_REF] Beenhouwer | Acceleration of GATE SPECT simulations[END_REF][START_REF] Beck | Analysis of spect including scatter and attenuation using sophisticated monte carlo modeling methods[END_REF][START_REF] Ljungberg | Attenuation and scatter correction in SPECT for sources in a nonhomogeneous object: a Monte Carlo study[END_REF], fixed forced detection (FFD) [START_REF] Cajgfinger | Fixed forced detection for fast SPECT Monte-Carlo simulation[END_REF] or Geometrical Importance Sampling (GIS) [START_REF] Haynor | The use of importance sampling techniques to improve the efficiency of photon tracking in emission tomography simulations: Use of importance sampling techniques to improve photon tracking[END_REF]. GPU-based approaches were also developed and implemented in softwares dedicated to speeding up the simulation of SPECT devices [15; 16]. In previous works [17; 18], a method was proposed to model complex particles phase-space by Generative Adversarial Network (GAN). The concept of GAN was proposed as a deep neural network architecture [START_REF] Goodfellow | Generative Adversarial Networks[END_REF] allowing to model multidimensional distributions. In particular, in [START_REF] Sarrut | Modeling Complex Particles Phase Space with Gan for Monte Carlo Spect Simulations: A Proof of Concept[END_REF], the authors modeled the distribution of particles exiting a patient or a phantom during Monte Carlo simulation of SPECT imaging devices. Once trained, one of the components of the GAN, a neural network called generator (G), serves as a compact and fast source of particles for the Monte Carlo simulation. This approach has several advantages. It allows reduced computation time compared to conventional Monte Carlo simulation, as there is no need to track the particles within the phantom or the patient. The file containing the GAN parameters is small (few MB). However, this approach requires to train a new GAN each time a parameter is modified, such as the energy of the emitted particles or the initial activity distribution inside the phantom or patient. Note that the concept of exploiting GAN within Monte Carlo simulations is also currently explored in the high energy physics community [START_REF] Carminati | Calorimetry with Deep Learning: Particle Classification[END_REF][START_REF] Erdmann | Generating and refining particle detector simulations using the Wasserstein distance in adversarial networks[END_REF][START_REF] Paganini | CaloGAN: Simulating 3D High Energy Particle Showers in Multi-Layer Electromagnetic Calorimeters with Generative Adversarial Networks[END_REF][START_REF] Paganini | Accelerating Science with Generative Adversarial Networks: An Application to 3D Particle Showers in Multi-Layer Calorimeters[END_REF]. For example, label conditioning using Auxiliary Classifier GAN (AC-GAN) was introduced in [START_REF] Odena | Conditional Image Synthesis With Auxiliary Classifier GANs[END_REF] and applied to Monte Carlo simulations of electromagnetic showers in [START_REF] Erdmann | Precise simulation of electromagnetic calorimeter showers using a Wasserstein Generative Adversarial Network[END_REF]. In this paper, we first extend the architecture proposed in [START_REF] Sarrut | Modeling Complex Particles Phase Space with Gan for Monte Carlo Spect Simulations: A Proof of Concept[END_REF] such that the GAN has to be trained only once on a dataset containing all the desired activity sources, for a given phantom. This is achieved thanks to conditional GAN (condGAN) [19; 26] that allows to model families of particle phase-space distributions. We also propose and evaluate a new parametrization of the training dataset. The structure of the paper is as follows. In section 2, we describe the Monte Carlo simulation setup, the training dataset, the proposed parametrization, the condGAN architecture, and the evaluation method. The experiments and their results are presented in sections 3 and 4. Finally, sections 5 and 6 are left for discussion and conclusion.

Materials and methods

Overview

We focus on the first step of a Monte Carlo simulation of a SPECT imaging device: the emission and transport of gammas in a phantom. In this step, the gammas emitted from a source distribution are tracked until they are absorbed or exit the phantom. The gammas exiting the phantom can be stored in a phase-space located on a spherical surface surrounding the phantom. Following the approach described in [17; 18], the goal is to create a "forward model", able to produce exiting gammas distributions from a user-defined source distribution, for a given phantom. We split the work into 2 steps. In a first step, a training dataset is obtained via a low statistics Monte Carlo simulation. In the second step, a condGAN is designed to learn a family of distributions of gammas exiting a given phantom, whatever the input source.

Monte Carlo generated training datasets

Without loss of generality, we considered a Monte Carlo simulation of the NEMA IEC phantom (see left image in figure 1 and detailed description in section 3) and of a 3D CT image of a patient (right image in figure 1). The first training dataset was obtained considering gammas that were emitted isotropically from a gamma source uniformly distributed within the NEMA IEC phantom. We considered gammas emitted with 8 energies in the range 113-637 keV, that corresponds to the principal gamma emission lines of 4 radionuclides commonly used in SPECT imaging ( 99m Tc, 111 In, 131 I, 177 Lu) as shown in table 1. The second training dataset was obtained considering gammas that were emitted isotropically from a gamma source uniformly distributed within the patient contour in the 3D CT image (right image of figure 1). We considered gammas emitted with an energy of 140.5 keV ( 99m Tc).

Table 1: Initial energies (E 0 ) considered for the Monte Carlo simulations of the sources emitting gammas within the phantom. These energies correspond to the principal gamma emission lines of 4 radionuclides [START_REF]Decay radiation database[END_REF] commonly used in SPECT imaging.

Radionuclide

99m Tc 111 In 131 I 177 Lu Initial energy E 0 (keV) 140.5 171. 3, 245.4 284.3, 364.5, 637.0 112.9, 208.4 The simulations for the two training datasets were performed with a phase-space scorer extension in GATE that records every gamma that reach the spherical surface surrounding the phantom that was added to the simulation. The emitted gammas were tracked in the phantom and potentially underwent Compton scattering until they were absorbed or exited the phantom. The gammas exiting the phantom and reaching the surface were stored in the phase-space training dataset described with 12 dimensions. The first 8 dimensions correspond to the parameters of the exiting gamma: kinetic energy E kin , 3D position P, 3D direction D and time t (time elapsed between emission and arrival to the phase-space surface). The information related to the gammas absorbed within the phantom that do not reach the surrounding spherical phase-space volume was artificially added in the phase-space. Their kinetic energy was set to 0, and 3D positions and time were obtained by projecting their trajectories on the spherical phasespace volume considering their initial 3D positions and directions. Note that the time value is not used in conventional SPECT systems because this value cannot be known. However, it has been included here to facilitate the evaluation (see subsection 2.6). The additional 4 dimensions are the initial conditions of the emitted gammas: the energy at emission E 0 and the 3D position of emission within the phantom P 0 .

Parametrization of the training dataset

We considered two parametrizations of the phase-space training datasets described previously. Parametrization 1 is the one described in the previous section and used in the previous work [START_REF] Sarrut | Modeling Complex Particles Phase Space with Gan for Monte Carlo Spect Simulations: A Proof of Concept[END_REF]. We also proposed a different way to describe the particles by replacing the exiting positions P by the so-called "ideal" positions of emission P ideal 0 within the phantom defined as follow:

P ideal 0 = P -c × t × D ( 1
)
where c is the speed of light in vacuum (we neglected the difference of the speed of light in water or other material). The computed P ideal 0 corresponds to the position in the phantom where the exiting gamma would have been emitted if it did not undergo any interaction. It means that P ideal 0 and P 0 are identical for unscattered gammas but will be different for Compton scattered gammas. The two parametrizations are strictly equivalent, there is no loss or creation of information. Exiting position can be retrieved from ideal position with the reverse equation: 2 illustrates those parametrizations. As a conclusion, the phase-space training dataset is thus still composed of gammas described with 12 dimensions: E kin , P ideal 0 or P, D, t, and 4 initial conditions: E 0 , P 0 .

P = P ideal 0 + c × t × D. The figure

Scattered gamma

Unscattered gamma Parametrization 2 Parametrization 1

Figure 2: Scheme of the phantom (in blue), the phase-space surface surrounding the phantom and the phase-space parametrization. Gamma sources within the phantom are represented by the initial energy and position of emission (E 0 , P 0 ). The 8 parameters of parametrizations 1 and 2 are shown. The "ideal" position of emission P ideal 0 is obtained via equation 1.

In experiment 1 (subsection 3.2), we will compare two different condGANs optimized with the parametrizations 1 and 2.

Conditional Wasserstein GAN with Gradient Penalty

Taking as input the previously described dataset, the goal of training the conditional GAN is to train a generator G which, given a gamma source distribution described as initial conditions (E 0 , P 0 ), is able to generate gammas following the distribution of the emitted gammas in the training dataset. This generator will hence replace the tracking of gammas inside the phantom.

Generative Adversarial Networks (GANs) is a well known concept of generative models, first introduced in [START_REF] Goodfellow | Generative Adversarial Networks[END_REF]. A GAN is composed of two adversarial networks, a generator G and a discriminator D, competing against each other. The goal of this adversarial training is to train a generator G to generate samples distributed similarly to the training data distribution. The generator takes as input a random noise vector z, sampled from an uniform or normal distribution, and produces a sample G(z) similar to real sample. The discriminator is trained to distinguish between the distributions of generated samples G(z) (P f ake ) and real samples x (P real ). It takes as input a real sample x or a fake sample G(z) probability of the sample to be real or fake (generated by G). The generator and discriminator are trained alternatively. During the training process, the generator improves using the feedback provided by the discriminator. However, the training process of a GAN tends to be unstable and hard to monitor.

Wasserstein GAN with Gradient Penalty (WGAN-GP) An alternative to GAN, called Wasserstein GAN with gradient penalty (WGAN-GP) [19; 28] was recently proposed and has several advantages. The training is, in general, more stable, leading to better quality gradients to train the generator and it does not suffer from mode collapse. It also provides a meaningful loss whose convergence correlates with an improvement in the generated samples, which makes the training much easy to monitor. In a Wasserstein GAN (WGAN), the Wasserstein metric is used to measure the similarity between the distribution of generated samples G(z) (P f ake ) and real samples x (P real ). The Wasserstein metric is an approximation of the Earth's Mover Distance (EMD), which measures the distance between two distributions. EMD also defines the cost of the optimal transport needed for moving a distribution (P f ake ) onto a target distribution (P real ). With this new cost function, the discriminator is no more bounded between 0 and 1, and instead of a sigmoid activation function in the output layer, a linear activation function is used. As the discriminator no longer classifies samples being real or fake, it is now called a critic C, which has to be 1-Lipschitz continuous (the norm of its gradients is at most 1 everywhere). Many strategies where proposed to enforce 1-Lipschitz continuity, such as weight clipping [START_REF] Arjovsky | Wasserstein GAN[END_REF] or gradient penalty (GP) [28; 30].

Conditional GAN (condGAN)

In a conditional WGAN-GP (condGAN) architecture, the vector of conditions is appended to the noise z, to the fake samples G(z) and to the real samples x. We note y=(E 0 , P 0 ) the four dimensional vector containing the initial conditions and x the 8 dimensional real sample from the Monte Carlo phase-space file. The losses of the critic and the generator can be written as:

L C = E[C(G(z|y))] -E[C(x|y)] + λE[( ∇ x|y C(x|y) 2 -1) 2 ] GP regularization (2)
L G = -E[C(G(z|y))] (3)
where C(G(z|y)) [C(x|y)] is the critic's score on generated [real] samples, E is the expected value and λ is a hyper-parameter used to scale the strength of the GP regularization. The critic outputs any real number, which represents a score that can be interpreted as how real a sample is. x corresponds to a mix of real and generated samples:

x = (x|y) + (1 -)(G(z|y)) (4)
The mix term in equation 4 ensures that the 1-Lipschitz constraint is enforced by sampling on straight lines between pairs of generated and real samples. The random sampling is obtained by sampling from a uniform distribution.

Training the condGAN

The training of the condGAN alternates the interdependent training of the generator G and the critic C into 2 steps. In step 1, G generates a fake sample G(z|y) and the critic gets both fake and real x|y samples as input. In this step, the generator's parameters are fixed, only the critic's parameters are updated by minimizing the loss in equation 2, as the critic is trying to maximize its score on real sample C(x|y) and minimize its score on generated sample C(G(z|y)). In step 2, G generates a fake sample and the critic gets only fake samples. The critic's parameters are fixed, only the generator's parameters are updated by minimizing loss in equation 3, as the generator tries to maximize the critic's score on fake samples. At the beginning of the training, the critic's loss in equation 2 should be negative, as the critic's score on real samples should be higher than the critic's score on fake samples. As the critic's loss converges toward 0, we should observe an improvement of the generated samples, because it means the critic's score on generated samples is getting closer to the critic's score on real samples.

Evaluation

In order to evaluate the proposed method with the NEMA IEC phantom, we considered the image representing the spatial distribution of ideal positions of emission, considering images of 128 3 voxels with size 2 3 mm 3 . This image was obtained considering the time at which the gammas reached the surface surrounding the phantom, along with the 3D direction D and position P. Taking the speed of light in vacuum, the ideal position of emission of the gammas within the phantom was obtained with equation 1 (see figure 2). This method is useful as it is very fast to compute and does not depend on the reconstruction method. It can be seen as a worst case scenario. Notice that the absorbed gammas with kinetic energy equal to 0 are removed before computing the image. In the image representing the spatial distribution of ideal positions of emission, the calibration factor (used to convert a number of counts into an activity concentration) is unknown, but identical for the condGAN and the reference Monte Carlo simulation. Thus, we considered a relative activity recovery (rAR, rAR = A GAN /A M C ) defined as the activity ratio between the condGAN (A GAN ) and the reference Monte Carlo simulation (A M C ). The activity (A GAN and A M C ) was estimated by summing the count values of voxels in the considered volume (sphere or cylinder).

For the evaluation with the 3D CT image of a patient and 3D activity images, we considered the simulation of a complete SPECT acquisition. The considered SPECT system was the imaging head of the GE Discovery 670 with NaI(Tl) crystal. The real camera is composed of two heads but four heads were considered here to speed up the simulation. The collimator used for 99m Tc was Low Energy High-Resolution (LEHR) with parallel-hole. The collimator holes diameter was 1.5 mm with a septal thickness of 0.2 mm, a length of 35 mm and the crystal thickness was 9.525 mm. The effect of the digitizer chain was modeled by applying an energy resolution of 10 %.

Experiments

Three experiments were conducted: we considered the NEMA IEC phantom in the two first experiments, and the 3D CT image in the last experiment. In experiment 1, we evaluated the differences between the two proposed parametrizations. In experiment 2, we evaluated the conditional generation capabilities of the trained condGAN. We considered several reference phase-space files, and compared them with the condGAN-generated phase-space files using as inputs the initial conditions of the reference phase-space files. In experiment 3, we evaluated the performance of the condGAN in the simulation of a SPECT acquisition and reconstruction. Before describing the experiments, we first detail the condGAN training parameters.

Detailed parameters of the condGAN training

The NEMA IEC phantom (left image of figure 1) is composed of 6 spheres of radius 10, 13, 17, 22, 28 and 37 mm, with outer shells made of plastic and filled with water, a central cylinder of radius 2.5 cm and height 204 cm composed of the lung ICRP material and with outer shell made of plastic. The first training dataset was generated with an uniform activity concentration of 0.06 kBq/mL everywhere during 1 second, corresponding to 3.3 × 10 7 emitted and stored gammas (the absorbed events are also stored). The second training dataset was generated with an uniform activity of 10 MBq everywhere inside the countour of the 3D CT image of the patient (right image of figure 1) corresponding to ∼ 10 7 emitted and stored gammas. CT and activity images were resampled to 4 3 mm 3 voxel size.

For both training datasets (and parametrizations 1 and 2), the condGAN was trained with batches of 2 × 10 4 gammas, for 10 5 iterations (59 epochs for the first training dataset and 199 for the second training dataset). For the critic, we used 2 hidden layers with 500 neurons and LeakyReLU activation. The weight of the gradient penalty λ (equation 2) was set to 10, as recommended in [START_REF] Gulrajani | Improved Training of Wasserstein GANs[END_REF], and the critic was updated 4 times for 1 generator update. For the generator, we used 3 hidden layers with 500 neurons with LeakyReLU activation and the dimension of the noise vector z was set to 15. The learning rate for the generator and the critic was set to 4 × 10 -4 . The inputs from the two training datasets and the initial conditions from the reference phase-space files were normalized (mean to 0 and standard deviation of 1), and the condGAN-generated phase-space distributions were denormalized. PyTorch [START_REF] Paszke | PyTorch: An Imperative Style, High-Performance Deep Learning Library[END_REF] with CUDA GPU acceleration was used. All Monte Carlo simulations were performed with Gate version 9.2 [START_REF] Sarrut | A review of the use and potential of the GATE Monte Carlo simulation code for radiation therapy and dosimetry applications: GATE for dosimetry[END_REF], using Geant4 version 11 [START_REF] Allison | Recent developments in Geant4[END_REF]. All experiments used the Geant4 physics list "standard electromagnetic option 4". Production cuts were set to 1 mm. Computations were performed on an Intel® Core™ i9-10900K CPU @ 3.70GHz with NVIDIA Quadro RTX (4000/PCIe/SSE2, 8 GB memory) and on the Jean Zay CNRS computing center (IDRIS, GENCI, Orsay, France).

Experiment 1: comparing the parametrizations

For the reference phase-space file of experiment 1, the gamma source distribution was located in the 6 spheres, the cylinder and with background activity (left image in figure 1). Gammas were emitted with 8 energies in the range 113-637 keV (see table 1). The activity concentration in the spheres and in the cylinder was 50 times higher than the background, for an activity of 1 kBq/mL. A total of 2.7 × 10 7 gammas were emitted and stored in the phase-space file.

Using the generators of two condGANs trained with parametrizations 1 and 2, we generated two phase-space files containing 2.7 × 10 7 gammas and compared the marginal distributions of the 8 parameters (E kin , P, D, t). We also compared correlation matrices, the proportions of absorbed and Compton scattered gammas at each energy, the spatial distribution of ideal positions of emission and the relative activity recovery. Note that Rayleigh scattering was included in the Monte Carlo simulation, but was not included in the analysis because it is negligible for SPECT imaging.

Experiment 2: conditional generation

In experiment 2, we illustrate that the condGAN can generate gammas corresponding to different gamma source distributions than the one used for training. Three tests were conducted. We first consider that only 3 of the six spheres (10, 17 and 37 mm) and the central cylinder were filled with activity. The activity was set to 50 kBq/mL with the principal gamma emission lines of 177 Lu, taking into account the branching ratios of the emitted gammas: 6.2% of 112.9 keV and 10.3% of 208.4 keV [START_REF]Decay radiation database[END_REF]. The activity concentration in the three spheres was 3 times higher than the one in the cylinder, and 20 times higher than the background (left image of figure 3). A total of 3.4 × 10 7 gammas were emitted and stored in the phase-space file.

For the second test, the activity was set to 10 kBq/mL with the principal gamma emission lines corresponding to 131 I: 6.1% of 284.3 keV, 82% of 364.5 keV and 7.2% of 637 keV [START_REF]Decay radiation database[END_REF]. The activity concentration in the 10 mm sphere was respectively 2, 4, 6, 8 and 10 times higher than the spheres of 13, 17, 22, 28 and 37 mm diameter, and 20 times higher than the background (middle image of figure 3). No activity was set within the cylinder. A total of 3.3×10 7 gammas were emitted and stored in the phase-space file.

For the third test, we considered six spherical activity sources of 12, 20, 30, 35, 40 and 45 mm diameter (red spheres in the right image of figure 3). The activity was set to 2 kBq/mL with 140.5 keV gammas ( 99m Tc). The activity concentration in the spheres and cylinder was 10 times higher than in the background. A total of 1.8 × 10 7 gammas were emitted and stored in the phase-space file. The activity spheres do not correspond to realistic sources and cannot be reproduced with a real IEC NEMA phantom, they are used to evaluate the genericity of the proposed condGAN. Left image: activity concentration in the 3 spheres (orange) was 3 times higher than in the cylinder (green), and 20 times higher than the background (blue). Middle image: the activity concentration in the 10 mm sphere was respectively 2, 4, 6, 8 and 10 times higher than the spheres of 13, 17, 22, 28 and 37 mm diameter, and 20 times higher than the background. Right image: the activity concentration in the 6 artificial spherical sources of 12, 20, 30, 35, 40 and 45 mm diameter (red) and cylinder was 10 times higher than the background. White areas correspond to the absence of gamma sources.

In order to evaluate the conditional generation properties of the proposed method, we generated 3 phase-space files with the condGAN trained with parametrization 2, containing 3.4 × 10 7 , 3.3 × 10 7 and 1.66 × 10 7 gammas. We compared the proportions of absorbed and Compton scattered gammas with respect to the initial energy E 0 , the spatial distribution of ideal positions of emission and the relative activity recovery.

Experiment 3: simulation of a complete SPECT acquisition

In experiment 3, we evaluated the condGAN trained with an uniform activity source in a 3D patient CT image (right image of figure 1). We considered an artificial 99m Tc source of activity composed of a large background area and 3 spheres of 52, 36 and 30 mm diameter with respective activity of 50, 20 and 40 times higher than the background. The sources were positioned in different regions (lung and soft tissues) in the CT image to obtain various attenuation conditions. CT and activity images were resampled to 4 3 mm 3 voxel size. From this activity source, we first generated a reference phase-space file containing 10 7 stored gammas and compared it with the condGAN generated gammas using the conditional information of this phase-space file, in order to evaluate the proportions of absorbed and Compton scattered gammas. Then, two complete SPECT acquisitions were simulated: the reference one with the CT image and the artificial source described previously, the second one with the condGAN. In both acquisitions, 120 SPECT projections over 360°were generated with 1 × 10 9 gammas per rotation angle. The four heads acquired incoming gammas for 30 seconds in order to create projection images with two channels for the peak energy window (126-154.55 keV) and the scatter energy window (114-126 keV). Gantry rotation was performed with a constant distance of 40 cm between the rotation center and the detector. The SPECT images were reconstructed with the RTK software [34; 35] using OSEM algorithm with 10 iterations, 15 subsets and 4.418 mm voxel size. Scatter correction was taken into account through the Double Energy Window method [START_REF] Koral | SPECT Dual-Energy-Window Compton Correction: Scatter Multiplier Required for Quantification[END_REF], attenuation and point spread function correction were performed through the method described in [START_REF] Zeng | A projector/backprojector with slice-to-slice blurring for efficient three-dimensional scatter modeling[END_REF]. We also computed the relative activity recovery for the 3 spheres, between the SPECT images reconstructed from projections obtained with the reference Monte Carlo simulation and from projections obtained with condGAN generated events. We only considered the condGAN trained with parametrization 2.

Results

Experiment 1

Figure 4 depicts the evolution of the critic's loss as a function of the epochs during the training of the condGAN with the first training dataset described in subsection 2.2 with parametrization 1. We observe the expected behavior described in subsection 2.5: at the beginning of the training, the critic's loss is negative and slowly converges towards 0. The top eight and bottom eight histograms of figure 5 depict the marginal distributions for the eight dimensions (E kin , P, D, t), obtained from the reference phase-space file of experiment 1 (blue) and generated by the generator G of both condGAN (green with parametrization 1 and orange with parametrization 2). Figures 7 shows slice of the images that represent the spatial distribution of ideal positions of emission from the reference phase-space file (blue) and from the condGAN-generated phasespace with parametrizations 1 (green) and 2 (orange). The horizontal and vertical lines are used for the profiles in the bottom panels. Despite the good agreement of the marginal distributions shown in the top eight histograms of figure 5, there is a visible discrepancy between the image obtained from the reference phase-space file (blue) and generated with the condGAN with parametrization 1 (green). In the top middle panel of figure 7, the image appears blurry and the three smallest spheres (10, 13 and 17 mm) are not clearly visible. It is also confirmed in the profiles. In contrast, the parametrization 2 (orange) leads to better agreement, clearly visible both on images and profiles. Quantitatively, table 2 depicts the relative activity recovery for both parametrizations. The relative activity recovery for the condGAN-generated phase-space in parametrization 1 is much lower than in parametrization 2, especially for the small spheres. Indeed, in the rest of this paper, we will only consider the condGAN trained with parametrization 2.

Table 2: Relative activity recovery computed between the reference phase-space file of experiment 1 and the condGANgenerated phase-space for parametrizations 1 and 2. The table 3 depicts detailed information about the gamma distributions in the bottom eight histograms of figure 5. Indeed, the distributions of the kinetic energy was in good agreement for the Compton scattered gammas, but we noticed some differences for the unscattered and absorbed gammas corresponding to the peaks. We observed, as expected, more Compton scattered gammas at low energies (between 113 and 364 keV) than at higher energy (637 keV peak). We considered that all gammas with kinetic energies below 15 keV were absorbed. To separate the Compton scattered gammas from the unscattered, we defined a threshold for each energy, and all the gammas with kinetic energy lower than the threshold but higher than 15 keV were considered as Compton scattered gammas. The threshold, the mean and the standard deviation for the 8 energy peaks are shown in the table 3. Overall, the condGAN manages to reproduce the main energy peaks within ≈ 4 keV, except for two peaks corresponding to the principal gamma emission lines of 131 I (284.3 and 637 keV).

Object

Table 3: Threshold, mean and standard deviation for the 8 peaks in kinetic energy distribution (top left panel in the bottom eight histograms of figure 5) for the gammas generated with the condGAN with parametrization 2. The gammas with kinetic energy lower than the threshold but higher than 15 keV were considered as Compton scattered gammas, while gammas with kinetic energies below 15 keV were considered as absorbed gammas.

Energy peak (keV) 112.9 140. The proportions of absorbed and Compton scattered gammas at each energies are summarized in figure 8 for both reference phase-space file (blue) and condGAN-generated phase-space with parametrization 2 (orange). For the reference phase-space file, the proportions were computed considering that gammas with zero kinetic energy were absorbed, and gammas with kinetic energy lower than the initial energy were scattered. For the condGAN-generated phase-space with parametrization 2, the proportions were estimated using the thresholds shown in table 3.

As expected, the proportions decrease when the initial energy increases. For the reference phase-space, the proportions were computed considering that gammas with zero kinetic energy were absorbed, and gammas with kinetic energy lower than the initial energy were scattered. For the condGAN, the proportions were computed using the thresholds shown in table 3.

Experiment 2

The proportions of Compton scattered and absorbed gammas with respect to the initial energy for the three reference and condGAN-generated (with parametrization 2) phase-space files of experiment 2 are shown in table 4.

Table 4: Proportions (%) of Compton scattered and absorbed gammas with respect to the initial energy for the first (112.9 and 208.4 keV), second (284.4, 364.5 and 637 keV) and third (140.5 keV) reference and condGAN-generated (with parametrization 2) phase-space files of experiment 2. The proportions were computed using the thresholds shown in table 3. 60.9 61.16 10.9 7.9

Compton

Figures 9, 10 and 11 show a slice of the images representing the spatial distribution of ideal positions of emission from the 3 reference phase-space (blue) and the condGAN-generated phase-spaces with parametrization 2 (orange). The horizontal and vertical lines are used for the profiles in the bottom panels. In tables 5, 6 and 7, we show the relative activity recovery computed between the 3 reference phase-space files and the condGAN-generated phase-spaces. 5.

Table 5: Relative activity recovery computed between the first reference phase-space file of experiment 2 and the condGAN-generated phase-space with parametrization 2. Images and profiles are shown in figure 9. Proportions of Compton scattered and absorbed gammas at the considered energies (112.9 and 208.4 keV) are shown in table 4. The proportions of Compton scattered and [absorbed] gammas for the reference and condGANgenerated phase-space files of experiment 3 were very close, respectively 66.4% [20.1%] and 67.1% [20.9%]. Figure 12 shows CT slices with reconstructed events from the reference Monte Carlo simulation and from the condGAN (with parametrization 2). The relative activity recovery (computed between the images reconstructed from the reference Monte Carlo simulation and from condGAN generated events) for the 52, 36 and 30 mm spheres was respectively 1.03, 0.98 and 0.96.

Object

Discussion and prospects

Training a condGAN. The training of the condGAN depends on a large number of parameters. We observed that gradient penalty was necessary, we tested several penalty weights and obtained the best results with λ = 10. The number of layers did not have a large influence on the results. We also obtained better results with 4 updates of the critic, and the learning rate of the generator and the critic was set to 4 × 10 -4 . We noticed that a large batch size (2 × 10 4 ) is required to learn the 8 dimensional distribution. The set of hyperparameters used in this work has been chosen according to the literature and leads to good results, but further exploration of the hyperparameter space could lead to better results. Notice also that slight differences were observed between two training datasets. The influence of the size of the training dataset was also important. The size of the two low statistics training datasets was empirically chosen in order to speed up the training while preserving the performance. For example, the training datasets described in subsection 3.1 were obtained in about 10 minutes using GATE. Additional experiments were performed with larger datasets but no improvement was obtained. The training time using GPU was about 6 hours for 10 5 iterations, longer training or larger batch size did not significantly improved the results. We observed that with larger dataset (∼ 10 8 gammas), the accuracy of the condGAN appears unchanged.

For the NEMA IEC phantom, we also tested a training dataset with non uniform gamma sources (not reported in this paper), with an activity concentration much higher in the spheres and cylinder than in the background, and similar statistics (∼ 3 × 10 7 gammas). This led to similar results compared to those presented in experiments 1 and 2. We also tested a training dataset with only one initial energy (140.5 keV, corresponding to 99m Tc). In that case, a better modeling of the peaks was observed (at 0 and 140.5 keV for the absorbed and unscattered gammas) in the kinetic energy distribution.

Parametrization. We showed that the parametrization of the position P has a strong impact on the quality of the image representing the spatial distribution of ideal positions of emission for a phase-space generated with the condGAN. Indeed, the parametrization 2 leads to much better agreement between images and profiles, as well as better relative activity recovery for the spheres and cylinder. We conjecture that this is related to a more homogeneous uncertainty in the GAN predicted gammas, as the condGAN learns the correlation between P ideal 0 and P 0 . This leads to a lower uncertainty in the predicted positions P and thus a better estimation of the spatial distribution of ideal positions of emission. Further experiments are needed to better understand it.

Modeling families of particle phase-space distributions. In experiment 1, using the condGAN trained with parametrization 2, we were able to generate gammas emitted from activity distributions within a NEMA IEC phantom, with good agreement with the reference phase-space file, for marginal distributions, correlation matrices, images of the spatial distribution of ideal position of emission, and relative activity recovery.

Conditional generation. In experiment 2, we were able to generate gammas emitted from gamma sources with various initial energies and position of emissions within the NEMA IEC phantom. We obtained a good agreement with the 3 reference phase-space files for images of the spatial distribution of ideal position of emission, relative activity recovery and the ratio of Compton scattered and absorbed gammas. Especially, in the third test of experiment 2, we considered positions and sizes of gammas sources that were not in the expected volumes in the training dataset, illustrating the genericity of the proposed condGAN. In experiment 3, we generated gammas emitted from sources with various initial position of emissions within the 3D CT image of a patient. We obtained a good agreement with the reference Monte Carlo simulation for reconstructed SPECT images, relative activity recovery and the ratio of Compton scattered and absorbed gammas.

Computational speed-up. The computational speed-up compared to Monte Carlo simulations depends on several considerations as detailed in [START_REF] Sarrut | Modeling Complex Particles Phase Space with Gan for Monte Carlo Spect Simulations: A Proof of Concept[END_REF]. Here, 8 × 10 5 particles per second (PPS) can be generated with the condGAN, while it is, in Monte Carlo simulation, 3 × 10 4 and 7 × 10 3 PPS, respectively for an analytic or a voxelized phantom. In voxelized phantom, the tracking time within the phantom dramatically increases as the voxel's size decreases (PPS close to 2 × 10 3 for 2mm size voxel). The GAN speedup is thus larger for complex simulations. The condGAN model was implemented in GATE, is open-source, and will be available in the next release. Once trained, the model can be used in two modes. In the first mode, the condGAN generates a NumPy array [START_REF] Harris | Array programming with NumPy[END_REF] of particles, outside GATE. The generated NumPy arrays can be written on disk (.npy standard binary file format) or processed on the fly by a python script, e.g. sent to ARF (Angular Response Function) detector modeling, as shown in [START_REF] Sarrut | Modeling Complex Particles Phase Space with Gan for Monte Carlo Spect Simulations: A Proof of Concept[END_REF]. The second mode directly operates in GATE, using the condGAN as a source, generating new Geant4 particles from the gammas properties (E kin , P, D, t). In this second mode, additional time is needed by the Geant4 engine to create the particles, which limits the computational speed-up. Indeed, in the simulation of a SPECT acquisition in experiment 3, the condGAN was used directly as a gamma source in GATE, PPS was around 1.2 × 10 5 using GPU and around 8 × 10 4 using CPU, against ∼ 7 × 10 3 PPS for the voxelized CT and activity images.

Comparison with other methods. Several other methods have been proposed to improve simulation efficiency for SPECT devices. Indeed, ARF [START_REF] Song | Fast modelling of the collimator-detector response in Monte Carlo simulation of SPECT imaging using the angular response function[END_REF][START_REF] Descourt | Implementation of angular response function modeling in SPECT simulations with GATE[END_REF][START_REF] Sarrut | Learning SPECT detector angular response function with neural network for accelerating Monte-Carlo simulations[END_REF] (Angular Response Function) models detector response according to gamma angles and energies, either by look-up-tables or neural networks. This method provides about ×20 speedup, by assigning to every gamma reaching the detector the probability to be detected in a given energy window. In the (Fixed) Force Detection (FD) methods [START_REF] Beenhouwer | Acceleration of GATE SPECT simulations[END_REF][START_REF] Beck | Analysis of spect including scatter and attenuation using sophisticated monte carlo modeling methods[END_REF][START_REF] Ljungberg | Attenuation and scatter correction in SPECT for sources in a nonhomogeneous object: a Monte Carlo study[END_REF][START_REF] Cajgfinger | Fixed forced detection for fast SPECT Monte-Carlo simulation[END_REF], the gamma is directed towards the detector for every interaction and is weighted according to the probability that a gamma with such direction would exist and is detected. This method does not require dividing the simulation into two parts. Reported speed-up varies but can reach ×100. However, for a complete 360 • tomographic simulation, still a large number of particles must be tracked within the phantom. It would be interesting to combine FD with GAN generated particles in order to potentially synergize the speedups. Also, it is envisioned to investigate the possibility for the condGAN to generate gammas with their corresponding likelihood of being detected and participate to the image formation. The computational speed-up of our approach is still far from GPU-based Monte Carlo code capable of generating 3200 × 10 6 gammas per seconds [START_REF] Rydén | Fast gpu-based monte carlo code for spect/ct reconstructions generates improved 177lu images[END_REF], combined with FD methods and adjusted cross-sections. But again, GAN approach may also be integrated in a GPU framework. It should be emphasised that the proposed approach is more general and can also be useful for other simulations, not only for SPECT systems.

Limitations. Without the condGAN modeling, the phase-space approach is practically intractable because billions of gammas would be needed in the file. Training the condGAN still requires an initial phase-space but of limited size (few tens of millions of particles, around 3 GB). As shown in the results of experiment 1, the condGAN did not perfectly reproduce the kinetic energy distribution. Indeed, the energy peaks (corresponding to unscattered gammas) were spread around the value of the initial energies (E 0 ). However, the standard deviations were lower than 4 keV and is sufficient regarding most of current SPECT acquisition systems. In experiments 1 and 2, the condGAN did not perfectly reproduce the ratio of absorbed and Compton scattered gammas. This means that the condGAN may slightly over or under estimate the attenuation at each energy, leading to an uncertainty on the estimated activity, especially in small spheres. This may explain why, in tables 2 and 6, the relative activity recovery was lower for the smaller spheres, especially the 10 mm sphere. Indeed in that case, several initial energies were considered with an overestimation of the proportion of absorbed gammas (see figure 8 andtable 4). It should be however noted that the evaluation method using the image representing the spatial distribution of ideal positions of emission can be viewed as a worst case bound on the quality of the results. Indeed, real SPECT reconstruction cannot use the time information and will depict a much higher uncertainty. Finally, while being generic for a family of source distributions, the condGAN is still specific for a given phantom. A different condGAN must be trained for each phantom. In a future work, we will investigate to extend the proposed method so that the condGAN is trained only once for a family of phantoms.

Conclusion

In this paper, we presented a method to model families of distributions of particles exiting a phantom during Monte Carlo simulation of SPECT imaging devices. The method extended and improved the one proposed in [START_REF] Sarrut | Modeling Complex Particles Phase Space with Gan for Monte Carlo Spect Simulations: A Proof of Concept[END_REF] using a conditional GAN (condGAN) and an optimized parametrization. Our results show that with a condGAN trained only once for a given phantom, particles with various energies and positions of emission within the phantom can be generated, illustrating the genericity of the proposed condGAN. The generation is fast, around 8 × 10 5 particles per second, and can be used to speed up the development of SPECT imaging systems. In the condGAN-generated particle distributions, sharp features, such as the energy peaks, are still imperfectly reproduced. The proposed concept is shown here for SPECT imaging with gamma emission but could, in principle, be extended to other types of simulations.

Figure 1 :

 1 Figure 1: Left: Scheme of the NEMA IEC phantom composed of a water container (blue), 6 spheres of 10, 13, 17, 22, 28 and 37 mm diameter (orange) with outer shells made of plastic and filled with water, and a central cylinder (green) of radius 2.5 cm and height 204 cm composed of the lung ICRP material and with outer shell made of plastic. Right: Slice of the CT image with patient contour (in red).

Figure 3 :

 3 Figure3: From left to right: schemes of the gamma sources within the NEMA IEC phantom for the first (left), second (center) and third (right) tests of experiment 2. Left image: activity concentration in the 3 spheres (orange) was 3 times higher than in the cylinder (green), and 20 times higher than the background (blue). Middle image: the activity concentration in the 10 mm sphere was respectively 2, 4, 6, 8 and 10 times higher than the spheres of 13, 17, 22, 28 and 37 mm diameter, and 20 times higher than the background. Right image: the activity concentration in the 6 artificial spherical sources of 12, 20, 30, 35, 40 and 45 mm diameter (red) and cylinder was 10 times higher than the background. White areas correspond to the absence of gamma sources.

Figure 4 :

 4 Figure 4: Critic's Wasserstein loss as a function of the epochs during the training of the condGAN with parametrization 1. The evolution of the critic's loss during the training with parametrization 2 is similar.

Figure 5 :

 5 Figure5: Marginal histograms obtained from the reference phase-space file of experiment 1 obtained in Monte Carlo simulation (blue) and generated from the generator G of the condGAN with parametrization 1 (green in top eight histograms) and parametrization 2 (orange in bottom eight histograms). Input conditions (E 0 , P 0 ) were taken from the reference phase-space. The peak at 0 keV represents the absorbed gammas. Positions X, Y and Z for the condGANgenerated phase-space in parametrization 2 were obtained from the ideal positions X ideal 0 , Y ideal 0

Figure 6 Figure 6 :

 66 Figure6illustrates the correlations between the pairs of parameters for the reference phase-space file obtained in Monte Carlo simulation and generated with both condGAN parametrizations.

Figure 7 :

 7 Figure 7: Top panels: slice of the spatial distributions of ideal positions of emission from the reference phase-space file of experiment 1 (left panel) and from the condGAN-generated phase-space, with parametrizations 1 and 2 (middle and right panel). The horizontal and vertical lines are used for the profiles. Bottom panels : vertical (left panel) and horizontal (right panel) profiles extracted from the top left image (blue) and the top middle and right image (green and orange).

Figure 8 :

 8 Figure8: Evolution of the proportions of absorbed (left) and Compton scattered (right) gammas as a function of the initial energy E 0 for the reference phase-space file of experiment 1 (blue) and the one generated with the condGAN in parametrization 2 (orange). For the reference phase-space, the proportions were computed considering that gammas with zero kinetic energy were absorbed, and gammas with kinetic energy lower than the initial energy were scattered. For the condGAN, the proportions were computed using the thresholds shown in table3.

Figure 9 :

 9 Figure 9: Top panels: slice of the spatial distribution of ideal positions of emission from the first reference phase-space file of experiment 2 (left panel) and from the condGAN-generated phase-space with parametrization 2 (right panel). The horizontal and vertical lines are used for the profiles. Bottom panels : vertical (left panel) and horizontal (right panel) profiles extracted from the top left image (blue) and the top right image (orange). Relative activity recovery is shown in table 5.

Figure 10 :

 10 Figure 10: Top panels: slice of the spatial distribution of ideal positions of emission from the second reference phase-space file of experiment 2 (left panel) and from the condGAN-generated phase-space with parametrization 2 (right panel). The horizontal and vertical lines are used for the profiles. Bottom panels : vertical (left panel) and horizontal (right panel) profiles extracted from the top left image (blue) and the top right image (orange). Relative activity recovery is shown in table 6.

Figure 11 :

 11 Figure 11: Top panels: slice of the spatial distribution of ideal positions of emission from the third reference phase-space file of experiment 2 (left panel) and from the condGAN-generated phase-space with parametrization 2 (right panel). The horizontal and vertical lines are used for the profiles. Bottom panels : vertical (left panel) and horizontal (right panel) profiles extracted from the top left image (blue) and the top right image (orange). Relative activity recovery is shown in table 7.

Figure 12 :

 12 Figure 12: Reconstructed events overlaid on patient CT slices. Left: reconstruction performed with the projections obtained from the reference Monte Carlo simulation. Right: reconstruction performed with the projections obtained from condGAN generated events with parametrization 2.

Table 6 :

 6 Relative activity recovery computed between the second reference phase-space file of experiment 2 and the condGAN-generated phase-space with parametrization 2. Images and profiles are shown in figure10. Proportions of Compton scattered and absorbed gammas at the considered energies (284.4, 364.5 and 637 keV) are shown in table4.

	Object	GAN parametrization 2
	Sphere 10 mm	0.68
	Sphere 13 mm	0.78
	Sphere 17 mm	0.86
	Sphere 22 mm	0.92
	Sphere 28 mm	0.91
	Sphere 37 mm	0.95

Table 7 :

 7 Relative activity recovery computed between the third reference phase-space file of experiment 2 and the condGAN-generated phase-space with parametrization 2. Images and profiles are shown in figure11. Proportions of Compton scattered and absorbed gammas at the considered energy (140.5 keV) are shown in table4.

	Object	GAN parametrization 2
	Sphere 12 mm	0.97
	Sphere 20 mm	0.95
	Sphere 30 mm	1.01
	Sphere 35 mm	1.02
	Sphere 40 mm	0.98
	Sphere 45 mm	0.98
	Cylinder	1
	4.3 Experiment 3:

Acknowledgments

This work was performed within the framework of the MOCAMED project (ANR-20-CE45-0025), the SIRIC LYriCAN Grant INCa-INSERM-DGOS-12563, the LABEX PRIMES (ANR-11-LABX-0063) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) and the POPEYE ERA PerMed 2019 project (ANR-19-PERM-0007-04). This work was granted access to the HPC resources of IDRIS under the allocation made by GENCI (Jean Zay computing center).