Study of secondary particles resulting along 50 – 150 MeV accelerated protons in liver tissues using MCNPX and abdomen DICOM images

Purpose – This paper aims to study the resulting Brag peak and secondary particles (neutrons, photons, deuterons, alpha, helium_3, and tritons) along protons ’ path in tissue. Design/methodology/approach – MATLAB program and MCNP code were used to read abdomen Digital ImagingandCommunicationsinMedicine(DICOM)imagesandbuilda3Dphantomtoliverinpurposetostudy resultingBraggpeakandsecondaryparticles(neutrons,photons,deuterons,alpha,helium_3andtritons)alongprotons ’ path. Findings – Fromthestudy,itwasfoundthatBraggpeakvariesfroma2cmdepthwithinthetissuefor50MeV protons to a 14.2 cm depth for 150 MeV protons; in the other hand, the total deposited energy decreases from 0.656 [MeV/g]/proton, at the depth 2 cm and 50 MeV protons, to the value 0.220 [MeV/g]/proton, at the depth 14 cm and 150 MeV protons. Originality/value – As for the flow rate of secondary neutrons and photons, the flow rate of secondary neutrons takes a maximum value (peak) in the middle of the proton path, i.e. when the energy of the protons dropstothevalueof30MeV,andthismaximumvalueoftheneutronsflowrateisaccompaniedbyamaximumvalueofthephotonflowrate,asfortherestofthesecondaryparticlesproduced(alphaparticles,deuterons,electrons,tritonsandtriplehelium),theydepositmostoftheirenergylocally.


