PATIENT-SPECIFIC DOSE ESTIMATES IN DYNAMIC COMPUTED TOMOGRAPHY MYOCARDIAL PERFUSION EXAMINATION

Abstract The study aimed to implement realistic source models of a computed tomography (CT) scanner and Monte Carlo simulations to actual patient data and to calculate patient-specific organ and effective dose estimates for patients undergoing dynamic CT myocardial perfusion examinations. Source models including bowtie filter, tube output and x-ray spectra were determined for a dual-source Siemens Somatom Definition Flash scanner. Twenty CT angiography patient datasets were merged with a scaled International Commission on Radiological Protection (ICRP) 110 voxel phantom. Dose simulations were conducted with ImpactMC software. Effective dose estimates varied from 5.0 to 14.6 mSv for the 80 kV spectrum and from 8.9 to 24.7 mSv for the 100 kV spectrum. Significant differences in organ doses and effective doses between patients emphasise the need to use actual patient data merged with matched anthropomorphic anatomy in the dose simulations to achieve a reasonable level of accuracy in the dose estimation procedure.


INTRODUCTION
Accurate evaluation of patient-specific dose in computed tomography (CT) is a challenging task due to highly varying characteristics in patient anatomy and scan parameters. Rough order-of-magnitude estimates of effective dose are possible using doselength product and body region-specific normalised effective dose conversion coefficients, acknowledging the limitations of the traditional computed tomography dose index (CTDI) formalism. (1) Size-specific dose estimates (SSDEs) can be used to increase the patient-specific accuracy of the mean physical dose values within the scan range. (1) However, these methods do not take into account patients' unique properties and anatomy in an accurate way. In high-dose and high-risk applications, more accurate estimates of patient exposure are needed to e.g. avoid detriment in repeated procedures and as a basis for radiological optimization. Monte Carlo simulations are conducted with real patient data complemented with matched anthropomorphic phantom model data (1) to enable accurate organ and effective dose estimates at a patient-specific level.
Many Monte Carlo software packages are available that can be used to simulate patient exposure to photons and electrons. However, CT scanner-specific information is needed for MC simulations. Not all needed information is readily available, but it is possible to measure the properties of missing parameters such as x-ray spectra and bowtie filter shape. Measured CT scanner radiation dose output is required to connect simulations to real patient dose. In this study, the CT scanner information is acquired with the procedure developed and tested in an earlier study, (2) where it was demonstrated that accuracy of dose (simulation vs. measurement) to within 10% is achievable in the applied anthropomorphic phantoms.
There are several non-invasive techniques that allow the evaluation of myocardial perfusion. Myocardial perfusion can be studied with ultrasound (stress echocardiography), magnetic resonance imaging (stress CMR) and nuclear medicine imaging techniques like positron-emission tomography and single-photon emission CT. (3)(4)(5) Recent advancements make myocardial perfusion imaging increasingly feasible also in CT. (4)(5)(6) CT-myocardial perfusion imaging (CT-MPI) can be performed in two ways: static (snapshot, conventional) CT-MPI and dynamic CT-MPI. (5,7) CT-MPI examination starts with intravenous iodinated contrast administration. In static mode, a single CT dataset of the left ventricular myocardium is acquired at the estimated time of myocardial enhancement peak. In dynamic mode, many CT datasets of the left ventricular myocardium are acquired. This enables the generation of time-attenuation curves. (5) In earlier studies, mean effective dose values between 3.8 and 12.8 mSv, with an average of 9.23 mSv, have been reported for dynamic CT-MPI examinations. (8) These values are based on the dose-length product and body region-specific normalised effective dose conversion coefficients. Variability is high due to the different CT scanners and imaging protocols used. For static examinations, significantly lower doses were reported, with an average of 5.93 mSv, (8) but a large variation between 1.9 and 15.7 mSv existed.
The main purpose of this work was to implement earlier equivalent source model measurement methods and Monte Carlo simulations to real patient data merged with ICRP voxel phantom (9) to examine the spread of patient-specific dose estimations under otherwise identical CT-MPI scan protocols with the same scanner.

