Development of anthropomorphic mathematical phantoms for simulations of clinical cases in diagnostic nuclear medicine

ABSTRACT We seek to propose a new generation of mathematically stylised phantoms with high levels of anatomical detail for the purposes of SPECT and PET simulations. Two basic phantoms were created: a mathematical model of the torso (MMT) and a mathematical model of the whole body (MMB). The phantoms can be manipulated by varying the position, size and shape of bones and organs to simulate patients with anatomical variants. Verification of the developed phantoms was performed by comparing simulated planar images for MMT and MMB with clinical data. Both phantoms did produce realistic planar images, which closely mimic real clinical planar images. Images obtained in computer simulation of myocardial perfusion SPECT with virtual patient (MMT) demonstrate the same errors and interpretation problems that exist in clinical routine. These phantoms can be used for simulations of clinical cases to help medical experts to understand the causes of artefacts and bias in PET and SPECT images.


Introduction
Nuclear medicine imaging techniques, such as single-photon emission computed tomography (SPECT) and positron emission tomography (PET) provide non-invasive information about the functional state of living organisms. Such information enables the detection of diseases at a very early stage, even before morphological changes may appear, which is of a particular importance for patient management in oncology and cardiology. SPECT and PET images describe the distribution of radiopharmaceuticals (aka tracers) inside the subjects, and, thereby exhibit comparatively little -if any -anatomical background information.
The quality of reconstructed SPECT and PET images depends on many factors, such as data acquisition protocol, image reconstruction algorithms and patient anatomy. Inaccurate image reconstruction may result in false interpretation of the diagnostic information and subsequently, in an inaccurate diagnosis and suboptimal patient management. Clinical methods are limited in assessment of image quality. Hypothetically, the assessment of image quality can be performed by comparison of reconstructed images with 'the true image' -the distribution of radiopharmaceutical in a patient body. While such comparisons are not possible in reality, they can be addressed by using physical or mathematical anthropomorphic phantoms.
Since 1920s different physical phantoms have been used in experimental research on radiation protection and dosimetry (Ramos et al. 2017). The phantoms were built from materials, which were radiologically close to human tissues. Today, specially manufactured physical phantoms are used in different applications, mainly in dosimetry and system testing. For example, the NEMA NU2-1994 physical phantom (https://www. supertechx-ray.com/NuclearMedicine/NEMA1994PETPhantom. php) was applied as a 'reference standard' for testing and comparative analysis of images obtained by using four different SPECT/CT systems (Seret et al. 2012). Another phantom from Data Spectrum (Hillsborough, USA) was used for SPECT quality assessment (Purser et al. 2008); it was gradually modified by using attachments mimicking tissue to create three different phantoms of the small, medium and large sizes. Nonetheless, applications of physical phantoms to study the dependency of reconstructed image quality on patient anatomy are limited given the cost associated with phantoms and their potentially limited availability.
Therefore, anthropomorphic mathematical phantoms are used that can be divided into three main categories: (1) stylised phantoms, (2) voxel phantoms, and (3) phantoms based on computer graphics (Kainz et al. 2019). The stylised phantoms are created by using constructive solid geometry (CSG) modelling techniques based on simple quadratic equations. Initially, rather rough stylised phantoms were developed for dosimetry to estimate organ doses from internal radioactive sources (Fisher and Snyder 1967). In 1990s, an improved stylised MCAT (Mathematical Cardiac-Torso) phantom was created for imaging research in nuclear medicine. MCAT was applied in computer simulations of myocardial perfusion SPECT imaging, including the variation of the size of myocardial left ventricle (LV) model and body constitution parameters in accordance to the real patient anatomical data (Tsui et al. 1993;La Croix et al. 2000). The main advantage of stylised phantoms is that they can be created by using simple equations and transformed easily to simulate patients with various anatomic structures.
The second generation of phantoms, known as 'voxel phantoms', was developed by using magnetic resonance imaging (MRI) and computed tomography (CT) patient images (Zankl et al. 1988), or by using colour photographs of cryo-sectioned slices of post-mortem subjects (Xu et al. 2000). Voxelized phantoms most closely match the tasks of personalised medicine and are used for radiation treatment planning and dose estimation (Zankl et al. 2002), but are not generally applied to PET or SPECT imaging studies.
The third generation of mathematical phantoms started with the development of non-uniform rational B-spline (NURBS)-based cardiac-torso phantom (NCAT) (Segars et al. 1999). The next step in development of anthropomorphic phantoms was the generation of whole-body 'hybrid' phantoms using Boundary REPresentation (BREP) techniques including NURBS or polygon mesh surfaces. Over the last decade, many hybrid phantoms and phantom populations were developed (Segars et al. 2009;Segars et al. 2018;Lee et al. 2007;Broggio et al. 2011;Gosselin et al. 2014;Segars et al. 2013;Segars et al. 2018;Bolch et al. 2010;Park et al. 2005;Abadi et al. 2018) and new phantoms are continuously being proposed.
This article presents a new generation of mathematical stylised phantoms with higher level of anatomical realism than in the previous stylised models. The phantoms were developed for computer simulations of clinical cases to understand the causes of artefacts and biases in PET and SPECT images. Two basic mathematical stylised phantoms are presented in this work. The first is the MMT (Mathematical Model of Torso) phantom. The second is the MMB (Mathematical Model of whole Body) phantom.

Basic assumptions in development of stylised phantoms
Patients undergoing SPECT or PET imaging, are administered a radiolabeled pharmaceutical (aka tracer), which is distributed in various organs of the body with varying concentration n. In the discrete representation, the tracer distribution will be denoted by n j , the voxels being indexed by j. The nuclear decay probability is described as: where n j ðtÞ is the radiopharmaceutical concentration at a time t in the j-th voxel, λ is the decay constant, n j ð0Þ is an initial radiopharmaceutical concentration, τ ¼ 1 λ is the decay time. The radiopharmaceutical activity f j describes the number of emitted gamma-photons from the voxel j during the data acquisition time Δt. The emission process is assumed to be isotropic. f j is a Poisson random field with the mean � f j : Assuming the acquisition time is much shorter than the decay time Δt < < τ: Thus, the mean activity � f j is determined by the initial distribution of radiopharmaceutical.
Mathematical phantoms simulate mean activity distribution in a patient body. The phantoms are developed assuming that the radiopharmaceutical has already accumulated in various organs and that this distribution remains constant during the imaging procedure. The shapes and forms of different organs are created using clinical SPECT and PET images and the Atlas of Human Anatomy (Sinelnikov et al. 2007). The technique of constructive solid geometry (CSG) is applied. This technique is based on algebraic equations, which describe 2D and 3D geometric figures, such as ellipsoids, paraboloids, planes and so on. A set of suitable figures is used to create models of each organ and each bone. The complex forms were compiled by using Boolean operators.
Two examples of simulation by using CSG technique are demonstrated in Figure 1. A 3D model of the heart left ventricle (LV) is shown in Figure 1(a). The LV shape was modelled by using two truncated ellipsoids built into each other. The external ellipsoid equation is described as: and the inner ellipsoid is: Figure 1(a) shows the developed 3D myocardial left ventricle model with a defect that simulates the ischaemic lesion in the anterior apical zone. The defect is modelled by using the secant plane equation. For comparison, a clinical 3D SPECT image of the LV is shown in Figure 1 The shapes of stylised LV model and clinical LV image are similar. By varying the parameters in Equations (4) and (5) LV models of various shapes with different wall thicknesses are created. One can also add defects of different locations and sizes. Such variations of the parameters were carried out in the works (Denisova and Ansheles 2018;Denisova et al. 2019) and have demonstrated the usefulness of stylised models for studying the causes of errors in clinical images. Figure 1(c) shows the CSG technique for modelling the rib cage. Ellipse equations and logical operators are used to create each slice of a rib cage. The resulting 3D rib cage model is depicted in Figure 1(d,e). All anatomical structures are presented in a discrete (voxel) form by using standard orthogonal grids in a Cartesian coordinate system.

MMT phantom
The MMT (Mathematical Model of Torso) phantom was originally conceived in (Denisova and Terekhov 2016). This phantom is intended for simulations of myocardial perfusion SPECT/CT imaging in nuclear cardiology. Here, MMT was improved by adding additional structures and organs, such as stomach, kidneys, shoulder joints, arm bones and scapulae. Stomach, kidneys and shoulder joints models are based on ellipsoid equations; stomach and kidneys shapes were modelled as combination of two ellipsoids of different sizes. Arm bones and scapulae models were developed by using planes, parabolic and cylinder equation (Figure 2).
The phantom simulates a patient position with arms up. The heart motion effect is not included into present MMT phantom. At the moment, the MMT model includes two variants of left ventricle (LV) form in systolic and diastolic phases. These forms were developed and studied in (Denisova and Ansheles 2018). A 3D elliptical model with anatomical apical thinning was developed to simulate the LV in diastole and a 3D LV model of conical form was developed to simulate the LV in systole. In systolic phase the apex wall is much thicker than in the diastole phase. Both forms were developed on the basis of the ellipsoid equations. The myocardial defects of different size and location can be added and studied in computer simulations. Since the phantom is defined using analytical formulas, it can be voxelized. In our numerical simulations, the MMT phantom was voxelized on grids of 256 × 256 × 256, 128 × 128 × 128 and 64 × 64 × 64.
The MMT phantom was used to create both the 3D 'source function' (the distribution of radiopharmaceutical in patient organs) and 3D 'attenuation map' to take into account the attenuation of gamma radiation in biological tissues. In this paper, the attenuation map was generated as a spatial distribution of attenuation coefficients in three different attenuating  media: bones, lungs (air) and soft tissues (water). The attenuation map was calculated for gamma quanta of 140 keV emitted by Tc-99 m (Patton and Turkington 2008). Figure 3 shows the transverse sections of the MMT-based attenuation map (a) and an example of clinical CT image (b) that was generated as the attenuation map in clinical SPECT/CT study.
In order to create 'the source function', a special method based on clinical SPECT/CT data was developed. The method included the following steps. First, shapes, sizes and localisations of thoracic organs were simulated using CT data from clinical SPECT/CT scans. Next, the relative activity values in different organs were simulated. These values can be obtained from clinical SPECT images. They can be varied in the MMT model to simulate a specific clinical case. For example, the following values of relative activity 99mTc-MIBI distributions were assumed for the MMT-phantom in Figure 4 (1). In the third step, noise-free projection data were calculated for the MMT with given relative activity values. The projection data were generated in accordance with standard protocol of SPECT data acquisition procedure (patient positioning, 'step and shoot' operation mode, number of projection angles, image matrix, total number of counts) (Verberne et al. 2015;Dorbala et al. 2018). In the fourth step, Poisson noise was simulated for the noise-free data by the von Neumann rejection method (Christopher 2006). The simulations were run for the given total number of registered quanta (Figure 4(b)). The scaling factor was calculated from the relations between the calculated noise-free data and the simulated data with Poisson noise. In the next step, the given relative activity values were scaled by that factor.

MMB phantom
The second phantom, the MMB (Mathematical Model of whole Body), is developed for simulations of SPECT and PET imaging in oncology ( Figure 5). In PET and SPECT oncological images, tumours do not necessarily coincide with anatomic organ forms. The anatomical structures are typically used for modelling the 'attenuation map'. In Figure 5, the MMB simulates an 'adult patient with a normal body type'. The phantom includes the main internal organs, skeletal system and soft tissues. Particular attention has been paid to modelling the skeletal system, since bone metastases are frequently observed in oncology scans (Table 1)  The MMB phantom includes three main segments: chest organs, abdominal organs and limbs. Modelling of each organ in phantoms was performed through separate procedures. This allows the user to change shapes and sizes of 'organs' more easily. In this paper, manipulation is demonstrated by means of MMB skeletal system transformation. The skeletal system includes rib cage, clavicles, scapulae, sternum, spine, and bones of the hands, legs and pelvic. Manipulation of skeletal system is the most time-consuming procedure compared to manipulations of individual organs. Figure 6(a) shows a 3D image of MMB skeletal system of an adult patient. The skeletal system was transformed into a 3D skeletal system of a teenager by changing the radius of a rib cage and the sizes of pelvic girdle. The teenager skeletal system is shown in Figure 6(b). Planar SPECT image calculated for the anterior view of the teenager skeletal phantom is demonstrated in Figure 6(c). The transformation of the skeletal system was performed to compare with clinical planar image shown in Figure 6(d). In numerical simulations, the MMB phantom was sampled on matrices of 256 × 256 × 768 and 128 × 128 × 384.

Simulations of raw data
In this work, verification of the developed MMT and MMB phantoms is performed by comparing simulated and clinical projection (raw) data. Projection data for the virtual patient were calculated by using home-made codes. The methodology of simulating SPECT data acquisition was developed in (Denisova and Terekhov 2016). It is assumed that the spatial distribution of emitted photon densityf ðx; y; zÞ is a stochastic field with the Poisson distribution and the average local value � f ðx; y; zÞ is proportional to the local activity nðx; y; zÞ. The mean of the projection data � g (noise-free data) can be represented in the discrete form: where A ij is the system matrix describing the probability that a photon emitted in the A ij -th voxel will be registered in the j-th pixel of the detector, θ k is the projection angle. The system matrix a) b) c)  A ij includes attenuation correction and collimator-detector response. The stochastic (raw) projection data g i ðθ k Þ are simulated using a Poisson distributed pseudo random number generator. The simulated data were generated by using the MMT phantom that models the spatial distribution of 99mTc-MIBI in thoracic organs (Figure 4(a)). A comparison between the simulated and clinical raw data is shown in Figures 4(b,c). In this simulation, MMT plays the role of a virtual patient.
The MMT was cropped according to the field-of-view of the SPECT system. An attenuation map was generated to consider attenuation effect. Figure 4 Simulated SPECT data were generated by using the MMB teenager phantom shown in Figure 6(b). Whole-body SPECT/CT data acquisition procedure was simulated. Planar SPECT image was calculated for the anterior view. This calculation simulates a standard bone scintigraphy. The result of simulation is demonstrated in Figure 6(c). For comparison, clinical bone scintigraphy image of a young patient with hip osteosarcoma is shown in Figure 6(d).