Introduction
The importance of this research came from the fact that it sheds light on a new method of radiation therapy for cancerous tumors used in some countries of the world in addition to the traditional methods of radiation therapy for these tumors (neutron therapy, photon therapy).Many research studies are focused on this method in radiotherapy such as; Niloofar et al. (2021) that simulated a proton beam in a 10 cm distance (from source-to-surface distance) was used to irradiate a hypothetical breast tumor at the MIRD phantom under conditions of variable beam radius and energy, and concluded to that the maximum scattered radiation occurs with the beam of 1 cm radius and 70 MeV energy, while the neutron spectra show that scattered radiation is predominantly at low energy; and Sayyed et al. (2014) where the Monte Carlo N-Particle eXtended (MCNPX) and a slab head phantom were used for the purpose of simulating proton therapy in brain tissue and determine the amount of energy escaped from the head by secondary neutrons and photons.It was found that for high energy proton beams, the amount of escaped energy by neutrons is almost 10 times larger than that by photons.
In most previous research studies, the mathematical phantom [Medical Internal Radiation Dose (MIRD)] was used and also they were just focused on photons and neutrons as a secondary particle resulting from protons' reaction.
In this research, we will try to study one of the main reasons behind the use of protons and charged particles in general radiation treatments, namely the Bragg peak, and calculate it accurately for a three-dimensional computer model of the liver area obtained by processing Digital Imaging and Communications in Medicine (DICOM) Images using software, like as, MATLAB program.
Damaged cells repair tissue by rebuilding themselves by mitosis.And this process is organized and tightly controlled, but for some reasons, this process gets out of control and regularity, and the cells continue to grow abnormally, forming a mass called a cancerous tumor.This tumor may appear anywhere on the body, skin surface, ear, brain, liver, stomach, intestine, blood, . . .etc. Despite the great scientific progress in oncology that occurred in recent years, cancer is considered at the forefront of serious diseases that affect societies, and it is considered a curable disease if the reasons are taken, the most important of which are early detection and treatment (Anjak, 2011;Akhtiar, 2011).
The principle of radiotherapy is based on exposing cancer cells to radiation energy to be destroyed, stunted, and prevented from multiplying and spreading.This type of treatment is accompanied by the destruction and damage of some healthy cells through which the incident radiation passes or is dispersed, and this is one of its most important disadvantages.From here started the researchers' attempts to find ways to focus the radiation dose in the tumor and reduce it as much as possible to the extent of safety on healthy cells and sensitive organs.So came the principle of proton therapy, which depends on depositing a large dose in the tumor area, and the lowest dose in healthy tissue.This is due to the Bragg curve of charged particles, which represents the change in the linear stopping power of a beam of charged particles emitted from a given source as a function of the distance traveled by the beam in the absorption medium.The Bethe-Bloch relationship shows that the charged particles deposit the most of their energy at the end of their path into matter or medium, forming a peak called the Bragg peak, according to the British scientist William Henry Bragg, the first to discover this property in 1903 for alpha particles, as shown in Figure 1.
The linear stopping power (in one dimension) is given by the Bethe-Bloch relationship according to the formula: where − dE dx : stopping power, i.e. the decrease in energy within a unit distance traveled within the material; Z t : atomic number of material; Z p : atomic number of charged particles (1 in the case of protons); m e : electron mass; < I >: mean ionizing energy; β: the ratio of the velocity of particles to the speed of light and is a dimensional parameter; C Zt : orbit correction factor and δ 2 : density correction factor.
Then the range of charged particles within the material is computed mathematically by the integral of the stopping power over the entire path of the particle beam (James, 2006): The stopping power eq.( 1) is the basis on which the MCNPX code is based in calculating the deposited energy, and thus calculating the radiation dose resulting from charged particles.AGJSR 41,1 2. Materials and method 2.1 Computed tomography CT In order to create a three-dimensional model of the liver region, imaging was performed using a TOSHIBA Asteion imaging device at X-ray tube current 5 250 mA, 120 kVp tube voltage, abdominal scan protocol "Protocol Name 5 N. ABDOMEN" and Slice Thickness 5 5 mm.This is to obtain medical images in DICOM dimensions 512 3 512 and a spatial resolution of 1.0244 pixel/mm for a 60-year-old patient, as shown in Figure 2.
Each material has a specific linear attenuation coefficient, denoted by, which expresses the probability that a photon in length unit will disappear from the material by one of the three previous effects.The linear attenuation coefficient depends on the atomic number of the material (electronic density) and on the energy of the incoming photon, and the value of the attenuation coefficient, and thus the intensity of the rays transmitted through a substance or tissue, varies according to the chemical composition of the material, and this property is used in the computerized imaging process.Mathematically, the ratio of the transmitted beam intensity to the intensity of the primary beam received on a target is expressed as: where μðcm −1 Þ: linear attenuation coefficient of absorbent material of thickness x (cm); IðxÞ: the intensity of the primary beam transmitted through the x thickness of the material in MeV/ (cm 2 .s)and I o : the intensity of the incoming primary beam MeV/(cm 2 .s).When the X-rays, penetrating across the body, fall on the array of detectors, electrical signals of different heights are generated according to the difference in the intensity of the transmitted photon beams.These signals are collected and processed to give a gradient image in the color level (gray) according to the intensity of the transmitted beam.The digital image consists of a number of two-dimensional elements, each of which is called a pixel, defined by the x and y coordinates, and each pixel contains important information (the third dimension) that represents the intensity of the effective beam that are translated as an electrical signal and then converted into a color level.Each pixel in the CT image corresponds to a volumetric component of the scanned object called a voxel.
Depending on the protocol used in the scanning process, the body is scanned according to slices of specific thickness, so that each slice covers very small elemental volumes of tissue (voxels).Each elemental volume corresponds to an average attenuation coefficient related to the thickness and structure of the elemental size.

MATLAB program (version R2016a)
The MATLAB platform was designed to solve engineering and scientific problems.This language is based on matrices, the language most capable of expressing computational mathematics (The MathWorks Inc., 1994-2019).The MATLAB program also contains a wide range of libraries for reading and processing medical images, libraries for power systems and renewable energies, libraries for signal and logic systems . . .etc., where the MATLAB program was used to read CT images, type DICOM and create a set of instructions within a file of type Text to be exported to the radiological simulation program used in this study MCNPX.Also this program can be used to obtain the weighting factors value and SOBP calculation (Niloofar et al., 2021).