MATERIALS AND METHODS
Organ and effective doses were estimated for 20 patients in dynamic CT myocardial perfusion examination. To estimate organ doses outside the imaged region and effective doses more accurately, patient data were merged with ICRP 110 voxel phantom data. (9) All data preprocessing, segmentation and data analysis were conducted with MATLAB (MathWorks, Natick, MA, USA).

Patient and phantom data
Patient data comprise 20 CT angiography (CTA) axial image volumes. These 20 patients were women aged 40-84 years. Female patients were chosen in order to include highly radiation-sensitive breast tissue in the study context. Furthermore, breasts are directly irradiated organs in cardiac CT examinations. Consistent with the patient data, a reference ICRP voxel phantom model representing an adult female was used in simulations. Images were taken in HUS Medical Imaging Center with Somatom Definition Flash dual-source CT scanner (Siemens Healthcare, Forchheim, Germany). The pixel size of the data varied between patient scans, but the slice thickness was 3 mm in all cases. The pixel size of the ICRP phantom was 1.775 mm, and the slice thickness was 4.84 mm. (9) The length of the imaged area varied, yielding 34-80 images per patient. In this study, we did not have information on patient weight or height, so we used the maximum effective diameter to get an estimation of patient size. The effective diameter was calculated as follows: where LAT refers to the lateral dimension and AP refers to the anterior-posterior dimension of the  1  313  2  331  3  320  4  272  5  276  6  292  7  245  8  256  9  297  10  251  11  322  12  316  13  342  14  312  15  313  16  329  17  332  18  304  19  322  20  354 patient. (10) The effective diameter was calculated slice by slice, and the maximum value was used. The effective diameters of the patients, calculated after the missing tissue values were added, are listed in Table 1.

Data preprocessing and phantom matching
The CTA image stack consists only of a small part of the human body. The applied dosimetric patient model should include the whole body to take into account the scatter dose component outside of the scanned region. Therefore, the ICRP 110 female voxel phantom data was merged with the patient CTA data. For most of the images, the outer adipose tissue did not exist on images because of the small field of view. Missing adipose tissue was added to patients' images if it was deemed necessary. Added adipose tissue had the shape of a circle segment in axial slice. However, one image stack contains only the inner part of the body, leaving most of the muscles, adipose tissue and breast out of the field of view. In this case, the outer part of the body was replaced by the ICRP phantom data. The ICRP 110 phantom data were transformed into phantom-based Hounsfield unit (HU) values. Values of −1000 HU for air, −700 HU for lungs, −100 HU for adipose tissue, 30 HU for soft tissue, 200 HU for blood and 1000 HU for bone were used. Cortical bone, spongiosa and medullary cavity were considered to be bone tissue. During the imaging, the patients were positioned arms up. Since in the ICRP phantom the arms and hands are next to the body, it was modified to an arms up position ( Figure 1). In most CT images, the patient support table is partially visible, whereas the ICRP phantom does not include the patient table. The table was removed from the CT images as the first preprocessing step. In the second step, the patient and ICRP phantom data were resampled to have the same pixel size of 0.98 mm. During the resampling process, the slice thickness of the ICRP phantom was also resampled to 3 mm. In the third step, corresponding slices between patient data and ICRP phantom data were identified. In the fourth step, the ICRP phantom data were rescaled in such a way that cross-sectional areas of the bodies were the same at the lower junction point. In all resampling and rescaling processes, the nearestneighbour interpolation method was used. In the fifth step, patient and ICRP phantom data were merged. In the sixth step, the patient table was added to the images. In the final step, the missing adipose tissue was added to the patient images if needed. Figure 2 shows schematic axial slice images of a patient in different preprocessing steps. Data preprocessing can take from a few minutes to an hour depending on the dataset. Most of the processes are automatically done, but finding the correct slices from the ICRP phantom to intersect with the patient and phantom data and adding adipose tissue still require manual help from the user.

