# Monte Carlo Simulation of CRAND Protons Trapped at Low Earth Orbits

Ritabrata Sarkar<EMAIL_ADDRESS>Abhijit Roy<EMAIL_ADDRESS>Indian Centre for Space Physics, 43 Chalantika, Garia Station Road, Kolkata
700084, W.B., India

###### Abstract

The Cosmic Ray Albedo Neutron Decay (CRAND) is believed to be the principal
mechanism for the formation of inner proton radiation belt – at least for
relatively higher energy particles. We implement this mechanism in a Monte
Carlo simulation procedure to calculate the trapped proton radiation at the
low Earth orbits, through event-by-event interaction of the cosmic ray
particles in the Earth’s atmosphere and their transportation in the
magnetosphere. We consider the generation of protons from subsequent decay of
the secondary neutrons from the cosmic ray interaction in the atmosphere and
their transport (and/or trapping) in the geomagnetic field. We address the
computational challenges for this type of calculations and develop an
optimized algorithm to minimize the computation time. We consider a full 3D
description of the Earth’s atmospheric and magnetic-field configurations using
the latest available models. We present the spatial and phase-space
distribution of the trapped protons considering the adiabatic invariants and
other parameters at the low Earth orbits. We compare the simulation results
with the trapped proton flux measurements made by PAMELA experiment at low
Earth orbit and explain certain features observed by the measurement.

###### keywords:

Trapped protons at low Earth orbit , CRAND protons , Monte Carlo simulation ,
Radiation belt

††journal: Advances in Space Research

##  1 Introduction

The study of the inner radiation belt is important, as this hostile radiation
environment impose crucial hazards to the satellites and space missions at low
and middle Earth orbits. While, the outer radiation belt is predominantly
composed of energetic electrons, the inner belt mainly consists of energetic
protons up to a few GeV, with a trace of e± and some other ions. The trapped
protons in the radiation belt may have long residence time (up to several
thousand years) and thus the accumulated radiation effect results from long
time inputs from different sources like solar protons and Galactic Cosmic Rays
(GCR). Several processes control the internal dynamics of the radiation belts.
Radial diffusion provides the inward transport of the solar protons injected
into the radiation belt at L $\gtrsim$ 2, during large solar proton events and
magnetic storms [12, 9]. Losses of radiation belt particles are imposed by the
processes like: ionization of the neutral atmosphere, energy loss to the
plasmasphere and ionosphere, inelastic nuclear scattering and adiabatic energy
and drift shell changes due to geomagnetic secular variations. For the
radiation intensity at energies $\lesssim$ 100 MeV and for L $\gtrsim$ 1.3 are
dominated by the solar protons [20]. Otherwise, the protons from the
$\beta$-decay of the free neutrons produced by the GCR interaction with the
Earth’s atmosphere, according to Cosmic Ray Albedo Neutron Decay (CRAND)
mechanism [21, 4], dominates the radiation belt.

There are some empirical trapped proton models available, to represent the
inner radiation belt, like: NASA AP8/AP9 model [18, 6], the model provided by
the Institute of Nuclear Physics of Moscow State University [5], PSB97 model
[8] etc. Selesnick et al. [20] came up with a theoretical model of the inner
proton radiation belt considering the physics and mechanism that populates the
radiation belt. There are several satellite experiments that contributed to
the understanding and modeling of the radiation belt like CRRES [7],
SAMPEX/PET [11], TIROS/NOAA series [10] etc. Lately, measurement of CR
radiation at Low Earth Orbit (LEO) by PAMELA mission [15], provides important
information on the high energy ($\gtrsim$ 70 MeV) part of the trapped
radiation at the lower region of the radiation belt.

In this work, we recreate the trapped proton distribution at LEO region,
considering realistic physical models of: Earth’s atmosphere, magnetosphere,
GCR distribution, particle interactions in the atmosphere and their transport.
We embedded all the physical mechanisms in a single Monte Carlo (MC)
simulation framework. However, the total simulation procedure is segregated
into different stages for convenience of the calculation and the outcome of
each stage of the simulation flows into the following to get the final
results. The details of the simulation procedure are discussed in Sec. 2. The
spatial and spectral distribution of the secondary particles produced from the
GCR interactions with the atmospheric nuclei depend on the rigidity cut-off of
the GCR at different locations determined by the geomagnetic field
distribution. The albedo neutrons produced in this way will have location
dependency. The protons generated from the subsequent $\beta$-decay get
trapped or escape according to their rigidity that depends on the energy and
location of their generation. To include all these effects we perform the
simulation starting right from the generation of GCR distribution at the
Earth’s orbit, their transportation into the Earth’s atmosphere depending on
their rigidity and subsequent interaction to produce secondary particles.
Then, we consider the neutrons which escape into the space and decay to
produce protons. The trajectories of these protons can be calculated
individually from the equation of motion in the magnetic field. However,
calculating the actual path due to Newton-Lorentz motion may lead to time
consuming computation. So, instead, we adopted the guiding center equation to
calculate the average track of the particles. In this work, we carry on with
some simplified assumptions, like neglecting geomagnetic secular variation and
considering an average heliospheric effect which, in principle, can limit some
aspects of the calculation and will be discussed later at proper places (viz.
Sec. 3.4). We concentrate on the high energy part $\gtrsim$ 50 MeV to compare
the simulation results with PAMELA measurement of trapped protons [1].

In the following Sec. 2, we outline the general simulation procedure to
describe the neutron generation from the GCR interaction with Earth’s
atmosphere, proton production from subsequent decay and their entrapment in
the geomagnetic field. In Sec. 3, we present the analysis procedure and
results from the calculations. We compare the simulation outcome with the
recent measurements of the trapped protons by the PAMELA mission in Sec. 4.
Finally we conclude the results from this work in Sec. 5.

##  2 Simulation Procedure

We considered an event-by-event simulation procedure using the MC technique to
study the problem of trapped protons from the CRAND mechanism. The general
scheme for the simulation is as follows. At first we calculate the neutron
production in the atmosphere due to the interaction of the CR particles. These
free neutrons, therefore decay to produce protons according to the
$\beta$-decay process. Finally, the protons either get trapped into the
radiation belt, transported off to the outer space or absorbed at the Earth –
according to their position of generation, rigidity, the direction of motion
(pitch angle) and the geomagnetic field distribution. In the following
subsections we briefly discuss these steps of the simulation procedure.

###  2.1 Neutron production in the atmosphere

To calculate the neutron production from the CR interaction in the atmosphere,
we considered a detailed MC simulation using Geant4 simulation toolkit [2].
The details of this simulation procedure can be found in Sarkar et al. [17].
However, for the sake of completeness we briefly mention the key aspects of
the simulation here. We implemented a full 3D model of the atmosphere and
magnetosphere in Geant4 simulation environment, with proper distribution of
primary GCR particles at the Earth’s orbit and using a list of interaction
processes suitable for the energies and phenomena considered in the
simulation.

The atmospheric matter distribution which depends on the location, time and
other conditions like solar activity, is realized using the NRLMSISE-00
standard atmospheric model [14]. The atmospheric distribution is limited up to
100 km from the Earth’s surface since most of the CR interactions effectively
take place below this level to produce neutrons, which is our prime concern
for this step of the simulation. The whole atmosphere is subdivided into 100
concentric spherical layers with the thickness equal in logarithmic scale of
altitude to implement the gradual material distribution. However, with an
approximation, to evade the complexity in the atmospheric geometry
construction, we did not implement the zonal and meridional distribution of
the atmospheric constituents in the current simulation. We neither considered
the temporal variation of the model parameters. Instead, we considered global
atmospheric parameters at each layer corresponding to the tropical region
(23∘N, 88∘E) and at a fixed time in the month of May in 2016 [17].

The magnetic field distribution is also embedded in the simulation environment
along with the atmospheric matter distribution which directly takes care of
the rigidity cut-off of the charged primary particles and defines the
trajectories of the primary and secondary charged particles. The inner
magnetic field due to the Earth’s inherent magnetism extends up to about 4
Earth radii ($R_{E}$) from the surface and is represented by the 12th
generation IGRF model [24] with proper input parameters at moderate solar
condition (during 2016 at declining phase of 24th solar cycle; vide Sarkar et
al. [17]). The external magnetic field depends on the inter-planetary magnetic
field and solar conditions, is considered up to 25 $R_{E}$ defined by using
the Tsyganenko model [25]. Here also, in defining the magnetic field model, we
did not consider the temporal variation due to solar activity or the secular
variation of the geomagnetic field, instead used the average field
distribution corresponding to the above mentioned time.

For the primary GCR flux distribution, we considered only the most abundant
components, i.e., protons (H) and helium nuclei (He) which constitute about
99% of the GCR flux. We used up-to-date spectral information of the primary CR
particle flux using the model for local interstellar spectrum given by Vos &
Potgieter [26]. This model is conceived by using the Voyager I data of the
post-heliopause CR observation complimented by PAMELA and AMS02 data. The
details of the primary CR flux distribution is given in Sarkar et al. [17]. To
maintain the efficiency of the simulation, we generated an isotropic
distribution of the charged particles from a geocentric spherical surface at
500 km above the Earth’s surface and backtracked them to reach an outer
surface at 25 $R_{E}$ in presence of the magnetic field distribution. We
select only those particles that reach the outer sphere as the allowed CR
particle tracks at 500 km and proceed with the simulation for atmospheric
interactions. Thus, we consider the modulation of the primary CR at the
Earth’s orbit by the geomagnetic field surrounding Earth till 25 $R_{E}$. We
restricted the primary generation in the energy range of 0.1–800.0 GeV/n. We
simulated 3.5 $\times$ 105 particles (both for H and He) that survived the
rigidity cut-off at 500 km above the Earth’s surface in the back-tracking
method, to interact with the atmosphere. To achieve this number of particles
to interact in the atmosphere, the simulation discarded about 106 H and 1.5
$\times$ 106 He particles. Simulation of these number of events was required
to validate simulation results with the AMS02 observations with comparable
significance [17].

To simulate the interaction of high energy particles with the atmosphere we
used the reference physics list (QGSP physics list with Binary Cascade model
QGSP_BIC_HP) provided in Geant4, which optimally covers the interaction
processes relevant here.

Here, we are particularly interested in the neutrons produced due to CR
interactions in the atmosphere, which move outwards from the atmosphere. So,
we calculated the secondary neutron distribution at the spherical surface at
100 km above the Earth’s surface, which are outgoing from the atmosphere into
the space (albedo neutrons). The spatial distribution of the secondary
neutrons depend on the rigidity cut-off of the primary CR particles and hence
on the geomagnetic field distribution. The spatial distribution (in geographic
latitude ($\theta$)/longitude ($\phi$)) of the neutrons at 100 km above the
Earth’s surface is shown in Fig. 1 along with vertical rigidity cut-off values
at different locations. The apparent drop-off in the neutron flux toward the
poles are due to the coordinate ($\theta$-$\phi$) we chose for the plotting.
(For a uniform positional distribution on a spherical surface $\theta$ is not
isotropic in [0, $\pi$] rather $\cos\theta$ is uniform in [0, 1].) The
overlaid contour plot is also calculated according to the chosen coordinates
and this note on the distortion of the plots in $\theta$-$\phi$ coordinate is
also valid for other subsequent plots (viz. Fig. 9 and 10). The positional
distribution, as well as the directional and spectral distributions of the
secondary neutrons are used to generate the CRAND protons described in the
next Sec. 2.2. We used the calculated positional distribution shown in Fig. 1
to generate random neutron positions in $\theta$/$\phi$ for the next step of
the simulation.

Figure 1: Positional distribution of the secondary neutron flux at the top of
the atmosphere (100 km from Earth’s surface) as obtained from the simulation,
due to primary cosmic rays (H and He). The contour plot of vertical rigidity
cut-off values at different locations along with the continental coast line
map for physical location guidance are also superimposed in the plot.