Monte Carlo N_particle eXtended (MCNPX) version
Monte Carlo methods are one of the mathematical statistics and probability methods (X-5 Monte Carlo Team, 2003).These methods consist of a set of computational algorithms that depend on repeating mathematical operations from initial random values until reaching the solution of the problem.The higher the number of iterations, the more accurate the results.MCNPX is one of the most famous computer codes that adopt Monte Carlo methods in calculating.This code is built in the Fortran language, where the MCNPX code can be used in the study of shielding (Marzieh et al., 2022), energy spectrum from various sources and radiation sources (For neutrons, photons, electrons), absorbed dose, radiation dose resulting from exposure to different types of radiation (Mehdi et al., 2022a), calculate the radioactivity of materials and elements and radiation detectors (Mehdi et al., 2022b), dairy and food irradiation (Mehdi et al., 2021(Mehdi et al., , 2022b)), in addition to many other issues that can be studied using the code.

3D Model (phantom) generation algorithm
Figure 3 below shows the algorithm for building a three-dimensional model of the liver area based on medical CT images, as this algorithm includes the following steps:

Study of secondary particles
(1) Read the CT images in MATLAB using the two instructions dicomread and dicominfo.
(2) Store the values of the color levels of pixels of each slice in a square matrix (512 3 512).
(3) Using the {for} iterations and {if} comparison relationships to replace the values of the CT numbers corresponding to a specific interval with the universe number, this number links each material, used in MCNPX input file, with a corresponding voxel element.
(4) The resulting array of component numbers (512 3 512) is stored, to be filled later within the corresponding layer in the three-dimensional model, Figure 4.
(5) Create a three-dimensional model in the form of a rectangular prism divided into a number of small elemental sizes each of them with a size (0.2 3 0.2 3 1 cm3) using the MCNPX code, the length and width of this model (on the x, y axes) is determined by the dimensions of the DICOM image, while the height of this model (on the z axis) is determined by the number of slices of the scanned area, Figure 5.
(6) Repeat the previous step for each slice.
(7) Writing an algorithm using MATLAB to fill each voxel element in the 3D model in the MCNPX input file with the special universe number value, previously stored, (in step 5).
As previously shown, each voxel element in the 3D model (phantom) with dimensions (0.2 3 0.2 3 1 cm3) corresponds to a pixel specified by (x, y) coordinates in the CT image or a specific slice.Since the dimensions of each image are (512 3 512), we need, to represent each image, to (512 3 512) a voxel element in the phantom which equivalent to (102.4 3 102.4 3 1 cm3) to cover each slice.
The height of the model on the z-axis is determined by the number of slices of the scanned area, and it varies according to the scan protocol used for each case.
Figure 5 shows the phantom of the abdominal area in which the internal organs used in this study are shown depending on reading of CT DICOM imagesthe figure on the left takes into account the air between the patient and the CT device, and the figure on the right excludes the previous air layers from the model.

Protons beam characteristics
In this study, a 1mA protons beam current was used, which contains a number of particles (protons) calculated by the relationship: where N: total number of particles in the beam; t: time in (s); e: charge of the electron, in coulombs, and I: particle current in (A).
But in order to avoid the phenomenon of range straggling resulting from the difference in the range of particles with the same initial energy, where it is difficult to determine the exact extent of the Bragg peak, the study was carried out on one proton, and then the final result is multiplied by the total number of protons calculated from the relationship (5).

Maximum range and deposited energy in the liver tissue
The phantom, in Figure 5, was used to estimate the maximum range of the Bragg peak resulting from accelerated protons of  MeV and the current of 1 mA.The study was done for one proton to avoid the phenomenon of energy.