Computer simulations of myocardial perfusion SPECT/CT
To apply the phantoms, challenging clinical cases were simulated. One of such cases is considered in this paper. In accordance to international medical recommendations (Verberne et al. 2015;Dorbala et al. 2018), for interpretation of myocardial perfusion SPECT images, both AC and non-AC (with attenuation correction and without attenuation correction) images should be used. However, it is not clear why non-AC images should be used, since those images (from the physical point of view) are incorrect. In order to address this question, a simulation of myocardial perfusion SPECT imaging procedure was performed. A case of a patient with the ischaemic lesion was studied by using the MMT phantom shown in Figure 7(a). The defect (ischaemic lesion) in the upper apical zone is marked by the arrow. Projection data were simulated in a 'step and shoot' operation mode, from 32 angular views over an arc extending from the right anterior oblique to the left posterior oblique. Effects of non-uniform attenuation, collimator-detector response and Poisson statistics were included in simulations. The standard Ordered Subset Expectation Maximisation (OSEM) algorithm (Hudson and Larkin 1994) was applied for image reconstruction from the generated raw data. The two imaging protocols were simulated: with attenuation correction (AC) and without attenuation correction (non-AC). AC and non-AC images were reconstructed (Figure 7(c,d)).
To estimate the severity of the defect in the reconstructed images quantitatively, we used the following method which is explained in the inset of Figure 8. The two regions were selected: the normal myocardium zone (the point 1) and the defect zone (the point 2). The white dashed line passes through the two selected zones. The intensity profiles for the exact phantom and the reconstructed AC and non-AC images along the white dashed line are shown in Figure 8. The degree of defect severity was assessed by calculating the ratio of the intensity I 2 in the defect zone to the intensity I 1 in the healthy myocardium zone: The values Δ were calculated for the exact profile and the reconstructed AC and non-AC profiles.