The propagation (or momentum) direction of the produced neutrons are also
important in this concern; since the pitch angle of the protons, generated
from the decay of these neutrons, depends on their direction. We have shown
this directional dependence in terms of $\theta_{Z}/\phi_{Z}$ angle in Fig. 2.
Where $\theta_{Z}$ is the angle between propagation direction and zenith
direction at the position of the neutron and $\phi_{Z}$ is the azimuthal angle
about the zenith direction and with respect to local north. The cut-off that
appeared in the $\cos\theta_{Z}$ distribution of the albedo neutron is due to
the fact that, we generated the primary H and He from a spherical surface at
500 km above the Earth’s surface while the extent of the atmosphere is only up
to 100 km. So, only the primaries inside a limited solid angle (about the
zenith direction) among all those isotropically generated from the spherical
surface, can interact with the atmosphere. While this particular configuration
of the primary generation was originally considered to serve some general
purpose to calculate space and atmospheric radiation due to cosmic rays as
discussed in Sarkar et al. [17], this may have some minor effect on the low
pitch angle trapped particle population as will be discussed later at proper
place (viz. Sec. 3.1 and 3.2).

(a)  (b)

Figure 2: (a) Neutron directional distribution as a function of
$\cos\theta_{Z}$. (b) Neutron directional distribution as a function of the
azimuthal angle ($\phi_{Z}$).

The energy distribution of the produced neutron flux is shown in Fig. 3. This
spectrum is fitted with an empirical function with three modified cut-off
power-laws (at three different energy intervals) of the form:

$\frac{dN}{dE}=\sum\limits_{i=1}^{3}A_{i}E_{i}^{-\alpha_{i}}exp\left(-\left(\frac{E_{i}}{B_{i}}\right)^{\beta_{i}}\right),$
(1)

where the energy $E_{i}$ is in the ranges $(\leftarrow:1.5\times
10^{-2}:2.0\times 10^{-2}:\rightarrow)$ GeV. The three sets of ($A_{i}$,
$\alpha_{i}$, $B_{i}$, $\beta_{i}$) values obtained from the fit are (4.06712
$\times$ 10+4, 3.03 $\times$ 10-1, 2.27 $\times$ 10-3, 6.14 $\times$ 10-1);
(7.77445 $\times$ 10+4, -3.07 $\times$ 10-1, 1.30 $\times$ 10-2, 5.67 $\times$
10-1); (1.45930 $\times$ 10+3, 1.75, 4.84 $\times$ 10-5, 1.68 $\times$ 10-1)
respectively for three energy ranges. We used this function to generate the
neutrons for subsequent generation of CRAND protons described in the next Sec.
2.2.

Figure 3: Albedo neutron flux distribution at the top of the atmosphere (at
100 km) due to cosmic ray (H and He) interactions in the atmosphere. The
spectrum is fitted (dashed line) with an empirical function given by Eq. 1.

It might be worth to mention here that actually the energy and directional
distribution of the secondary neutrons may depend on the location (with
geomagnetic latitude through rigidity cut-off and also to some extent due to
atmospheric distribution). So, to achieve a more realistic neutron generation
model in this phase of the simulation, one might incorporate the location
dependent energy and directional distribution. In the current simulation, for
the sake of simplicity, we adopted the averaged distribution in energy and
direction over the whole area. However, we confirmed that at different
geomagnetic latitude the shape of the energy spectra remains more or less
identical. The variation in amplitude of the omnidirectional energy spectrum
is taken care by the location dependent particle distribution shown in Fig. 1.
We also observed only a nominal deviation in $\theta_{Z}$ distribution at the
equatorial region (0–30∘ geomagnetic latitude) with respect to the average
plot shown in Fig. 2a. We estimated about 20% lower relative flux at the
higher $\theta_{Z}$ (or lower $\cos\theta_{Z}$, i.e., away from zenith) end in
this comparison.

###  2.2 Protons from albedo neutron decay

To generate the protons from the CRAND source, we used the spatial,
directional and spectral distribution of albedo neutrons shown in Fig. 1, 2
and 3. At first we generated random neutrons according to these distributions
on a spherical surface at 100 km above Earth’s surface. Then protons were
generated according to free neutron $\beta$-decay scheme. We randomly sampled
the neutron decay time ($t_{D}$) for each neutron from the decay time
distribution function:
$\frac{t_{D}}{\tau_{n}}\,exp\left(\frac{-t_{D}}{\tau_{n}}\right)$, where,
$\tau_{n}$ = 881.5 s is the neutron mean life time. Then we calculated the
position of generation of the proton from the distance travelled by the
neutron before its decay as: $v_{n}\,\gamma\,t_{D}$, where, $v_{n}$ is neutron
velocity (calculated from the kinetic energy and track direction) and $\gamma$
is the relativistic time dilation factor. For the sake of simplicity, we
preserved the direction and kinetic energy of the protons same as the
neutrons. In the original simulation calculation we considered all the protons
generated inside the magneto-pause (10 $R_{E}$) boundary. However, considering
our current purpose of calculating trapped protons and comparing with PAMELA
measurement [1], we enforced two conditions limited by the measurement – take
account of only those protons which: (i) have kinetic energies in the range of
50 MeV to 4 GeV and (ii) are generated below 2 $R_{E}$. The first limit is
obvious due to the energy threshold of the instrument and maximum energy of
the trapped particles found by the measurement. The second limit imposed here
considering the fact that PAMELA can observe the stably trapped protons only
for L-shell values in between $\sim$1.18 $R_{E}$ up to $\sim$1.7 $R_{E}$ [1].

###  2.3 Proton transportation in the geomagnetic field

The trajectory of the protons produced from the decay of neutrons are guided
by the geomagnetic field according to their momentum and position. These
particles will then be trapped in the magnetic field or deflected (into the
outer space or towards Earth) depending on the mentioned parameters and their
tracks can be traced using numerical calculations [22, 23].

However, the calculation of the actual charged particle tracks due to its
Newton-Lorentz motion in the magnetic field can be computationally expensive.
So, instead, we considered the guiding center motion of the particles, while
the directional information of the original track is preserved by the pitch
angle information. An offhand calculation (with 1000 proton tracks of 50 MeV
energy) shows a reduction of the tracking time to $\sim$ 37%.

The equations of motion of the guiding center are given by [13]:

$\frac{d\mathbf{R}}{dt}=\frac{\gamma
mv^{2}}{2qB^{2}}\left(1+\frac{{v_{\parallel}}^{2}}{v^{2}}\right)\hat{\mathbf{b}}\times{\mathbf{\nabla}}B+v_{\parallel}\hat{\mathbf{b}},$
(2)

$\frac{dv_{\parallel}}{dt}=-\frac{\mu}{\gamma^{2}m}\hat{\mathbf{b}}\cdot{\mathbf{\nabla}}B.$
(3)

Here $\mathbf{R}$ is the guiding center position of the particle track, which
moves in the local magnetic field with magnitude $B$ and directional unit
vector $\hat{\mathbf{b}}$. $m$, $q$, $\gamma$, $v$ and $v_{\parallel}$ are the
particle (proton) mass, charge, Lorentz factor, speed and component of
velocity parallel to the magnetic field direction respectively.
$\mu=\frac{\gamma^{2}m{v_{\perp}}^{2}}{2B}$ (with $v_{\perp}$ the
perpendicular component of the particle velocity), is the magnetic moment of
the current generated by the circular motion of the charged particle in the
magnetic field.

We calculated the guiding center trajectories of the particle motion in the
geomagnetic field by integrating the Equations 2 and 3. We considered the
adaptive integration technique using boost::numeric::odeint library [3] with
Cash-Karp algorithm (runge_kutta_cash_karp54) for stepper. The initial values
supplied to the integration method are: the position of the generated proton
from the neutron decay, the velocity component of the particle parallel to the
magnetic field and initial stepping time depending on the particle velocity to
travel 500 km distance along the guiding center.

However, we record only the stably-trapped particle trajectories satisfying
the conditions that, while drifting through the full rotation around the Earth
($360^{\circ}$ change of azimuthal angle of the trajectory; we are using
geocentric coordinate system) it does not reach: (i) within 40 km from Earth’s
surface (absorbing atmospheric limit) and (ii) the magnetopause limit at 10
$R_{E}$. Both of these two conditions corresponds to the quasi-trapped or un-
trapped particles and are not considered in the present work.

The termination of the stably-trapped trajectory integration depends on the
morphology of the drift trajectory. Keeping in mind, the latitude/longitude
grid resolution of the order of a degree, during the spatial analysis of the
trajectory data, we intend to keep an azimuthal distance of $\sim 1^{\circ}$
between two successive crossing points of the trajectory with the latitude
(say at magnetic-equator). Then, we calculate the number of full drifted
rotation of the trajectory required to cross the equator at least 360 times.
So the trajectories with higher drift speed (with larger separation between
successive crossings) are traced for longer time. This is done for better
interpolation of the particle location on the drift shell. The overall
information of each trajectory is saved for the analysis purpose described in
the next section.

We simulated a total of 106 events in which protons generated from the neutron
decay are trapped in the geomagnetic field, satisfying the aforementioned
trapping conditions. To get this number of trapped particles, the simulation
discarded another 2793176 events which did not satisfied the the trapping
conditions and correspond to quasi-trapped or untrapped particles and not been
considered further.

##  3 Data Analysis and Results

###  3.1 Physical parameters

At first, we check the three most important parameters at the generation point
of the proton trajectories: position of generation, direction of propagation
and kinetic energy of the particle — to look into the facts of delay in
neutron decay, coupling of the proton with the geomagnetic field and
energetics of the particle distribution. The initial position distribution of
the trapped protons is shown in Fig. 4a, which shows the distance of proton
generation from the Earth’s center (in this simulation we considered this
distance only up to 2 $R_{E}$). It is evident from this figure that the number
of CRAND generated protons are increasing beyond 2 $R_{E}$, but, as we have
already mentioned earlier, this limit is imposed in the prospect of comparing
the results with PAMELA measurement. We show the overall distribution of the
initial pitch angle ($\alpha_{i}$; i.e., at the position of the proton
generation) and the equatorial pitch angle ($\alpha_{0}$) of the trapped
protons in Fig. 4b. The $\alpha_{0}$ values are calculated directly from the
magnetic field information of the trajectory at the mirror points ($B_{m}$)
and magnetic equator ($B_{0}$) as $\alpha_{0}$ =
$\sin^{-1}\sqrt{\frac{B_{0}}{B_{m}}}$, averaged over the whole drift
trajectory. The higher population near lower pitch angle (i.e., higher
probability of particles having velocity nearly along or opposite to the
magnetic field lines) can be explained from the following facts. The neutrons
are generated in majority near the polar regions (where the inclination of the
field lines are larger), as evident from Fig. 1. In addition, more particles
are directed towards the zenith as suggested by the distribution of the
momentum directions (Fig. 2). However, the scarcity of the protons with
extremely low values of initial pitch angle (parallel or anti-parallel to the
magnetic field) is due to the effects both from position distribution of the
neutron generation and simulation artifact arose from primary particle
generation surface and atmospheric distribution consideration as mentioned in
Sec. 2.1.

(a)  (b)

Figure 4: (a) Trapped proton distribution as a function of the distance of
initial position of the trajectories. (b) Trapped proton distribution as the
function of cosine of the initial pitch angle (black) and the equatorial pitch
angle (gray).

In Fig. 5a, we show the trapped proton distribution with respect to the
kinetic energy and the distance of its generation. The overall energy spectrum
is shown In Fig. 5b. While the total entrapment efficiency of the CRAND
protons in the considered energy range and altitude limit can be estimated
from total number of simulated neutrons and entrapped protons given in Sec.
2.3 which is about 26%, a relative calculation of the energy dependent
entrapment efficiency was also done. This shows about 40% lesser trapping
efficiency at the lower energy end in comparison to that at the higher end and
the efficiency increases almost monotonically towards the higher energy end.

(a)  (b)

Figure 5: (a) Trapped proton distribution as the function of kinetic energy
and distance of initial position of the trajectory. (b) Integrated flux
distribution of the trapped protons.

We verified the trapped proton trajectories to fulfill the adiabatic
conditions, particularly checking the hierarchy of the temporal scales:

$\tau_{gyro}\ll\tau_{bounce}\ll\tau_{drift},$ (4)

where $\tau_{gyro}$, $\tau_{bounce}$ and $\tau_{drift}$ are time scales
associated to gyration, bounce and drift motions, respectively. $\tau_{gyro}$
($=\frac{2\pi m_{p}}{qB}$; where $m_{p}$ proton mass, $q$ charge, $B$ local
magnetic field) is calculated from the magnetic field value at the initial
position of the proton trajectory. While, $\tau_{bounce}$ and $\tau_{drift}$
are calculated by averaging the bounce periods and drift periods respectively,
over the whole trajectory, obtained from the simulated trajectory information.
Figure 6 shows the $\tau_{gyro}$, $\tau_{bounce}$ and $\tau_{drift}$
distribution for all the simulated trajectories.