Study of secondary particles
From Figure 6 we concluded that with the increase in the protons' energy, its range (distance traveled) within the tissue increases, as the total energy deposited over the entire proton path within the tissue increases, but on the other hand, the amount of energy deposited in the Bragg peak (at the end of the proton path) decreases.This is explained by the fact that the cross-section of the nuclear interaction of the particles, which represents the probability of the particle interacting with the nucleus, is inversely proportional to protons' energy, at the low energy of the protons, the probability of interaction is high, so when a proton enters the target material with low energy, it interacts directly to deposit the largest possible amount of energy, while at higher energies, the proton travels a greater distance within the tissue to lose its energy until it reacts, depositing the rest of its energy, which is low at the end of the path.Since with the increase in the energy of the proton, its range within the tissue increases, that is, it travels a greater distance, thus as a result of scattering the nuclei of the medium along its path, it deposits a quantity of energy along its entire path until reaching the Bragg peak region to deposit the remainder of its energy suddenly, thus the greater the path of the proton within the tissue, the greater the amount of energy that proton loses in its path as a result of collision with the molecules of the medium and the amount of energy it deposits at the end of its path at the Bragg peak decreases.
However, the total energy lost by the proton along its entire path will increase with the increase in the energy of the incoming proton, Table 1 this is shown in , which shows the variations of the range, the energy deposited at Bragg peak [MeV/g] and the total energy deposited over the whole path [MeV/g] for each source proton.Figure 7 represents the change in the value of the energy deposited at the Bragg peak (end of the proton path) for each source proton, and Figure 8 represents the changes of the total energy deposited over the entire proton path for each source proton.
It should be noted that the energy deposited from all the particles can be obtained by multiplying each value by 6:241509745 3 10 12 Â Proton s Ã , the number of protons in 1mA.

Secondary particles flux resulting from the interaction of protons with liver tissue
The flux of secondary particles resulting from the interaction of protons with the nuclei and atoms of the medium represented by the liver tissue was studied over the entire path of the protons, where the flux within each cell (voxel 0.2 3 0.2 3 1 cm 3 ) of the phantom was calculated.In the beginning, we review the flux of protons with the previous energies in the liver tissue cells for each source proton, as shown in Figure 9.

Study of secondary particles
From Figure 9, it is noticed that the protons flux is greatest at the beginning of the path and then gradually decreases with the crossing of the protons within the tissue due to the gradual decrease in the energy of the protons; therefore, the possibility of these protons interacting with the nuclei of the medium increases until it suddenly disappears at the end of path.As a result of these interactions (the interaction of low-energy protons with the nuclei of the medium), a group of secondary particles (neutrons, protons, electrons, deuterons, helium, tritons and alpha particles) in addition to photons (gamma rays, and characteristic X-rays), and the production of these particles increases as the energy of protons is decreased because of scattering from the nuclei of the medium.This statement is evident through Figure 10, which represents the flux of secondary neutrons (that is, resulting from the interaction of particles with the medium) within the liver tissue as function of the depth within the tissue for each source proton, as well as Figure 11 which represents the secondary photons flux (resulting from the interaction of particles with the medium) within the liver tissue as a function of the depth within the tissue for each source proton.
From Figure 10 it is noticed that the secondary neutrons flux is low at the beginning of the path, but with increasing depth, the secondary neutrons increases, forming a peak to return and decrease again.The reason for this is that the cross-section of the interaction of neutrons production with protons, shown in Figure 12, has a threshold that starts from the energy protons are about 7 MeV, Then it increases, forming a peak at the energy of protons E 5 30 MeV, then returning and gradually decreasing with the increase in the energy of  the protons.Thus, when a beam of high-energy protons reaches the tissue, the cross-section of the neutron production reaction is low, and when the proton beam travels a distance within the tissue, it slows down and the cross-section of the neutron production gradually increases until it reaches the proton energy of 30 MeV, then the cross-section of the neutron production and thus the neutron flux is greater, from Figure 10 for protons of 50 MeV energy need a distance of approximately 1 cm until their energy becomes about 30 MeV, and at this distance, a peak of the total neutron product of the proton interaction is formed.
As for the photons (Figure 10), most of these photons result from the interaction of secondary neutrons with the nuclei of the medium in addition to those resulting from the interaction of protons and other particles, so the production rate of the photons dramatically followed to the rate of the secondary neutrons interaction.