Segmentation
Patient organ segmentation was performed partially manually and partially based on HU values. For the ICRP phantom data, the original segmentation was used. In the manual segmentation of the CT image, the organ was delineated slice by slice. In the HU method, voxels in certain areas with predetermined value ranges were associated with that organ. The manual method and the HU-based method were used together for most of the datasets: target tissue/organ was first manually delimited roughly and then in that area the correct HU values for the particular tissue/organ were chosen. Lymphatic nodes could not be segmented from the patient data due to their small size. Therefore, only the ICRP phantom results were used to estimate absorbed dose for lymphatic nodes. This segmentation method is time-demanding due to the manual segmentation. Depending on the size of the image stack and the quality of segmentation required, it can take up to a few hours per patient.

Source models
Source models for the Somatom Definition Flash were used in this study. Bowtie filter thicknesses for tube A and tube B were acquired on site. COBRA formalism (characterisation of bowtie relative attenuation) was used to determine bowtie filter shape, and x-ray spectra were calculated using aluminium attenuation measurements together with SpekCalc software tool and iterative algorithms. (2,11,12) The bowtie attenuation characteristics with inherent filtering were measured for the 80 kV spectrum for tube A and tube B and Al-equivalent bowtie filters were calculated, and the resulting thickness is used in the simulation presented in Figure 3.

Simulations
Monte Carlo simulations were performed with ImpactMC software (CT Imaging GmbH, Erlangen, Germany). ImpactMC is designed especially for Monte Carlo simulations of the 3D dose distribution in CT diagnostics. ImpactMC needs a CT image stack and scan parameters as input. Simulations were performed for the 20 patient data cases, including the scaled and matched ICRP phantom and the regular ICRP phantom. For ICRP phantom data simulation, the blood HU value was changed to 500 HU to simulate blood with iodinated contrast agent.
To simulate shuttle mode imaging, two simulations were made per image stack, one simulation for tube A and another simulation for tube B. The resulting dose maps from the two simulations were merged to get the total absorbed dose map. In dynamic study, several images are taken and starting angle can vary between image acquisitions; therefore, two full rotations were used in the simulations to average the starting angle variations. The beam collimation was 38.4 mm.
Rotations were overlapping 10% such that the table increment was 34.56 mm (i.e. pitch = 0.9) and the total length of the imaged area was ∼73 mm. The estimated absorbed dose was multiplied by the factor 95/360 since the shuttle mode uses an actual rotation angle of only 95 • . Simulations were made for 80 and 100 kV spectra. Based on Siemens instructions, 370 mAs exposure was used for 80 kV simulations and 300 mAs exposure for 100 kV simulations. The distance from the focus to the centre of rotation was 595 mm, and the fan angle for tube A was 0.7955 rad and for tube B 0.5441 rad. Values for the Air KERMA free-in air at the iso-centre were measured (2) and used for the simulation: 5.66 mGy/(100 mAs) was applied for 80 kV simulations and 10.35 mGy/(100 mAs) for 100 kV simulations. The bowtie filter model for the 80 kV spectrum was used in the simulations. Fifteen points on the time-attenuation curve were simulated, and thus, the simulated doses were multiplied by the factor 15.
In each simulation, 15 × 10 9 photons were simulated. Large numbers of photons were simulated to obtain better statistical accuracy in the low dose area. The values for the number of interactions and the cut-off energy were set to 10 and 10 keV, respectively. The number of projections per rotation was 36. Simulation parameters were chosen based on the ImpactMC manual to have good speed and accuracy for simulation. Simulations took less than a half-hour using GPU acceleration with NVIDIA GeForce GTX 1060 6 GB GPU. Size of the image stack affects the duration of the simulation. Different iodine concentrations on blood forced to make many inputs to material conversion tables. Material conversions are presented in Table 2. The maximum blood-iodine mixture HU value was determined as follows: (1) If >10% of the blood had HU values >501 HU, the maximum blood HU value was the value containing 90% of the blood voxels.
(2) If the maximum blood HU value was <501 HU, the maximum blood HU value was the maximum value.
(3) Otherwise, the maximum blood HU value was 499 HU.
Iodine concentration in the blood was calculated based on the HU value. The HU value of blood without iodine was assumed to be 40 HU. The relationship between the iodine concentration and the contrast enhancement was taken from the previous study by Bae et al. (13) For 80 kV, 41.12 HU/(mg iodine/ml) and for 100 kV 31.74 HU/(mg iodine/ml) were used. Table 3 shows the mass percentage of blood and iodine in all blood-iodine mixtures.