Figure 6: Distribution of temporal scales for three types of motions of all
the simulated trapped protons under the effect of geomagnetic field.

###  3.2 Adiabatic invariants

Here, in this work, we have not considered the gradual energy loss, nuclear
scattering or radial diffusion of the trapped particles. So the adiabatic
invariants of the trapped particles are conserved. The versions of the
adiabatic invariants considered are as follows [16]:

$M=\frac{\gamma^{2}mv_{\perp}^{2}}{2B},$ (5)

$K=\int_{s_{m}}^{s^{\prime}_{m}}\left(B_{m}-B(s)\right)^{1/2}ds,$ (6)

$L=\frac{2\pi\mu_{E}}{R_{E}\Phi}.$ (7)

The first invariant $M$ is calculated using the particle velocity and magnetic
field value at the initial position of the track. The second invariant $K$ is
calculated by integrating the magnetic field along the guiding center path
between two successive mirror points. The average $K$ value over the whole
trajectory has been used as the second invariant of the trajectory. We
consider $L$ value as the third invariant. It is directly calculated from
distance of the trajectory from the Earth’s center at the geomagnetic equator,
where the track is at the highest position. We subsequently calculate $\Phi$,
the magnetic flux inside the drift shell from Eq. 7, where $R_{E}$ is Earth’s
radius and $\mu_{E}$ Earth’s magnetic dipole moment (arbitrarily taken for
2008 A.D. to compare with the results from PAMELA measurements).

The distribution of the three adiabatic invariants, for all the simulated
particles are shown in Fig. 7. The drop in population at extreme low $M$
values may be attributed to the scarcity of trapped particles at the lowest
values of initial pitch angle as shown in Fig. 4b. The cut-offs apparent in
the $K$ plot at around $10^{-1}$ $T^{1/2}R_{E}$ and in $L$ plot at 2 $R_{E}$
are due to the limit used in the simulation for the particle generation
distance. However, we calculated the trapped particle trajectories beyond this
distance (up to the average magnetopause distance at 10 $R_{E}$).

(a)  (b)  (c)

Figure 7: Distribution of adiabatic invariants (a) M, (b) K and (c) L values
for all the simulated trajectories.

###  3.3 Spatial distribution of trapped protons

A stably trapped particle track creates a drift shell surrounding the Earth
and can be detected from any position in the continuous shell given a
sufficiently long entrapment time. So, from the simulation point-of-view, we
can spatially distribute one event to contribute all over the positions
throughout the drift shell, improving the particle statistics without using
extra computational resource for simulating more events. However, we must
interpolate the trajectory to represent the drift shell. In Fig. 8a, we show,
for example, the guiding center trajectory of a 288.7 MeV proton trapped in
the geomagnetic field started from a position at 1.66 $R_{E}$ altitude and
with initial pitch angle 135.9∘. The trajectory chosen mere arbitrarily just
to present the analysis procedure. For the analysis purpose, we replotted the
trajectory in the geocentric coordinate system ($\phi$, $\theta$ and elevation
from Earth’s center), with $1^{\circ}\times 1^{\circ}$ grid in $\phi$/$\theta$
and dividing the elevation range of 1–2 $R_{E}$ in 320 bins (giving $\Delta
Alt$ $\sim$ 20 km). The trajectory plot is shown in Fig. 8b which represents
the drift shell for the corresponding track. Then we can calculate the
distribution of the drift shell interceptions at each layer of different
elevation (measurement altitude) and interpolate the interception locations.
This is shown in Fig. 8c, for example at the measurement surface at 1.7
$R_{E}$. We adopted this method to ensure that, for the simulated
trajectories, where the azimuthal gap between two subsequent bouncing tracks
are larger (than the considered phi grid, i.e., 1∘ here), we don’t have to
rotate the trajectory around Earth for too many times, which can be
computationally expensive. Instead, we fill the gaps (between subsequent
locations of track interceptions) by interpolation as shown in Fig. 8c.

(a)

(b)

(c)

Figure 8: (a) Example of a trapped particle track (guiding center) for a 288.7
MeV proton started at 1.66 $R_{E}$ with initial pitch angle 135.9∘. (b) Same
track represented in the grid view of geocentric coordinate system. (c)
Distribution of the drift shell interception locations with a spherical
surface at 1.7 $R_{E}$.

Then we calculate the overall contribution of all the simulated trapped
trajectories at different measurement altitudes. In Fig. 9, for example, we
show the spatial distribution of proton trajectories integrated over four
different altitude ranges at LEO: 300–385, 385–470, 470–555 and 555–640 km. It
is evident from these plots that the high energy trapped protons in the
considered energy range and with all pitch angles are largely concentrated at
north-east of the geographic south pole below the African peninsula.

Figure 9: Distribution of trapped proton trajectories from CRAND in
geographical latitude/longitude at four different altitude ranges. Geomagnetic
field contours at the corresponding heights has been drawn overlaid on the
proton distribution plots.

###  3.4 Flux normalization

The trapped proton intensity at some fixed location (on a satellite orbit) and
at a given time (of measurement) can be expressed as a function of kinetic
energy ($E$) and time, in the following form [19]:

$j(E,t)=v\int_{t_{0}}^{t}S(E^{\prime},t^{\prime})\exp\left(-\int_{t^{\prime}}^{t}\left(\frac{\partial}{\partial
E^{\prime\prime}}\frac{dE}{dt}+\frac{1}{\tau}\right)dt^{\prime}\right),$ (8)

where $v$ is the speed of the particle corresponding to the kinetic energy,
$S$ is the source rate, $\tau$ is the trapping lifetime and $\frac{dE}{dt}$
represents the continuous energy change of the particle. At the initial time
$t$ = $t_{0}$, $j$ = 0, when $E$ corresponds to maximum energy of the trapped
proton. However, in this work, we consider crude approximation of no gradual
energy loss (through ionization; to free electrons or adiabatic energy change
caused by secular geomagnetic variation) and no trapped particle loss through
nuclear scattering or due to the effect from geomagnetic storms. So,
$\frac{dE}{dt}$ = 0 and $\tau$ = $\infty$; we only have the time integration
over the source term, which, in this present case, is the GCR contribution
through CRAND mechanism. This approximation has severe implication over the
absolute flux calculation of the trapped particles but reveals the relative
flux distribution anyway.

The normalization factor of the simulation events to calculate the trapped
proton flux can be given by:

$NF=\frac{\int_{T_{s}}\int_{A_{s}}\int_{\Omega_{s}}\int_{E_{s}}\frac{dN}{dE}dEd\Omega
dadt}{N_{sim}\int_{A_{m}}\int_{\Omega_{m}}d\Omega da},$ (9)

where $\frac{dN}{dE}$ is the source function of the albedo neutron energy
spectrum given by Eq. 1 and $N_{sim}$ (= 106) is the total number of simulated
events giving trapped particles. Integration of the source function is done
over the energy range $E_{s}$ = [50 MeV, 4 GeV], over 2$\pi$ solid angle
($\Omega_{s}$) and over the area ($A_{s}$) of albedo neutron generation
surface (spherical surface at 100 km from Earth’s surface). Integration for
the part at the denominator, corresponds to the measurement of the trapped
protons and is carried over 4$\pi$ solid angle ($\Omega_{m}$) and over the
spherical surface area $A_{m}$ at the measurement altitude ($\sim$ 500 km)
considering the average PAMELA orbit. The time limit $T_{s}$ , which indicates
the accumulation time of the trapped particles, is undefined here. Since, we
are not considering the loss of the trapped particles, the average value of
$T_{s}$ can only be inferred using the information of the measured flux. In
this calculation, considering a gross comparison with the PAMELA measurement
of the trapped proton flux, we get $T_{s}$ $\sim$ 0.1 year. This value is
somewhat lower in comparison to the trapped proton residence time range
mentioned in Selesnick et al. [20] as less than a year to $\sim$ 4000 years.
But this can be justified for the high energy particles (with smaller drift
time scale) we are considering here and since we are ignoring the particle
loss effects.

##  4 Comparison with PAMELA Measurements

We compare the simulation results with the measurement of trapped proton
fluxes at LEO by the PAMELA experiment [1]. PAMELA measured the trapped proton
fluxes at $350\div 610$ km altitude between July, 2006 and September 2009. The
proton trajectories were reconstructed using tracing programs based on
numerical integration method [22, 23] to select the stably trapped
trajectories used for the analysis. A total of 7.3 $\times$ 106 stably trapped
particles were sampled for the analysis out of total 9 $\times$ 106 detected
events (including quasi- and un-trapped trajectories). We considered similar
order of event number in the current simulation to get comparable statistics
in the results.

In Fig. 10, we show the distribution of trapped proton flux integrated over
the pitch-angle range covered by PAMELA, at different altitudes and energy
ranges. We refer to the Fig. 1 of Adriani et al. [1] for the comparison of the
distribution. We rebinned the $\theta/\phi$ distribution in $2^{\circ}\times
2^{\circ}$ grid for the comparison. We also added an uncertainty of
$\sigma=3^{\circ}$ to the simulated track positions (in $\theta/\phi$), to
account the uncertainty in payload orientation measurement and angular
resolution of the detector. We sampled the detector location over all the
region at different altitude (but limited by the satellite orbit inclination
of 70∘) and obtained the trajectories within the acceptance angle of the
detector, by calculating the angle between the detector direction (towards the
local zenith) and trajectory direction (obtained from the local pitch angle of
the trajectory) at the location of the detector. However, this approach
averages the calculation over the gyro-phase angle and neglects the east-west
effect due to atmospheric absorption at larger gyro-radius. The simulated
result of the trapped proton distribution map precisely reproduces the
measured distribution by PAMELA which is concentrated around the South
Atlantic Anomaly (SAA) region. Relation of the area coverage with the altitude
and energy range is also apparent from the plot. In particular, the feature
reported in Adriani et al. [1], that the proton flux is more concentrated in
the southeast part of SAA is also apparent here. To make this feature more
apparent we also indicate the peak points of the distribution on the same
plot. The explanation of the peaks shifted towards the southeast can be
visible from Fig. 9, which shows that the trapped protons, irrespective of
PAMELA acceptance angle, are concentrated near the south-polar region below
the African peninsula. We have already mentioned in Sec. 3.4, the absolute
value of the integral flux depends on the accumulation time of the trapped
particles, which require to include the loss effects to calculate it. Since,
in this work we did not consider the loss effects, to present the integral
flux we simply scaled the value by varying the $T_{s}$ value to grossly match
the PAMELA measurement.

Figure 10: Simulated trapped proton integral flux distribution, integrated
over the pitch angle range covered by PAMELA, in geographical coordinates, at
different altitudes (rows) and energy ranges (columns). The peak of the
distributions are indicated by the plus signs.

For additional comparison with the measured data, we also plot the integral
proton flux with other parameters than spatial distribution. Figure 11a shows
the proton integral fluxes as a function of second and third adiabatic
invariants ($K$ and $\Phi$), for the particles in PAMELA field of view, at
different kinetic energy ranges. We plotted the result in $K^{1/2}$ and
$\log_{10}\Phi$ coordinates for better resolution (and for the sake of direct
comparison with PAMELA results). Figure 11b shows the same fluxes as a
function of equatorial pitch angle and $L$ value. These results can be
compared with the results from the PAMELA measurements shown in Fig. 2 and 3
of Adriani et al. [1]. While the overall distribution agrees well with the
observed data, the comparisons show some differences in the extent of the
energy dependent data patches in $K$ and $\alpha_{0}$. These differences can
be attributed, as discussed in Adriani et al. [1], to the spacecraft orbit
which constraint the observation of the equatorial mirroring protons only for
L values up to $\sim$ 1.18 RE. However, one thing should be mentioned
regarding the simulation aspect — here in this work we optimized the algorithm
to gain statistically in the spatial distribution of the trapped particles
where one trajectory can contribute to the entire drift shell region. But,
since this optimization has no effect on increasing the number of
trajectories, the analysis concerning the parameters attributed to a
trajectory has little gain from the algorithm. We need to simulate more
particles for a better comparison of these aspects shown in Fig. 11.

(a)  (b)