Results
New stylised phantoms MMT and MMB (Figures 2 and 5) were created with higher level of anatomic realism than that in the previous generations of stylised models. Simulations of planar images using MMT and MMB have demonstrated that the both phantoms did produce realistic planar images, which closely mimic real clinical planar images (Figure 4(b,c) and Figure 6(c,  d), respectively). These results allowed us to conclude that the developed phantoms can be applied in simulations of clinical cases. Myocardial perfusion SPECT imaging was simulated. The MMT phantom modelled a virtual patient in the supine position, with arms away from the field of view, and the projection data generation code played the role of a virtual gammacamera. The obtained results are demonstrated in Figure 7(c, d). Zoomed AC and non-AC left ventricle reconstructed images in a selected slice are shown. For the sake of comparison, the exact LV image in the same slice is shown in Figure 7(b). The results of simulations have shown that both AC and non-AC images demonstrate a well-visible defect. The defect size in the reconstructed AC and non-AC images is approximately the same as the defect size in the exact phantom. However, visual analysis does not yield an estimate of the significance of the reconstructed defect. In clinical practice, a perfusion defect is assessed as the perfusion intensity reduction compared to the most normal segment. In our simulations, the degree of defect severity was similarly assessed by using Equation (7). The severity of the simulated true defect in the exact phantom wasΔ ¼ 40%. The defect severity in the reconstructed images was Δ ¼ 26% for the AC-image and Δ ¼ 33% for the non-AC image. Comparing the defect severity assessments, we can formulate the results of the numerical simulation as follows: 1) the severity of the defect is more correctly estimated in non-AC image than in AC image; 2) the severity of the disease is exaggerated in AC images. These results look as the paradoxical ones, but they are in a good agreement with clinical studies (). In the editorial paper (Apostolopous and Savvopoulos 2016), the authors stated that AC is a valuable innovation, but 'findings in the anterior, anteroseptal, anterolateral wall and the apex should be interpreted with reference to non-AC rather than AC images'. In (Ansheles 2014), it was concluded that the same LV defect can be assessed as more severe on an AC-image than on a non-AC image. Thus, we can conclude that, in the numerical simulation, we obtained the same ambiguous results as those that were obtained in clinical images. Our results explain why it is non-AC images are recommended for usage in interpretations. On the other side, the profiles presented in Figure 8 imply that the solution without attenuation correction is absolutely wrong. These controversial results are discussed in the next section.