Absorbed and effective doses
Absorbed doses for different tissues were calculated from simulations as a mean dose of segmented voxels. However, the absorbed dose of red bone marrow (rbm) was estimated in the patient data area using bone segmentation and in the ICRP phantom area using the spongiosa segmentation that contains rbm. The amount of spongiosa in the patient data area was estimated to be the same as in the ICRP phantom, so the number of segmented bone voxels was multiplied by a correction coefficient when calculating the absorbed dose for rbm. The mass-energy absorption coefficient method (MEAC method) was used to transform the bone dose in patient data to the rbm dose. The estimated absorbed dose for rbm was calculated as follows: where D rbm is the absorbed dose for the rbm and D bone is the absorbed dose for the bone, including all compositional structures, and μen ρ the 28 MEAC. (14) MEACs were calculated using the equation below:

PATIENT-SPECIFIC DOSE ESTIMATES
where E refers to photon energy, m refers to elemental materials, R m is the relative portion of elemental material on tissue x and R E is the relative portion of energy on spectrum. MEACs were interpolated from data with a 1 keV interpolation interval using cubic spline interpolation. Ratios of mean MEACs were calculated for the incident 80 and 100 kV spectra (i.e. before photons reach the patient) using rbm and bone elemental composition and MEACs acquired from the NIST database. (15) The rbmto-bone ratio of the mean MEACs for the 80 kV spectrum was 0.2761 and for the 100 kV spectrum 0.3014.
The effective dose was calculated as follows: where w i is the weighting factor for tissue or organ i based on the ICRP 103 publication and D i the corresponding mean absorbed dose (16) .

RESULTS
Patient-specific absorbed doses and effective doses are presented in Tables 4 and 5 for 80 and 100 kV simulations, respectively. Figure 4 shows simulated absorbed dose distribution with the 80 kV spectrum for modified ICRP phantom. Patients can be separated into two groups: those whose effective diameter is less than or equal to 32 cm and those whose effective diameter is >32 cm. The grouping is based on the Siemens instructions recommending that 80 kV tube voltage with 370 mAs exposure is suitable for patients with a thorax diameter of no >32 cm. Based on the same instructions, 100 kV tube voltage with 300 mAs is suitable for patients with a thorax diameter of >32 cm. Selected patients for the 80 kV group are therefore patient numbers 1, 3, 4, 5, 6, 7, 8, 9, 10, 12, 14, 15 and 18. The remainder comprise patients for the 100 kV group. Table 6 shows the mean effective dose and mean absorbed dose ± standard deviation for tissues/organs with mean absorbed dose >5 mGy in 100 kV simulations. Results are calculated based on all 80 and 100 kV simulations and also based on selected patients for the 80 and 100 kV groups. The highest mean absorbed doses were found for the heart and then for the lungs and breasts. However, the breasts have the highest standard deviation and the variability in breast doses is much higher than in doses in other tissues/organs. When comparing mean effective doses between 80 and 100 kV simulations in the case where all patients are included, the effective dose is 75% higher in 100 kV examinations. However, if the groups with the selected patients are compared, the effective dose is only 36% higher in 100 kV examinations. The combined mean effective dose for selected 80 and 100 kV examinations is 9.2 ± 2.5 mSv.