Figure 11: (a) Proton integral fluxes as a function of K and $\Phi$ at
different kinetic energy bins in the pitch angle range covered by PAMELA. (b)
The similar distribution as a function of L-shell parameter and equatorial
pitch angle.

##  5 Conclusions

We conducted a thorough simulation to calculate the trapped proton radiation
contribution at LEO, through the CRAND mechanism. This is the principal
contributor to the proton radiation at this region for proton energy $\gtrsim$
50 MeV. For this purpose we developed a simulation framework to calculate the
secondary particles from the atmospheric interaction of the GCRs. The protons,
produced from the $\beta$-decay of the albedo neutrons from these GCR
interactions, are then geomagnetically trapped and their trajectories have
been calculated using the guiding center equation of the particles. We studied
the spatial as well as the phase-space distribution of the trapped protons and
compared the results with the observation by PAMELA mission, giving
satisfactory agreement. This study also helps to understand certain features
of the measured flux distribution. However, in this calculation we considered
only the source term, neglecting the loss of trapped particles over time, due
to various scattering and energy loss processes. Consideration of the loss
term is important for the complete understanding of the radiation distribution
and will be added in the future extension of the work.

While, in this work we have particularly concentrated on the trapped proton
distribution at LEO and to compare the results with PAMELA measurement, the
same simulation framework can be used, to study the overall inner radiation
belt by investigating the spatial distribution, trapping limits etc. We can
also consider the contributions from other sources to the radiation belt like
solar proton events with some modification to the source particle generation
and the results can be compared to the recent empirical models like AP9 [6]
and that given by Selesnick et al. [20].

The simulated result agrees well with the observation which is evident from
Fig. 10. This flux distribution is mainly due to the trapped trajectories with
L values less than 1.6. However, a closer inspection of the overall flux
distribution shown in Fig. 9 reveals that there are two different populations
of trajectories that comprise the distribution. The stably trapped particles
with L $<$ 2 are visible near the SAA region. Trajectories in the other part
comprising the horns and deposition near the south pole mainly show higher L
values whose entrapment are theoretically perplexed for such high energy
particles. The origin and nature of these “trapped” trajectories could be
better revealed by a more detailed investigation of the morphology of these
trajectories which has not been done yet in this work and is pending for the
future extension of this work.

##  6 Acknowledgments

We express our sincere gratitude to M. Boezio (INFN, Trieste) and A. Bruno
(NASA/GSFC) for the helpful discussions. This work has been done under partial
financial support from the Science and Engineering Research Board (SERB,
Department of Science and Technology, Government of India) project no.
EMR/2016/003870. We also thank the Higher Education department of West Bengal,
for a Grant-In-Aid which allowed us to carry out the research activities at
ICSP.

## References

  * Adriani et al. [2015] Adriani, O., Barbarino, G. C., Bazilevskaya, G. A. et al. (2015).  Trapped proton fluxes at low Earth orbits measured by the PAMELA experiment.  Astrophys. J. Lett., 799:L4, 7–14. doi:10.1088/2041-8205/799/1/L4.
  * Agostinelli et al. [2003] Agostinelli, S., Allison, J., Amako, K. et al. (2003).  Geant4 - a simulation toolkit.  Nucl. Instrum. Methods Phys. Res., Sect. A, 506/3, 250–303. doi:10.1016/S0168-9002(03)01368-8.
  * Boost C++ libraries [2020] Boost C++ libraries (2020).  https://www.boost.org/doc/libs/1_66_0/libs/numeric/odeint/doc/html/index.html.  Accessed: 10 February 2020.
  * Farley & Walt [1971] Farley, T. A., & Walt, M. (1971).  Source and loss processes of protons of the inner radiation belt.  J. Geophys. Res., 76, 8223–8240. doi:10.1029/JA076i034p08223.
  * Getselev et al. [1991] Getselev, I. V., Gusev, A. N., Darchieva, L. A. et al. (1991).  Model of spatial-energetic distribution of charged particles (protons and electrons) fluxes in the Earth’s radiation belts.  INP MSU Preprint, MGU-91-37/241.
  * Ginet et al. [2013] Ginet, G. P., O’Brien, T. P., Huston, S. L., Johnston, W. R., Guild, T. B., Friedel, R., Lindstrom, C. D., Roth, C. J., Whelan, P., Quinn, R. A., Madden, D., Morley, S., & Su, Y.-J. (2013).  AE9, AP9 and SPM: New Models for Specifying the Trapped Energetic Particle and Space Plasma Environment.  Space Sci. Rev., 179, 579–615. doi:10.1007/s11214-013-9964-y.
  * Gussenhoven [1993] Gussenhoven, E. G., M. S.; Mullen (1993).  Space Radiation Effects Program - An overview.  IEEE Trans. Nucl. Sci., 40/2, 221–227. doi:10.1109/23.212345.
  * Heynderickx et al. [1999] Heynderickx, D., Kruglanski, M., Pierrard, V. et al. (1999).  A low altitude trapped proton model for solar minimum conditions based on SAMPEX/PET data.  IEEE Trans. Nucl. Sci., 46, 1475–1480. doi:10.1109/23.819110.
  * Hudson et al. [2004] Hudson, M. K., Kress, B. T., Mazur, J. E. et al. (2004).  3D modeling of shock-induced trapping of solar energetic particles in the Earth’s magnetosphere.  J. Atmos. Sol. Terr. Phys., 66, 1389–1397. doi:10.1016/j.jastp.2004.03.024.
  * Huston et al. [1996] Huston, S. L., Kuck, G. A., & Pfitzer, K. A. (1996).  Low Altitude Trapped Radiation Model Using TIROS/NOAA Data.  In J. F. Lemaire, D. Heynderickx, & D. N. Baker (Eds.), Radiation Belts: Models and Standards (p. 119).  Washington, DC USA: American Geophysical Union.
  * Looper et al. [1996] Looper, M. D., Blake, J. B., Cummings, J. R. et al. (1996).  SAMPEX observations of energetic hydrogen isotopes in the inner zone.  Radiat. Meas., 26/6, 967–978. doi:10.1016/S1350-4487(96)00098-4.
  * Lorentzen et al. [2002] Lorentzen, K. R., Mazur, J. E., Looper, M. D. et al. (2002).  Multisatellite observations of MeV ion injections during storms.  J. Geophys. Res. A: Space Phys., 107(A9), 1231. doi:10.1029/2001JA000276.
  * Öztürk [2012] Öztürk, M. K. (2012).  Trajectories of charged particles trapped in the Earth’s magnetic field.  Am. J. Phys., 80 (5), 420–428. doi:10.1119/1.3684537.
  * Picone et al. [2002] Picone, J. M., Hedin, A. E., Drob, D. P., & Aikin, A. C. (2002).  NRL-MSISE-00 Empirical Model of the Atmosphere: Statistical Comparisons and Scientific Issues.  J. Geophys. Res. A: Space Phys., 107, 1468. doi:10.1029/2002JA009430.
  * Picozza et al. [2007] Picozza, P., Galper, A. M., Castellini, G. et al. (2007).  PAMELA: A payload for antimatter matter exploration and light-nuclei astrophysics.  Astropart. Phys., 27/4, 296–315. doi:10.1016/j.astropartphys.2006.12.002.
  * Roederer [1970] Roederer, J. G. (1970).  Dynamics of geomagnetically trapped radiation.  Springer-Verlag New York Heidelberg Berlin.  doi:10.1007/978-3-642-49300-3.
  * Sarkar et al. [2020] Sarkar, R., Roy, A., & Chakrabarti, S. K. (2020).  Simulation of cosmic rays in the Earth’s atmosphere and interpretation of observed counts in an X-ray detector at balloon altitude near tropical region.  Adv. Space Res., 65, 189–197. doi:10.1016/j.asr.2019.09.046.
  * Sawyer & Vette [1976] Sawyer, D. M., & Vette, J. I. (1976).  AP-8 trapped proton environment for solar maximum and solar minimum.  NASA Tech. Memo, TM-X-72605.
  * Selesnick [2001] Selesnick, R. S. (2001).  Simulation of the anomalous cosmic ray radiation belt with atmospheric production and decay.  Geophys. Res. Lett., 28, 3417–3420. doi:10.1029/2001GL013383.
  * Selesnick et al. [2007] Selesnick, R. S., Looper, M. D., & Mewaldt, R. A. (2007).  A theoretical model of the inner proton radiation belt.  Space Weather, 5, S04003. doi:10.1029/2006SW000275.
  * Singer [1958] Singer, S. F. (1958).  Trapped albedo theory of the radiation belt.  Phys. Rev. Lett., 1, 181–183. doi:10.1103/PhysRevLett.1.181.
  * Smart & Shea [2000] Smart, D. F., & Shea, M. A. (2000).  Final Report, Grant NAG5-8009.  Center for Space Plasmas and Aeronomic Research, The University of Alabama in Huntsville.
  * Smart & Shea [2005] Smart, D. F., & Shea, M. A. (2005).  A review of geomagnetic cutoff rigidities for earth-orbiting spacecraft.  Adv. Space Res., 36, 2012–2020. doi:10.1016/j.asr.2004.09.015.
  * Thébault et al. [2015] Thébault, E., Finley, C., & Toh, H. (2015).  Special issue “International Geomagnetic Reference Field: the 12th generation”.  Earth Planets and Space, 67, 79–98. doi:10.1186/s40623-015-0313-0.
  * Tsyganenko & Andreeva [2016] Tsyganenko, N. A., & Andreeva, V. A. (2016).  An empirical RBF model of the magnetosphere parameterized by interplanetary and ground-based drivers.  J. Geophys. Res. A: Space Phys., 121, 10786–10802. doi:10.1002/2016JA023217.
  * Vos & Potgieter [2015] Vos, E. E., & Potgieter, M. S. (2015).  New modeling of galactic proton modulation during the minimum of solar cycle 23/24.  Astrophys. J., 815, 119. doi:10.1088/0004-637X/815/2/119.


# Early Diagnosis of Chronic Obstructive Pulmonary Disease from Chest X-Rays
using Transfer Learning and Fusion Strategies

Ryan Wang  Department of Computer Science, National Tsing Hua University  Li-
Ching Chen  Department of Computer Science, National Tsing Hua University
Lama Moukheiber  Institute for Medical Engineering and Science, Massachusetts
Institute of Technology  Mira Moukheiber  Institute for Medical Engineering
and Science, Massachusetts Institute of Technology  Dana Moukheiber  Institute
for Medical Engineering and Science, Massachusetts Institute of Technology
Zach Zaiman  Department of Computer Science, Emory University  Sulaiman
Moukheiber  Tess Litchman  Department of Medicine, Beth Israel Deaconess
Medical Center  Kenneth Seastedt  Department of Medicine, Beth Israel
Deaconess Medical Center  Hari Trivedi  Department of Radiology, Emory
University  Rebecca Steinberg  Department of Medicine, Division of General
Medicine, Emory University School of Medicine  Po-Chih Kuo  Department of
Computer Science, National Tsing Hua University  Judy Gichoya  Department of
Radiology, Emory University  Leo Anthony Celi  Institute for Medical
Engineering and Science, Massachusetts Institute of Technology  Department of
Medicine, Beth Israel Deaconess Medical Center

(30 September 2022)

###### Abstract