Discussion
Clinical methods are limited in studying the causes of artefacts and bias in PET and SPECT images associated with mathematical aspects of reconstruction algorithms and methods. Such studies can be performed using mathematical phantoms and computer simulations. Unfortunately, at the moment, the number of such works is very limited because of the complexity of mathematical modelling in nuclear medicine.
Voxel and hybrid phantoms describe anatomical structures in great detail, which is necessary for simulations in medical dosimetry and CT imaging. However, SPECT and PET images exhibit relatively few anatomical details. From the mathematical point of view, the tracer distribution in a patient's body is 'a source function'. For example, in cardiologic SPECT images, the source function includes only the myocardial left ventricle, but the right ventricle is not visible. In PET images, the source function, describing a tumour in an organ, does not coincide with the anatomic organ form. Accurate quantitative reconstruction of the source function is the main objective of SPECT and PET imaging. In order to bind the source function to the anatomical structures, the SPECT and PET systems are supplemented by CT or MRI. The anatomical structures of a patient play the role of 'attenuation map' including the effects of gamma radiation attenuation and scatter. In PET and SPECT, the attenuation map distinguishes between soft tissues, bones, and lungs (air) and, as such, is a simplification of the highresolution, high-contrast corresponding CT image. Therefore, more simplified phantoms can be applied to simulate PET and SPECT imaging than in simulations in medical dosimetry and CT imaging.
In the present work, we developed the stylised phantoms MMT and MMB, specifically for simulations in SPECT and PET. The MMT was used as a virtual patient in simulation of myocardial perfusion SPECT imaging. In the numerical simulation, we have obtained the same ambiguous results as those that were obtained in clinical images: for the assessment of the defect, it is helpful to refer to the non-AC images. However, analysing the profiles, presented in Figure 8, one can see that the solution without attenuation correction is absolutely wrong. These paradoxical results can be explained as follows: As was shown in our previous work (Sinelnikov et al. 2007), a decrease of intensity appears in the rounded upper apical region on both AC and non-AC images due to limitations of the reconstruction method. The decrease is due to the specific source function (LV) form and to a limited number of projection views. But in the non-AC image, this error is then masked by the increase in intensity since the rounded upper apical region is located closer to the body boundary and an attenuation effect in this zone is less pronounced. As a result, the non-AC image seems more correct for interpretation, but any interpretation based on the results of two mutually compensating errors cannot be reliable. It should be also noted that non-AC images can be applied for interpretation only in a relative form or thanks to renormalisation. Further research is needed to develop a new approach to myocardial perfusion SPECT image interpretation that should be based on AC images.
Finally, it should be noted that modern research trends in the field of high technologies are associated with the development of powerful software packages. At the moment, great efforts are directed towards developing universal software packages to simulate SPECT and PET images formation that include computer generated population of digital anthropomorphic phantoms by using computer graphics methods and simulation of data acquisition process by using Monte Carlo method. However, this approach is computer expensive, and simulations of a single case may take a lot of time (tens of hours) and huge computer resources. In the present paper, an approach is developed which includes simple stylised phantoms with improved anatomy and uses simplified method of data acquisition simulation. This approach is not universal and has its limitations, however, in computer simulations of SPECT imaging process it has demonstrated good agreement with clinical observations and, what is especially important, calculation of one variant takes minutes and significantly less computer resources than Monte Carlo calculations. We plan to make the phantom creation code open access. Any student and scientist can easily change and improve these phantoms for their tasks.

Conclusion
A new generation of stylised anthropomorphic phantoms is proposed and validated. These phantoms can help study clinically-relevant challenges in SPECT and PET imaging, and, thus, promote an understanding of the basic causes of errors and uncertainties in clinical images.

Disclosure statement
No potential conflict of interest was reported by the author(s).