DISCUSSION
Patient-specific dose estimate method can be used to evaluate accurately patient dose in CT. In this method, real patient data complemented with the matched anthropomorphic phantom model is used as input volume in Monte Carlo dose simulations. Taking into account patients' unique anatomical properties, this method provides more accurate patient dose estimates at the individual level than SSDEs or dose-length product and effective dose conversion coefficients. For Monte Carlo simulations, bowtie filter shape, x-ray spectra and radiation dose output are needed. Methods to acquire these parameters have been developed and tested in an earlier study. (2) In this study, Monte Carlo simulations and previously developed methods to acquire necessary scan parameters were used to evaluate patient-specific dose estimates for 20 patients in dynamic CT myocardial perfusion examination.
Li et al. (18)(19)(20)(21) have earlier evaluated patientspecific dosimetry and applied it to paediatric patients. However, no earlier studies exist on patientspecific dosimetry in myocardial perfusion examination. Dose estimations for myocardial perfusion examinations in various studies are based on doselength product and effective dose conversion coefficients. Such conversion coefficients are only orderof-magnitude level estimates provided for rough anatomical regions such as the coefficients given in the EU RP 154 document. (22) Danad et al. (8) collected information about radiation dose in CT myocardial perfusion examinations. For dynamic studies, the effective dose varied from 3.8 to 12.8 mSv, with a mean value of 9.23 mSv. Fujita et al. (23) investigated dose reduction in dynamic CT stress myocardial perfusion imaging and compared 80 kV with 370 mAs and 100 kV with 300 mAs protocols. Estimated effective dose for the 80 kV group was 6.1 ± 1.1 mSv and for the 100 kV group 10.7 ± 1.9 mSv. Compared with our results, these effective dose estimates are lower. However, patient characteristics were different from our study, and most importantly, the average number of acquisition phases was lower in their study, which has a proportional effect on patient dose.
The breasts and lungs have the largest contribution to the effective dose. However, marked variability in breast dose exists between patients. If the breast tissue is mostly in the scanned area, the absorbed dose will be high. This can be seen in the ICRP phantom data simulation where a large amount of breast tissue is in the scanned area. However, if only a small part of the breast tissue is located within the scan area, the mean breast dose will be substantially lower. This demonstrates the consequence of variable patient anatomy on the actual radiation doses to the patient, demonstrating the need for patient-specific dose estimates also in cardiac CT perfusion imaging.
When calculating the effective dose, correct organ segmentation is very important. In most clinical CT scans, as appropriate according to clinical indication and general optimization principles, the acquired CT scan data of the patient should cover only a part of the    body. When applied to chest CT data, as in this study, or especially when applied to perfusion scan data of the heart, the missing part of the patient anatomy has to be approximated with an anthropomorphic model. The model should match the actual patient anatomy as closely as possible. In this study, the female patients were merged with the female ICRP phantom model to supplement the missing individual anatomy outside the patient chest CT scan data. Organ size, shape, orientation and location in the body are important parameters when estimating the absorbed radiation dose to organs. Individual sizes of organs vary significantly (24,25) so an anthropomorphic phantom model, even a matched one, provides only a rough estimate of the real patient anatomy. Fortunately, the primary patient scan data represent the directly exposed anatomical area (excluding the helical overranging, which can be minimised by using dynamic helical collimation), and the supplementary phantom data are mainly used to determine the scatter contribution. This makes effective dose estimate more accurate than a situation where there is only phantom data. Manual segmentation is time-consuming. If personalised dosimetry is intended for general use for all patients after CT examinations, segmentations need to be done automatically. Deep learning and convolutional neural network might be useful for this purpose. (26,27) As stated by Dong et al., (27) automatic segmentation with a deep-learning method can be performed for multiple organs accurately within a few seconds, whereas manual segmentation usually requires 30 minutes or more.
In this study, source models and patient data from Somatom Definition Flash were used, which limits effective dose estimates to that device alone. However, the method is generally applicable to any CT scanner and any scanned region, not limited to myocardial perfusion imaging.