Chronic obstructive pulmonary disease (COPD) is one of the most common chronic
illnesses in the world and the third leading cause of mortality worldwide but
is often underdiagnosed. COPD is difficult to diagnose early and often is not
diagnosed until later in the disease course. Spirometry tests are the gold
standard for diagnosing COPD but can be difficult to obtain, especially in
resource-poor settings such as low-income countries. Chest X-rays (CXRs),
however, are readily available and may serve as a screening tool to identify
patients with COPD who should undergo further testing. Currently, no research
applies deep learning (DL) algorithms that use large multi-site data and
multi-modal data to detect early-stage COPD patients, and this is the first
study to address this problem using data fusion and model fusion techniques
while also evaluating fairness across demographic subgroups. We, therefore,
sought to develop DL algorithms to detect COPD using CXRs to create a
screening tool for patients where spirometry may be difficult to access. We
use three CXR datasets in our study, CheXpert to pre-train models, MIMIC-CXR,
and Emory-CXR to develop and validate our models. The CXRs from patients in
the early stage of COPD and not on mechanical ventilation are selected for
model training and validation. We apply transfer learning for training on the
base models and found the Xception base model performs better than other base
models, including DenseNet121, ResNet50V2, and MobileNetV2. We visualize the
Grad-CAM heatmaps of the true positive cases on the Xception base model for
both MIMIC-CXR and Emory-CXR test datasets. We further propose two fusion
schemes, (1) model-level fusion, including bagging and stacking methods using
MIMIC-CXR, and (2) data-level fusion, including multi-site data fusion using
MIMIC-CXR and Emory-CXR, and multi-modal fusion using MIMIC-CXRs and MIMIC-IV
EHR, to improve the overall model performance. Fairness analysis is performed
to evaluate if the fusion schemes have a discrepancy in the performance among
different demographic groups (race-ethnicity, sex, and age). The model-level
and multi-modal fusion schemes perform consistently on the MIMIC-CXR internal
test data, and the multi-site data fusion scheme improves the model
generalizability on the Emory-CXR test data compared to the other fusion
strategies. The results demonstrate that DL models can detect COPD using CXRs
only which can facilitate early screening, especially in low-resource regions
where CXRs are more accessible than spirometry. Further studies on using CXR
or other modalities to predict COPD ought in future work. All codes for the
study are available at this link: https://github.com/Ryan-RE-Wang/COPD-project
and the supplementary can be found at
https://docs.google.com/document/d/1N0VJtglD5DvtwbqjYiH8bdvt_LxNnswt3ZtzDyZpj1A/edit?usp=sharing.

##  1 Introduction

Chronic obstructive pulmonary disease (COPD) remains one of the leading causes
of mortality and morbidity [25]. Early intervention can mitigate disease
progression and decrease healthcare expenditure [8, 28]. However, due to the
lack of routine screening, COPD is still widely underdiagnosed [14]. In
clinical practice, pulmonary function tests (PFT) and spirometry tests are
necessary diagnostic tools for COPD diagnosis [45]. Nevertheless, spirometry
tests may fail to capture the presence of early COPD before it becomes
detectable. Therefore asymptomatic patients are less likely to be tested [2,
23, 43]. Further, the current screening tools, such as spirometry, are
expensive and scarce in low-to-middle-income countries, leading to a delay in
diagnosis [4, 26]. Instead, the chest radiograph (CXR) is ubiquitous and cheap
[27]. The ability to create an early diagnostic screening tool for COPD from
CXRs offers a significant clinical opportunity in developing a tool for COPD
detection and thus targeting individuals for early interventions such as lung
cancer screening and smoking cessation programs.

Recent studies have demonstrated the ability of deep learning (DL) models,
such as convolutional neural networks (CNN), to classify radiographic findings
from CXRs and achieve human-like performance [24, 32, 41, 46], while also
processing large amounts of images at high speed. Previous studies have mostly
focused on classifying and detecting the 14 common chest radiographic findings
seen on CXRs obtained from radiology reports using publicly sourced databases
[32, 34, 48, 50]. Despite the success of DL in pulmonary disease
classification using CXRs, a very limited number of studies have explored the
potential of DL techniques in COPD diagnosis using CXRs only.

To the best of our knowledge, there is little attempt to develop a COPD
prognostic model to identify the presence of COPD in early stages using both
CXR and EHRs. In our study, we propose to build deep learning techniques using
model-level and data-level fusion strategies to diagnose COPD using real-world
data such as CheXpert, Medical Information Mart for Intensive Care (MIMIC-
IV)[20], MIMIC-CXR-JPG [22] and Emory data. Identifying COPD from CXRs and
EHRs for early COPD screening could prompt early diagnosis, prevention, and
treatment.

Our contributions are three-fold:

  * •

We implement model-level fusion strategies by combining four base models,
including DenseNet121, ResNet50V2, Xception, and MobileNetV2 to build a
primary screening model for COPD using real-world CXRs. To our knowledge, we
are among the first to focus on model building for early screening of COPD.

  * •

We use two data-level fusion strategies, including a multi-modal data
strategy. (MIMIC-CXRs and MIMIC-IV EHR) to incorporate relevant clinical
factors with CXRs and multi-site data strategy (MIMIC-CXRs and Emory CXRs) to
improve the performance and generalizability of the base models.

  * •

We assess model bias by evaluating the model performance of the data-level and
model-level strategies across different demographic subgroups, including race-
ethnicity, sex, and age.

##  2 Related Work

###  2.1 COPD diagnostic models

DL strategies, such as convolutional neural network (CNN) classifiers, have
become a dominant approach in classifying radiographic findings from CXRs.
Recent studies have sought to adapt various DL techniques to predict COPD;
however, these studies present various limitations. They are limited to a very
small number of COPD patients in both the training and testing set, only one
data type [30], and a single-institution [35]. Further, the studies are
limited in their consideration of the mechanical ventilation status of
patients and thus pose uncertainty about the stage of COPD for early
screening.

[9] demonstrated that the artificial intelligence system would learn spurious
features as shortcuts to aid their classification rather than learning
pathologically relevant features. For this concern, in [35]’s study, the
intubation may become a shortcut for models to determine the stage of COPD
because the intubation would be obvious in the CXRs and indicative of
respiratory status. To the best of our knowledge, the latent ability to
recognize patients with early COPD and without mechanical ventilation solely
depending on CXRs has not yet been investigated.

###  2.2 Model-level and data-level fusion

Two types of fusion strategies have shown to be promising in prediction tasks.
Model-level fusion [10] is one approach which combines base-model predictions
to form a composite prediction [1, 6, 12, 31] and thus leverages the strength
of each model. Bagging and stacking are two popular model-fusion techniques
[12]. Bagging integrates a set of predictions from several independent base
models into a single prediction using an averaging or majority voting scheme.
Stacking is an approach where the meta-model learns how to combine the output
of the base models best to provide an optimal prediction. Model-level fusion
schemes have been implemented in different medical imaging tasks and have been
shown to perform better and more stable than single model architectures [1, 3,
7, 19, 33, 39].

Data-level fusion schemes integrate multi-modal features or multi-site data
into a single predictor [15, 16, 42]. Various studies have demonstrated the
consistent improvement in model performance using data-level fusion strategies
[15, 16, 29, 42, 35]. With the proliferation of publicly available medical
imaging datasets and electronic health records (EHR), it is important to
utilize pixel-level data and clinical patient data to obtain better feature
representations during training. It is also important to use data from
multiple centers to account for variability in data acquisition and processing
to increase the validity of the proposed fusion method [11].

###  2.3 Fairness evaluation across different subgroups

Fairness is a rising concern in DL applications in the medical field. [38]
demonstrated that the state-of-the-art DL model consistently underdiagnose
under-served patients. In [37]’s study, the researchers found disparities in
true positive rate (TPR) for pulmonary disease diagnosis across different
demographic subgroups. Because TPR requires a binary threshold which would
significantly affect the calculation of TPR, we use AUC as our fairness
evaluation metric instead.

Yet, there are only a few studies [47] investigating the evaluation of DL
model on different patient subgroups and comparing it using real-world
multimodal data. Our work aims to asses bias across different demographic
groups to examine if our models yield biased prediction results. We evaluate
the models across different demographic subgroups and calculate the AUC’s
standard deviation (SD) for each, respectively, to quantify the discrepancies
in the AUC across the models. We further compare the SD of the AUC for the
Xcpetion base model to that of the fusion strategies to see if the fusion
strategies would decrease the SD of the AUC in order to mitigate the
discrepancies across the models.

##  3 Method

###  3.1 Data

####  3.1.1 Chest X-Ray Datasets

This study uses frontal CXRs from three datasets: CheXpert, MIMIC-CXR-JPG, and
Emory-CXR. CheXpert is a retrospective dataset from the Stanford Hospital
consisting of 224,316 CXRs of 65,240 patients [18]. MIMIC-CXR-JPG is a
publically available dataset containing 377,110 JPG images and its structured
labels corresponding to 227,835 radiographic studies sourced from the Beth
Israel Deaconess Medical Center between 2011 and 2016 [13, 21, 22]. Both
CheXpert and MIMIC-CXR are labeled with the same 14 cardiopulmonary disease
labels. Emory-CXR is a private dataset collected from five hospitals across
the Emory Healthcare system between 2019 and 2020. This dataset is acquired
from inpatient and outpatient hospitals and contains 226,640 CXRs from 57,909
patients.

####  3.1.2 Electronic Health Record (EHR) databases

We use the MIMIC-IV database (version: 1.0), which is a publicly available EHR
database maintained by Beth Israel Deaconess Medical Center from 2008 to 2019,
consisting of more than 200,000 emergency department (ED) admissions and more
than 60,000 intensive care units (ICU) stays [13, 20]. We also use the private
EHR from Emory’s institutional database, collected from five hospitals across
the Emory Healthcare system.

Table 1: Summary of the datasets used for the COPD prediction experiments. | MIMIC | Emory | CheXpert
---|---|---|---
Number of patients (number of images) | 52,804 (201,748) | 57,909 (226,640) | 48,285 (140,921)
Asian | 1,930 (3.7%) | 2,035 (3.5%) | 7,422 (15.4%)
Black | 8,923 (16.9%) | 26,729 (46.2%) | 3,016 (6.2%)
Latino | 3,340 (6.3%) | 1,006 (1.7%) | 2,322 (4.8%)
Others | 4,720 (8.9%) | 3,199 (5.5%) | N/A
White | 33,891 (64.2%) | 24,940 (43.1%) | 34,657 (71.8%)
Female | 27,473 (52%) | 30,236 (52.2%) | 21,603 (44.7%)
Male | 25,331 (48%) | 27,673 (47.8%) | 26,182 (55.3%)
0-40 | 8,216 (15.6%) | 11,691 (20.2%) | 5,880 (12.4%)
40-60 | 16,433 (31.1%) | 16,909 (29.2%) | 13,565 (28.6%)
60-80 | 19,493 (36.9%) | 22,343 (38.6%) | 18,339 (38.7%)
80+ | 8,662 (16.4%) | 6,966 (12%) | 9,633 (20.3%)

####  3.1.3 Patient Cohort and Data Pre-processing

We use the MIMIC-IV database to extract the patient cohort with COPD and those
without COPD from both the ED and ICU based on the ICD diagnoses. We exclude
images from patients who were on mechanical ventilation for two reasons.
First, we did not want the algorithm to learn from the CXRs that patients
intubated with a mechanical ventilation device are indicative of COPD. Second,
our focus is to build an early-screening tool for patients with COPD, and so
we exclude COPD patients who are on mechanical ventilation as they are more
likely to have an advanced disease. We used images from patients on room air,
oxygen, or high-flow nasal cannula. We create binary labels for the CXRs using
the ICD-9 and ICD-10 diagnosis codes from the MIMIC-IV and Emory EHR listed in
Supplementary Table A1 to indicate the presence or absence of COPD.

The patient’s demographic extracted from MIMIC-IV and Emory EHR includes self-
reported race-ethnicity, sex, and age. The summary of the datasets we use in
our experiments is shown in Table 1. All images are processed with histogram
equalization, resized to (256, 256). The image pixel values are normalized
from 0 to 1. We split the imaging data into training (64%), validation (16%),
and testing (20%) based on the patient’s ID, and we ensure that patients in
the training set are not included in the validation and the testing set to
avoid data leakage.

###  3.2 Transfer learning

This study uses four CNN-derived state-of-the-art models, including
DenseNet121, ResNet50V2, MobileNetV2, and Xception equipped with ImageNet pre-
trained weights obtained from the TensorFlow Keras library. In the pre-
training stage, we train the four base models to classify six radiographic
labels using the CheXpert dataset to allow the models to learn pulmonary
features. We select top six radiographic labels, including ’Atelectasisa’,
’Cardiomegaly’, ’Edema’, ’Lung Opacity’, ’Pleural Effusion’, and ’Support
Devices’, out of fourteen common labels based on their prevalence in the
CheXpert dataset. In the fine-tuning stage, each pre-trained model is fine-
tuned on the MIMIC-CXR dataset for the COPD detection task. To accelerate the
fine-tuning process, we freeze the learned weights for the first 30% of the
layers and fine-tune the remaining layers. To address the imbalance in the
number of cases of COPD, we use the class weights parameter to give different
weightage to the COPD and non-COPD classes. Since the class weights regularize
the loss function, the model would give a higher weight for the loss of the
COPD class while a lower weight for the loss of the non-COPD class. We
calculate the class weights for the COPD and non- COPD class, which are
defined as follows:

$w_{COPD}=\frac{\text{total \\# of CXRs}}{2*(\text{\\# of COPD CXRs)}}$ (1)
$w_{non-COPD}=\frac{\text{total \\# of CXRs}}{2*(\text{\\# of non-COPD
CXRs)}}$ (2)

Figure 1: Overview of our proposed framework.

We apply 1.91 and 0.68 weights to the loss of the COPD and non-COPD classes,
respectively. For the pre-training and fine-tuning stages, the learning rate
is set to 0.001 and decays 5% for every epoch, and the batch size is set to
32. Adam optimizer is used, and binary cross entropy is used as the loss
function. We stop the model training when the validation loss does not improve
for three consecutive epochs. All hyperparameters and random seeds are
identical for all training processes to ensure reproducibility.

###  3.3 Fusion strategies

To improve the performance of the fine-tuned models, we use the fine-tuned
models as the base models to implement model-level fusion and data-level
fusion strategies. In terms of model-level fusion, we implement bagging and
stacking methods. For the bagging methods, we implement unweighted and
weighted averaging. The stacking method contains the meta-models trained using
the predictions of the base models. For data-level fusion, we experiment with
two different techniques. First, we merge the MIMIC-CXRs and Emory-CXRs to
create a multi-site data fusion strategy, which we use as our training set to
train the best-performing base model (Xception). Second, we implement a multi-
modal data fusion strategy ( joint fusion type II) [15], where we concatenate
the feature representations of the CXRs from the penultimate Xception model
layer with demographic variables from EHR. Our proposed approaches are shown
in Figure 1.

####  3.3.1 Model-level fusion

##### Unweighted average bagging

We average the predictions of the four fine-tuned base models with equal
weighting to obtain the final predictions. We later evaluate whether the
unweighted average bagging strategy outperform the single base models.

##### Weighted average bagging

We calculate the cosine distance of the validation set predictions obtained
from the four base models and construct a dendrogram using the Unweighted Pair
Group Method with an Arithmetic mean (UPGMA) algorithm, an hierarchical
clustering approach [44] defined as follows:

$d(u,v)=\sum_{ij}{}\frac{d(u[i],v[i])}{|u|*|v|}\text{,}$ (3)

where $|u|$ and $|v|$ are the cardinalities of clusters u and v, respectively
for all points i and j. We average the test predictions of the base models
based on the dendrogram as shown in Figure 2 to obtain the final prediction.
For example, in Figure 2, the predictions of the DenseNet121 and the
ResNet50V2 are averaged first, followed by the predictions of the MobileNetV2,
and finally, the predictions of the Xception. This is equivalent to assigning
weights of the predictions of the DenseNet121, ResNet50V2, MobileNetV2, and
Xception as $\frac{1}{8}$, $\frac{1}{8}$, $\frac{1}{4}$, and $\frac{1}{2}$,
respectively. We give the prediction of the base model a larger weight since
the length of the dendrogram leg for the base model is longer.

Figure 2: Dendrogram showing pairwise cosine distance for predictions across
the four base models.

##### Base model stacking

We use the training data predictions of the base models as inputs to train the
meta-model. The meta-model would optimize predictions based on the test data
predictions from the base models. We implement four machine learning models as
the meta-models, including logistic regression, k-nearest neighbors, XGBoost,
and random forest.

####  3.3.2 Data-level fusion

##### Multi-site data fusion

We merge two datasets, MIMIC-CXR and Emory-CXR, and use the combined dataset
to fine-tune and evaluate the best-performing base model (Xception). We
examine if data from different sources enhances the performance and
generalizability of the multi-site data fusion strategy.

##### Multi-modal data fusion

We construct a data-level fusion strategy with the CXRs and their
corresponding EHRs as input. We use an Xception base model and a multi-layer
perceptron model to create a multi-modal strategy (joint fusion—type II) [15].
The learned features representations of the CXRs from the penultimate fine-
tuned Xception layer are encoded into 64 hidden nodes concatenated with ten
input nodes representing the demographic variables —five for race-ethnicity
(Asian, Black, Hispanic, Others, and White), four for age (0-40, 40-60, 60-80,
and 80+), and one for sex (Female or Male).

###  3.4 Evaluation metrics

We report the point estimates and the 95% confidence interval (CI) of the
evaluation metrics: area under the receiver operating characteristic curve
(AUC), precision, recall, F1-score, the area under the precision-recall curve
(AUPRC), and balanced accuracy [5] for all the experimented and proposed
methodologies. 95% (CI) is calculated by bootstrapping the metric scores over
1000 runs.

###  3.5 Model interpretation

To understand which CXR regions the base models focus on, we visualize the hot
spots for the true positive by using the Gradient-weighted Class Activation
Mapping (GradCAM) algorithm [36].

###  3.6 Fairness analysis

To examine whether there is a discrepancy in the models performance across
demographic subgroups, we evaluate the average AUC and the standard deviation
(SD) of the AUCs among different race-ethnicity, age groups, and sex
separately on the MIMIC and Emory datasets. That is, the model would have
separate results for each demographic group. The average AUC shows how the
models perform in the demographic group, and the SD shows the deviation from
the average AUC across demographic subgroups.

##  4 Results

###  4.1 Transfer learning

We present the results of four different state-of-the-art CNN models tested on
MIMIC-CXR and Emory-CXR in Table 3 and Table 4 of Appendix A. The Xception
base model yields the best AUC score (0.76), F1-score (0.54), AUPRC (0.54),
and balanced accuracy (0.69) during the internal testing on MIMIC-CXR;
however, the performance of the four base models is relatively similar.
Further, when performing external testing on Emory-CXR, the Xception base
model attains the best performance. Additionally, the performance of the
Xception base model on both MIMIC-CXR and Emory-CXR is very similar,
demonstrating that our fine-tuned model does not overfit on the training data.
The detailed results with precision and recall are listed in Supplementary
Table B1 to B2.

###  4.2 Fusion strategies

We present the model testing results of our fusion strategies in Table 2 for
MIMIC-CXR and Table 5 of Appendix B for Emory-CXR. For testing on MIMIC-CXR,
the unweighted and weighted bagging models have the best AUC score (0.77),
F1-score (0.54), AUPRC (0.56), and balanced accuracy (0.70). Since the
Xception model performs the best among the base models, we use the Xception
model in the data-level fusion strategy. The stacking model using the random
forest as a meta-model performs the best among the other stacking models;
however, it does not outperform the Xception basemodel.

The multi-site data-level fusion strategy results show that merging the Emory-
CXR dataset with MIMIC-CXR for training does not improve the performance when
testing on MIMIC-CXR. Nonetheless, it performs better when testing on Emory-
CXR than the Xception base model. The multi-modal data-level fusion strategy
results show that incorporating demographic features does not improve the
performance of both testing datasets. We use the chi-square test to evaluate
the independence of the race-ethnicity, age, sex variables, and the 14
radiographic labels between the COPD and non-COPD cohorts. A lower p-value
($p<0.01$) demonstrates that the distributions of these variables are
statistically significant between the COPD and non-COPD patient cohorts. The
detailed results with precision and recall are listed in Supplementary Table
B3 to B4.

Table 2: Evaluation metrics averaged over 1000 epochs ±95% for the various fusion strategies on MIMIC-CXR test data compared to the Xception base model. To validate our implementation, we shuffled the images’ labels and used the randomly labeled images to train a model as a permutation test. The results of the permutation test indicated our implementation is sound. Fusion Strategies | AUC | F1-score | AUPRC | Balanced Accuracy
---|---|---|---|---
Base Xception Model | 0.76 [0.75-0.76] | 0.53 [0.52-0.54] | 0.54 [0.53-0.55] | 0.69 [0.68-0.70]
Base Models Bagging |  |  |  |
Unweighted Average | 0.77 [0.76-0.77] | 0.54 [0.53-0.55] | 0.56 [0.55-0.57] | 0.70 [0.69-0.70]
Weighted Average | 0.77 [0.76-0.77] | 0.54 [0.53-0.55] | 0.56 [0.55-0.57] | 0.70 [0.69-0.70]
Base Models Stacking |  |  |  |
Logistic Regression | 0.76 [0.75-0.76] | 0.53 [0.52-0.54] | 0.54 [0.53-0.55] | 0.69 [0.68-0.69]
XGBoost | 0.75 [0.75-0.76] | 0.53 [0.52-0.54] | 0.54 [0.53-0.55] | 0.69 [0.68-0.69]
K-nearest Neighbors | 0.75 [0.74-0.75] | 0.52 [0.52-0.53] | 0.52 [0.51-0.53] | 0.68 [0.68-0.69]
Random Forest | 0.76 [0.76-0.77] | 0.54 [0.53-0.55] | 0.55 [0.54-0.56] | 0.69 [0.69-0.70]
Data-level Fusion |  |  |  |
Multi-site | 0.75 [0.75-0.76] | 0.53 [0.52-0.54] | 0.54 [0.53-0.55] | 0.69 [0.68-0.69]
Multi-modal | 0.76 [0.75-0.77] | 0.54 [0.53-0.54] | 0.55 [0.53-0.56] | 0.69 [0.69-0.70]
Permutation test | 0.51 [0.50-0.52] | 0.31 [0.30-0.31] | 0.19 [0.18-0.19] | 0.51 [0.51-0.52]

##  5 Discussion and Conclusion

Our findings support our hypothesis that DL models trained using CXRs can
identify patients with early-stage COPD. CXRs are a relatively accessible and
affordable imaging modality. Therefore, our approach could provide a low-cost
and widely available screening opportunity for COPD. We demonstrate that the
Xception base model detects COPD with reasonable performance compared to the
other base models as shown in Table 3 and Table 4 of Appendix A. The results
are consistent for both MIMIC and Emory test datasets across the fusion
strategies. The multi-site data-level fusion scheme perfroms better than the
other schemes on the Emory-CXR test data indicating that our proposed strategy
is generalizable. We use Grad-CAM to evaluate the Xception model’s
explainability. Looking at the true positive cases for both MIMIC-CXR and
Emory-CXR datasets, we can see that the Grad-CAM of the Xception base model
focuses on the upper right lung area.

The limitations of this study include reliance on ICD9/10 coding as labels and
the lack of Pulmonary Function Testing (PFT) labels. PFT is an effective tool
for the diagnosis of COPD. However, our study uses ICD codes to identify COPD
diagnosis because concurrent PFT data was unavailable for many of the CXRs in
our datasets. Instead, our study is utilized as a proof of concept that deep
learning can be applied to patient data even if gold-standard PFT findings are
not available. This has widespread implications because the current diagnostic
criteria for COPD leads to a large population of underdiagnosed patients [33,
49]. Additionally, clinicians often only order PFTs when COPD progresses and
becomes symptomatic, indicating that many early-stage diagnoses are often
missed [40]. In the future, we will incorporate PFT findings to improve the
performance of our model. We encourage hospitals to build their COPD screening
model using their database based on our open-source code. Our proposed
approach is advantageous for the COPD screening as it only requires CXRs,
which are routinely collected and easily accessible.

###  5.1 Explainability with Grad-CAM

We visualize the Grad-CAM heatmaps of the true positive cases on the Xception
base model for both MIMIC-CXR and Emory-CXR test datasets. Table 6 of C.1
shows the heatmaps of the true positive cases of the Xception model. We
further average the heatmaps with predicted probabilities higher than 0.9 to
obtain an averaged heatmap (Mean). We can see that the Xception base model
mainly focuses on the upper right side of the lungs between the upper lung and
the trachea, illustrating that COPD is associated with the upper airway [17].

###  5.2 Fairness analysis

We evaluate the performance of the fusion strategies on each demographic group
separately in both datasets. Table 7 and Table 8 show the results of the
fairness analysis of each model on MIMIC-CXR and Emory-CXR. The average AUC
scores show the average performance of the models among each demographic
group, and the standard deviations show the performance gap between each
demographic group. The star signs in the table represent the highest standard
deviation, indicating the model’s unfairness. The baseline and weighted
average bagging models are the only two that do not have a star sign. That is,
both models are reasonably fair across demographic groups. More importantly,
since the weighted average bagging model has a higher average AUC, we can
conclude that it is the best model that combines both reasonable performance
and fairness evaluation across various demographic groups compared to other
fusion strategies. The detailed results with precision and recall are listed
in Supplementary Table B5 to B18.

## References

  * [1] F. Ahmad, A. Farooq, and M. U. Ghani.  Deep Ensemble Model for Classification of Novel Coronavirus in Chest X-Ray Images.  Computational Intelligence and Neuroscience, 2021:e8890226, Jan. 2021.  Publisher: Hindawi.
  * [2] E. Andreeva, M. Pokhaznikova, A. Lebedev, I. Moiseeva, O. Kuznetsova, and J.-M. Degryse.  Spirometry is not enough to diagnose COPD in epidemiological studies: a follow-up study.  npj Primary Care Respiratory Medicine, 27(1):1–9, Nov. 2017.  Number: 1 Publisher: Nature Publishing Group.
  * [3] A. Barragán-Montero, U. Javaid, G. Valdés, D. Nguyen, P. Desbordes, B. Macq, S. Willems, L. Vandewinckele, M. Holmström, F. Löfman, S. Michiels, K. Souris, E. Sterpin, and J. A. Lee.  Artificial intelligence and machine learning for medical imaging: A technology review.  Physica Medica, 83:242–256, Mar. 2021.
  * [4] D. Beran, H. J. Zar, C. Perrin, A. M. Menezes, and P. Burney.  Burden of asthma and chronic obstructive pulmonary disease and access to essential medicines in low-income and middle-income countries.  The Lancet Respiratory Medicine, 3(2):159–170, Feb. 2015.
  * [5] K. H. Brodersen, C. S. Ong, K. E. Stephan, and J. M. Buhmann.  The Balanced Accuracy and Its Posterior Distribution.  In 2010 20th International Conference on Pattern Recognition, pages 3121–3124, Aug. 2010.  ISSN: 1051-4651.
  * [6] Y. Cao, T. A. Geddes, J. Y. H. Yang, and P. Yang.  Ensemble deep learning in bioinformatics.  Nature Machine Intelligence, 2(9):500–508, Sept. 2020.  Number: 9 Publisher: Nature Publishing Group.
  * [7] T. B. Chandra, K. Verma, B. K. Singh, D. Jain, and S. S. Netam.  Coronavirus disease (COVID-19) detection in Chest X-Ray images using majority voting based classifier ensemble.  Expert Systems with Applications, 165:113909, Mar. 2021.
  * [8] M. Decramer and C. B. Cooper.  Treatment of COPD: the sooner the better?  Thorax, 65(9):837–841, Sept. 2010.  Publisher: BMJ Publishing Group Ltd Section: Review.
  * [9] A. J. DeGrave, J. D. Janizek, and S.-I. Lee.  AI for radiographic COVID-19 detection selects shortcuts over signal.  Nature Machine Intelligence, 3(7):610–619, July 2021.  Bandiera_abtest: a Cg_type: Nature Research Journals Number: 7 Primary_atype: Research Publisher: Nature Publishing Group Subject_term: Computational science;Predictive medicine;Radiography;SARS-CoV-2 Subject_term_id: computational-science;predictive-medicine;radiography;sars-cov-2.
  * [10] K.-L. Du and M. N. S. Swamy.  Combining Multiple Learners: Data Fusion and Ensemble Learning.  In K.-L. Du and M. N. S. Swamy, editors, Neural Networks and Statistical Learning, pages 737–767. Springer, London, 2019.
  * [11] L. Flynn.  The benefits and challenges of multisite studies: Lessons learned.  AACN Advanced Critical Care, 20(4):388–391, 2009.
  * [12] M. A. Ganaie, M. Hu, A. K. Malik, M. Tanveer, and P. N. Suganthan.  Ensemble deep learning: A review.  Engineering Applications of Artificial Intelligence, 115:105151, Oct. 2022.  arXiv:2104.02395 [cs].
  * [13] A. L. Goldberger, L. A. Amaral, L. Glass, J. M. Hausdorff, P. C. Ivanov, R. G. Mark, J. E. Mietus, G. B. Moody, C.-K. Peng, and H. E. Stanley.  PhysioBank, PhysioToolkit, and PhysioNet: components of a new research resource for complex physiologic signals.  circulation, 101(23):e215–e220, 2000.  Publisher: Am Heart Assoc.
  * [14] K. Hill, R. S. Goldstein, G. H. Guyatt, M. Blouin, W. C. Tan, L. L. Davis, D. M. Heels-Ansdell, M. Erak, P. J. Bragaglia, I. E. Tamari, R. Hodder, and M. B. Stanbrook.  Prevalence and underdiagnosis of chronic obstructive pulmonary disease among patients at risk in primary care.  CMAJ: Canadian Medical Association journal = journal de l’Association medicale canadienne, 182(7):673–678, Apr. 2010.
  * [15] S.-C. Huang, A. Pareek, S. Seyyedi, I. Banerjee, and M. P. Lungren.  Fusion of medical imaging and electronic health records using deep learning: a systematic review and implementation guidelines.  npj Digital Medicine, 3(1):1–9, Oct. 2020.  Number: 1 Publisher: Nature Publishing Group.
  * [16] S.-C. Huang, A. Pareek, R. Zamanian, I. Banerjee, and M. P. Lungren.  Multimodal fusion with deep neural networks for leveraging CT imaging and electronic health record: a case-study in pulmonary embolism detection.  Scientific Reports, 10(1):22147, Dec. 2020.  Number: 1 Publisher: Nature Publishing Group.
  * [17] J. R. Hurst.  Upper airway. 3: Sinonasal involvement in chronic obstructive pulmonary disease.  Thorax, 65(1):85–90, Jan. 2010.
  * [18] J. Irvin, P. Rajpurkar, M. Ko, Y. Yu, S. Ciurea-Ilcus, C. Chute, H. Marklund, B. Haghgoo, R. Ball, and K. Shpanskaya.  Chexpert: A large chest radiograph dataset with uncertainty labels and expert comparison.  In Proceedings of the AAAI conference on artificial intelligence, volume 33, pages 590–597, 2019.  Issue: 01.
  * [19] H. Jia, Y. Xia, Y. Song, W. Cai, M. Fulham, and D. D. Feng.  Atlas registration and ensemble deep convolutional neural network-based prostate segmentation using magnetic resonance imaging.  Neurocomputing, 275:1358–1369, Jan. 2018.
  * [20] A. Johnson, L. Bulgarelli, T. Pollard, S. Horng, L. A. Celi, and R. Mark.  MIMIC-IV, 2020.  Version Number: 2.0 Type: dataset.
  * [21] A. E. Johnson, T. J. Pollard, S. J. Berkowitz, N. R. Greenbaum, M. P. Lungren, C.-y. Deng, R. G. Mark, and S. Horng.  MIMIC-CXR, a de-identified publicly available database of chest radiographs with free-text reports.  Scientific data, 6(1):1–8, 2019.  Publisher: Nature Publishing Group.
  * [22] A. E. W. Johnson, T. J. Pollard, N. R. Greenbaum, M. P. Lungren, C.-y. Deng, Y. Peng, Z. Lu, R. G. Mark, S. J. Berkowitz, and S. Horng.  MIMIC-CXR-JPG, a large publicly available database of labeled chest radiographs.  arXiv:1901.07042 [cs, eess], Nov. 2019.  arXiv: 1901.07042.
  * [23] L. E. Labonté, W. C. Tan, P. Z. Li, P. Mancino, S. D. Aaron, A. Benedetti, K. R. Chapman, R. Cowie, J. M. FitzGerald, P. Hernandez, F. Maltais, D. D. Marciniuk, D. O’Donnell, D. Sin, and J. Bourbeau.  Undiagnosed Chronic Obstructive Pulmonary Disease Contributes to the Burden of Health Care Use. Data from the CanCOLD Study.  American Journal of Respiratory and Critical Care Medicine, 194(3):285–298, Aug. 2016.  Publisher: American Thoracic Society - AJRCCM.
  * [24] P. Lakhani and B. Sundaram.  Deep Learning at Chest Radiography: Automated Classification of Pulmonary Tuberculosis by Using Convolutional Neural Networks.  Radiology, 284(2):574–582, Aug. 2017.  Publisher: Radiological Society of North America.
  * [25] D. M. Mannino and A. S. Buist.  Global burden of COPD: risk factors, prevalence, and future trends.  The Lancet, 370(9589):765–773, Sept. 2007.
  * [26] J. Meghji, K. Mortimer, A. Agusti, B. W. Allwood, I. Asher, E. D. Bateman, K. Bissell, C. E. Bolton, A. Bush, B. Celli, C.-Y. Chiang, A. A. Cruz, A.-T. Dinh-Xuan, A. El Sony, K. M. Fong, P. I. Fujiwara, M. Gaga, L. Garcia-Marcos, D. M. G. Halpin, J. R. Hurst, S. Jayasooriya, A. Kumar, M. V. Lopez-Varela, R. Masekela, B. H. Mbatchou Ngahane, M. Montes de Oca, N. Pearce, H. K. Reddel, S. Salvi, S. J. Singh, C. Varghese, C. F. Vogelmeier, P. Walker, H. J. Zar, and G. B. Marks.  Improving lung health in low-income and middle-income countries: from challenges to solutions.  The Lancet, 397(10277):928–940, Mar. 2021.
  * [27] F. A. Mettler, M. Mahesh, M. Bhargavan-Chatfield, C. E. Chambers, J. G. Elee, D. P. Frush, D. L. Miller, H. D. Royal, M. T. Milano, D. C. Spelic, A. J. Ansari, W. E. Bolch, G. M. Guebert, R. H. Sherrier, J. M. Smith, and R. J. Vetter.  Patient Exposure from Radiologic and Nuclear Medicine Procedures in the United States: Procedure Volume and Effective Dose for the Period 2006–2016.  Radiology, 295(2):418–427, May 2020.
  * [28] C. Moretz, Y. Zhou, A. D. Dhamane, K. Burslem, K. Saverno, G. Jain, G. Devercelli, S. Kaila, J. J. Ellis, G. Hernandez, and A. Renda.  Development and Validation of a Predictive Model to Identify Individuals Likely to Have Undiagnosed Chronic Obstructive Pulmonary Disease Using an Administrative Claims Database.  Journal of Managed Care & Specialty Pharmacy, 21(12):1149–1159, Dec. 2015.
  * [29] J. G. Nam, H.-R. Kang, S. M. Lee, H. Kim, C. Rhee, J. M. Goo, Y.-M. Oh, C.-H. Lee, and C. M. Park.  Deep Learning Prediction of Survival in Patients with Chronic Obstructive Pulmonary Disease Using Chest Radiographs.  Radiology, page 212071, June 2022.  Publisher: Radiological Society of North America.
  * [30] A. Pyrros, J. R. Fernandez, S. M. Borstelmann, A. Flanders, D. Wenzke, E. Hart, J. M. Horowitz, P. Nikolaidis, M. Willis, A. Chen, P. Cole, N. Siddiqui, M. Muzaffar, N. Muzaffar, J. McVean, M. Menchaca, A. K. Katsaggelos, S. Koyejo, and W. Galanter.  Validation of a deep learning, value-based care model to predict mortality and comorbidities from chest radiographs in COVID-19.  PLOS Digital Health, 1(8):e0000057, Aug. 2022.  Publisher: Public Library of Science.
  * [31] M. Ragab, K. Eljaaly, N. A. Alhakamy, H. A. Alhadrami, A. A. Bahaddad, S. M. Abo-Dahab, and E. M. Khalil.  Deep Ensemble Model for COVID-19 Diagnosis and Classification Using Chest CT Images.  Biology, 11(1):43, Jan. 2022.  Number: 1 Publisher: Multidisciplinary Digital Publishing Institute.
  * [32] P. Rajpurkar, J. Irvin, K. Zhu, B. Yang, H. Mehta, T. Duan, D. Ding, A. Bagul, C. Langlotz, K. Shpanskaya, M. P. Lungren, and A. Y. Ng.  CheXNet: Radiologist-Level Pneumonia Detection on Chest X-Rays with Deep Learning, Dec. 2017.  arXiv:1711.05225 [cs, stat].
  * [33] E. A. Regan, D. A. Lynch, D. Curran-Everett, J. L. Curtis, J. H. M. Austin, P. A. Grenier, H.-U. Kauczor, W. C. Bailey, D. L. DeMeo, R. H. Casaburi, P. Friedman, E. J. R. Van Beek, J. E. Hokanson, R. P. Bowler, T. H. Beaty, G. R. Washko, M. K. Han, V. Kim, S. S. Kim, K. Yagihashi, L. Washington, C. E. McEvoy, C. Tanner, D. M. Mannino, B. J. Make, E. K. Silverman, and J. D. Crapo.  Clinical and Radiologic Disease in Smokers With Normal Spirometry.  JAMA internal medicine, 175(9):1539–1549, Sept. 2015.
  * [34] J. Rubin, D. Sanghavi, C. Zhao, K. Lee, A. Qadir, and M. Xu-Wilson.  Large Scale Automated Reading of Frontal and Lateral Chest X-Rays using Dual Convolutional Neural Networks, Apr. 2018\.  arXiv:1804.07839 [cs, stat].
  * [35] J. D. Schroeder, R. Bigolin Lanfredi, T. Li, J. Chan, C. Vachet, R. Paine III, V. Srikumar, and T. Tasdizen.  Prediction of Obstructive Lung Disease from Chest Radiographs via Deep Learning Trained on Pulmonary Function Data.  International Journal of Chronic Obstructive Pulmonary Disease, 15:3455–3466, Jan. 2021.
  * [36] R. R. Selvaraju, M. Cogswell, A. Das, R. Vedantam, D. Parikh, and D. Batra.  Grad-cam: Visual explanations from deep networks via gradient-based localization.  In Proceedings of the IEEE international conference on computer vision, pages 618–626, 2017.
  * [37] L. Seyyed-Kalantari, G. Liu, M. McDermott, I. Y. Chen, and M. Ghassemi.  Chexclusion: Fairness gaps in deep chest x-ray classifiers.  In BIOCOMPUTING 2021: proceedings of the Pacific symposium, pages 232–243. World Scientific, 2020.
  * [38] L. Seyyed-Kalantari, H. Zhang, M. McDermott, I. Y. Chen, and M. Ghassemi.  Underdiagnosis bias of artificial intelligence algorithms applied to chest radiographs in under-served patient populations.  Nature medicine, 27(12):2176–2182, 2021.
  * [39] I. Sirazitdinov, M. Kholiavchenko, T. Mustafaev, Y. Yuan, R. Kuleev, and B. Ibragimov.  Deep neural network ensemble for pneumonia localization from a large-scale chest x-ray database.  Computers & Electrical Engineering, 78:388–399, Sept. 2019.
  * [40] A. Sood, H. Petersen, C. Qualls, P. M. Meek, R. Vazquez-Guillamet, B. R. Celli, and Y. Tesfaigzi.  Spirometric variability in smokers: transitions in copd diagnosis in a five-year longitudinal study.  Respiratory research, 17(1):1–10, 2016.
  * [41] Y.-X. Tang, Y.-B. Tang, Y. Peng, K. Yan, M. Bagheri, B. A. Redd, C. J. Brandon, Z. Lu, M. Han, J. Xiao, and R. M. Summers.  Automated abnormality classification of chest radiographs using deep convolutional neural networks.  npj Digital Medicine, 3(1):1–8, May 2020.  Number: 1 Publisher: Nature Publishing Group.
  * [42] A. Tariq, L. A. Celi, J. M. Newsome, S. Purkayastha, N. K. Bhatia, H. Trivedi, J. W. Gichoya, and I. Banerjee.  Patient-specific COVID-19 resource utilization prediction using fusion AI model.  npj Digital Medicine, 4(1):1–9, June 2021.  Number: 1 Publisher: Nature Publishing Group.
  * [43] US Preventive Services Task Force (USPSTF).  Screening for Chronic Obstructive Pulmonary Disease: US Preventive Services Task Force Recommendation Statement.  JAMA, 315(13):1372–1377, Apr. 2016.
  * [44] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  Nature Methods, 17:261–272, 2020.
  * [45] C. F. Vogelmeier, G. J. Criner, F. J. Martinez, A. Anzueto, P. J. Barnes, J. Bourbeau, B. R. Celli, R. Chen, M. Decramer, L. M. Fabbri, P. Frith, D. M. Halpin, M. V. López Varela, M. Nishimura, N. Roche, R. Rodriguez-Roisin, D. D. Sin, D. Singh, R. Stockley, J. Vestbo, J. A. Wedzicha, and A. Agusti.  Global Strategy for the Diagnosis, Management and Prevention of Chronic Obstructive Lung Disease 2017 Report: GOLD Executive Summary.  Respirology, 22(3):575–601, Apr. 2017.
  * [46] C. Wang, A. Elazab, J. Wu, and Q. Hu.  Lung nodule classification using deep feature fusion in chest radiography.  Computerized Medical Imaging and Graphics, 57:10–18, Apr. 2017\.
  * [47] R. Wang, P. Chaudhari, and C. Davatzikos.  Machine learning models are not necessarily biased when constructed properly: Evidence from neuroimaging studies, 2022.
  * [48] X. Wang, Y. Peng, L. Lu, Z. Lu, M. Bagheri, and R. M. Summers.  ChestX-ray8: Hospital-scale chest x-ray database and benchmarks on weakly-supervised classification and localization of common thorax diseases.  In 2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR). IEEE, jul 2017.
  * [49] P. G. Woodruff, R. G. Barr, E. Bleecker, S. A. Christenson, D. Couper, J. L. Curtis, N. A. Gouskova, N. N. Hansel, E. A. Hoffman, R. E. Kanner, E. Kleerup, S. C. Lazarus, F. J. Martinez, R. Paine, S. Rennard, D. P. Tashkin, and M. K. Han.  Clinical Significance of Symptoms in Smokers with Preserved Pulmonary Function.  New England Journal of Medicine, 374(19):1811–1821, May 2016.  Publisher: Massachusetts Medical Society.
  * [50] L. Yao, E. Poblenz, D. Dagunts, B. Covington, D. Bernard, and K. Lyman.  Learning to diagnose from scratch by exploiting dependencies among labels, Feb. 2018.  arXiv:1710.10501 [cs].