Study of secondary particles
From Table 2, we conclude that protons, neutrons and photons have the largest value for the total flux, followed by electrons (beta particles), then alpha particles, deuterons, and finally tritons.
Table 3 shows the percentage of the total flux of secondary particles in addition to the photons excluding the protons (incoming and outgoing) for each incoming proton along for protons at energies 50,60,70,80,90,100,110,120,130,140 and 150 MeV.From Table 3, we conclude that at the lower energies of the protons E 5 50 MeV, the photons flux is the highest which is about 67.33%, while the flux of neutrons 30.93%, electrons 1.16%, alpha particles 0.27%, deuterons 0.27%, tritons 0.02% and helium 0.01%.When the energy of protons increases to 150 MeV, the percentage of the total flux of neutrons, deuterons, helium and tritons within the tissue increases to 61.77%, 1.04%, 0.04% and 0.07%, respectively, in contrast, the percentage of the alpha particles flux, electrons flux and photons flux decreases to 0.15%, 0.85% and 34.43%.
Returning to Table 2, we find that with an increase in the protons energy, the flux of all particles increases, but the ratio of neutrons, deuterons, helium and tritons flux is higher than the rest particles, meaning that the percentage of neutrons, deuterons, helium and tritons flux increases in conjunction with decreasing the flux of alpha particles, electrons and photons.
Figure 13 shows the percentage changes of secondary particles flux as a function of the energy of the incoming proton, where it is noticed that at low energies, the proportion of the The total secondary particle flux ((Particle/ cm2.s)/(proton)) for each incoming proton along the proton path for protons at energies 50,60,70,80,90,100,110,120,130,140,150 MeV Table 3.The percentage of the total flux of secondary particles in addition to the photons excluding the protons (incoming and outgoing) for each incoming proton along for protons at energies 50,60,70,80,90,100,110,120,130,140,150 MeV AGJSR 41,1 resulting photons is about 67%, while the proportion of neutrons does not exceed 31%.The matter is reversed at high energies, where the proportion of neutrons reaches 63% and the proportion of photons decreases to 34% of the total flux of secondary particles produced.

Conclusion
In this paper, a study of the range of Bragg peak in liver tissue resulting from accelerated protons of energies 50,60,70,80,90,100,110,120,130,140 and 150 [MeV] was carried out using the MCNPX program and the MATLAB program to create a three-dimensional model of liver region, where it was found that Bragg peak varies from a 2 cm depth within the tissue for the energy of incoming protons E 5 50 MeV to a depth of 14.2 cm for 150 MeV protons, while the deposited energy in Bragg region decreases from 0.656 [MeV/g.proton]at 50 MeV protons and 2 cm depth within the tissue, to 0.220 [MeV/g.proton]at 150 MeV protons and 14.2 cm depth.
The focus has been on the secondary neutrons and photons because they contribute to increase the dose of healthy tissues surrounding the tumor area, as they, relatively, are able to penetrate far away from their areas of production, while the rest of secondary particles (alpha particles, deuterons, electrons, tritons and triple helium) produced deposited most of its energy locally, i.e. close to its production areas, and it contributes to increase the dose introduced to tumor area without damaging the healthy tissue.
Finally, the use of protons in radiation therapy is a very effective method compared to traditional methods (treatment with photons and neutrons) in terms of delivering a large amount of energy to the tumor area compared to the dose provided to the healthy tissue.
Figure 1.Bragg peak for alpha particles in air Figure 2. CT images of the abdomen area showing the liver of a 60-yearold patient

Figure 3 .
Figure 3. Algorithm for building a 3D model of liver based on DICOM images Figure 4. Equivalent voxel model using MCNPX code for each CT slice Figure 5.The phantom of the abdominal area in which the internal organs used in this study are shown depending on reading of CT DICOM images

Figure 7 .
Figure 7.The changes in deposited energy at the Bragg peak (end of the proton path) vs proton's energy for each source proton

Figure 8 .
Figure 8.The changes of the total energy deposited over the entire proton path for each source proton

Figure 10 .
Figure 10.The flux of secondary neutrons as function of the depth within the tissue for each source proton

Figure 11 .
Figure 11.The secondary photons flux as a function of the depth within the tissue for each source proton Figure 13.The percentage changes of secondary particles flux as a function of the energy of the incoming proton