Uncertainties
There are several sources of uncertainty in the simulated doses. Sources of uncertainties are related to the simulation software, equivalent source models, input to material conversion, patient data and voxel phantom data. In simple phantom simulations with ImpactMC, when the Al-equivalent bowtie filter and spectra are determined and applied, as in this study, the maximum difference between measurement and simulation results has been ∼15% and the typical difference some 5%. This level of accuracy is achievable for simple cases with a medium-sized thorax CT phantom (CIRS 007TE-17). Simulations with patient data are more complex, which can lead to higher uncertainties. Each individual set of patient data contains a large number of different tissues. Due to the basic elemental material properties and physical interactions, it is anticipated and also observed that the CT-number (HU) contrast scale of the tissues can overlap significantly. This is also reflected in the application of the material conversion table in the simulation. For studies where iodine contrast agent is used, the overlap of iodine contrastenhanced blood and bones with calcium presents an overlap in the CT-number (HU) range. Therefore, some bone tissue with a low HU value was considered to be enhanced blood, and also some high HU value blood was considered to be bone tissue.
Patient CT image data may also incorporate different types of artefacts. Artefacts change HU values in an undesirable way, which may lead to errors in material conversion. If the artefacts are significant, HU values might need further correction in order to reach accurate results in dose simulations. Fortunately, the patient data in this study had only small metal artefacts, which had a minor effect on the results.
Obviously, the segmentation of patient data needs to be done accurately, particularly for organs that are only partially in the primary x-ray beam. Unfortunately, in this study, lymphatic nodes could not be segmented from the patient data. That is why the ICRP phantom simulations were used to estimate absorbed dose of lymphatic nodes. However, lymphatic nodes are distributed over the body, with only a small mass fraction of nodes in the imaged region. Therefore, lymphatic nodes have only a minor effect on the effective dose. Furthermore, rbm could not be segmented from CT images. Therefore, the MEAC method was used to estimate the absorbed dose.
The scaling method used to scale the ICRP phantom to the same size as the patient was simple. The method was selected because it was easy to implement. However, it has some shortcomings that should be taken into consideration when viewing the absorbed dose results for organs and tissues. The method scales all organs with the same scaling factor depending on the patient's size relative to the ICRP phantom size. Thus, it fails to take into account the patient's body composition. For a large patient, it is more probable that organs are scaled too much relative to the organs' actual size. For patients with a size nearer to that of the ICRP phantom, the scaling method most likely works better than for larger patients. Also, organ locations on the ICRP phantom are fixed, and these locations might not correspond to the actual locations on a patient's body.
The nominal x-ray beam collimation width was applied in the simulations in this study. However, the actual total collimation is slightly wider than the nominal collimation width. This can lead to an underestimation in dose estimates, especially for organs that are partly in the x-ray beam. Lack of the 100 kV Alequivalent bowtie filter shape forced us to use the 80 kV Al-equivalent bowtie filter also for 100 kV simulations. However, based on our measurements and a study by Yang et al., (17) there should be only small differences between the 80 and 100 kV Al-equivalent bowtie filter shapes. Usage of the 80 kV Al-equivalent bowtie filter shape in all simulations should therefore be an adequate solution.

CONCLUSION
Monte Carlo simulations with a realistic CT scanner model and phantom data-enhanced patient data were used to provide absorbed dose maps. Patient-specific organ doses and effective doses were determined from these dose maps. Significant differences in organ doses and effective doses between patients were found and can be explained by differences in patient data. This emphasises the need to use actual patient data merged with matched anthropomorphic anatomy in the dose simulations to achieve a reasonable level of accuracy in the dose estimation procedure.