##  Appendix A Results of Base Models

In this section, we list the performance metrics on the base models for both
MIMIC-CXR and Emory-CXR data in 4.1.

Table 3: Base models tested on MIMIC-CXR data. Model | AUC | F1-score | AUPRC | Balanced Accuracy
---|---|---|---|---
DenseNet121 | 0.74 [0.74-0.75] | 0.52 [0.51-0.53] | 0.52 [0.51-0.53] | 0.68 [0.67-0.69]
MobileNetV2 | 0.75 [0.74-0.76] | 0.53 [0.52-0.53] | 0.53 [0.52-0.54] | 0.68 [0.68-0.69]
ResNet50V2 | 0.74 [0.73-0.75] | 0.52 [0.51-0.52] | 0.52 [0.51-0.53] | 0.67 [0.67-0.68]
Xception | 0.76 [0.75-0.76] | 0.53 [0.52-0.54] | 0.54 [0.53-0.55] | 0.69 [0.68-0.70]
Table 4: Base models tested on Emory-CXR data. Model | AUC | F1-score | AUPRC | Balanced Accuracy
---|---|---|---|---
DenseNet121 | 0.73 [0.72-0.73] | 0.43 [0.42-0.44] | 0.40 [0.38-0.41] | 0.66 [0.66-0.67]
MobileNetV2 | 0.72 [0.72-0.73] | 0.42 [0.41-0.43] | 0.39 [0.38-0.40] | 0.66 [0.65-0.67]
ResNet50V2 | 0.72 [0.71-0.73] | 0.42 [0.41-0.43] | 0.39 [0.38-0.40] | 0.66 [0.65-0.67]
Xception | 0.74 [0.73-0.75] | 0.43 [0.43-0.44] | 0.40 [0.39-0.42] | 0.67 [0.67-0.68]

##  Appendix B Results of Fusion Strategies

We include the results of a different test dataset, MIMIC-CXR, for the same
experiments presented in Table 5 of 5.2.

Table 5: Evaluation metrics averaged over 1000 epochs ±95% CI for the various fusion strategies on Emory-CXR test data compared to the Xception base model. Fusion Strategies | AUC | F1-score | AUPRC | Balanced Accuracy
---|---|---|---|---
Base Xception Model | 0.74 [0.73-0.75] | 0.43 [0.43-0.44] | 0.40 [0.39-0.42] | 0.67 [0.67-0.68]
Base Models Bagging |  |  |  |
Unweighted Average | 0.75 [0.74-0.75] | 0.44 [0.43-0.45] | 0.42 [0.41-0.43] | 0.68 [0.67-0.68]
Weighted Average | 0.75 [0.74-0.75] | 0.44 [0.43-0.45] | 0.42 [0.41-0.43] | 0.68 [0.67-0.68]
Base Models Stacking |  |  |  |
Logistic Regression | 0.74 [0.73-0.74] | 0.44 [0.43-0.45] | 0.41 [0.39-0.42] | 0.67 [0.66-0.68]
XGBoost | 0.73 [0.73-0.74] | 0.44 [0.43-0.45] | 0.40 [0.39-0.42] | 0.67 [0.66-0.68]
K-nearest Neighbors | 0.72 [0.72-0.73] | 0.43 [0.42-0.44] | 0.39 [0.38-0.40] | 0.67 [0.66-0.67]
Random Forest | 0.74 [0.74-0.75] | 0.44 [0.43-0.45] | 0.41 [0.40-0.43] | 0.68 [0.67-0.68]
Data-level Fusion |  |  |  |
Multi-site | 0.76 [0.75-0.76] | 0.45 [0.44-0.45] | 0.43 [0.42-0.44] | 0.68 [0.68-0.69]
Multi-modal | 0.74 [0.73-0.74] | 0.44 [0.43-0.45] | 0.41 [0.39-0.42] | 0.67 [0.67-0.68]

##  Appendix C Additional Experimental Results

###  C.1 Explainability with Grad-CAM

The model is highlighting areas commonly affected by COPD, such as the apex of
the lung fields or the areas of the mediastinum affiliated with the pulmonary
vasculature. The lung apices can often become emphysematous in COPD, and the
pulmonary vessels may become more prominent as pulmonary hypertension worsens
in the later stages of COPD. This is likely why the models are identifying
these regions to make predictions. We will add the interpretation of the
heatmap to the appendix so that the reviewers can understand the heatmap from
a clinical perspective.

Table 6: The first and second rows display the Grad-CAM heatmaps of the true positive cases on the Xception base model using the MIMIC-CXR test data and the Emory-CXR test data, respectively. Pred: 0.996 |  Pred: 0.996 |  Pred: 0.994 |  Pred: 0.992 |  Mean
---|---|---|---|---
|  |  |  |
Pred: 0.996 |  Pred: 0.996 |  Pred: 0.994 |  Pred: 0.992 |  Mean
|  |  |  |

###  C.2 Fairness Analysis

Table 7: Average AUC ±SD for model-level and data-level fusion strategies on MIMIC-CXR test dataset across various demographic subgroups. Fusion Strategies | Race-ethnicity | Sex | Age
---|---|---|---
Base Xception Model | 0.78 [0.026] | 0.76 [0.010] | 0.73 [0.036]
Unweighted Average Bagging | 0.79 [0.031] | 0.77 [0.010] | 0.74 [0.030]
Weighted Average Bagging | 0.79 [0.029] | 0.77 [0.010] | 0.74 [0.032]
Base Models Stacking (Random Forest Classifier) | 0.78 [0.030] |  0.77 [0.015] $\star$ | 0.73 [0.029]
Multi-site Data Fusion | 0.78 [0.025] | 0.76 [0.005] |  0.71 [0.046] $\star$
Multi-modal Data Fusion |  0.78 [0.032] $\star$ | 0.76 [0.010] | 0.72 [0.035]
Table 8: Average AUC ±SD for model-level and data-level fusion strategies on Emory-CXR test dataset across various demographic subgroups. Fusion Strategies | Race-ethnicity | Sex | Age
---|---|---|---
Base Xception Model | 0.72 [0.024] | 0.74 [0] | 0.70 [0.030]
Unweighted Average Bagging | 0.72 [0.038] |  0.75 [0.005] $\star$ | 0.71 [0.025]
Weighted Average Bagging | 0.73 [0.030] | 0.75 [0] | 0.71 [0.027]
Base Models Stacking (Random Forest Classifier) | 0.71 [0.038] |  0.74 [0.005] $\star$ | 0.70 [0.029]
Multi-site Data Fusion |  0.74 [0.050] $\star$ | 0.76 [0] |  0.72 [0.031] $\star$
Multi-modal Data Fusion | 0.72 [0.038] | 0.74 [0] | 0.69 [0.026]