Synergy theory for murine Harderian gland tumours after irradiation by mixtures of high-energy ionized atomic nuclei

Experimental studies reporting murine Harderian gland (HG) tumourigenesis have been a NASA concern for many years. Studies used particle accelerators to produce beams that, on beam entry, consist of a single isotope also present in the galactic cosmic ray (GCR) spectrum. In this paper synergy theory is described, potentially applicable to corresponding mixed-field experiments, in progress, planned, or hypothetical. The “obvious” simple effect additivity (SEA) approach of comparing an observed mixture dose–effect relationship (DER) to the sum of the components’ DERs is known from other fields of biology to be unreliable when the components’ DERs are highly curvilinear. Such curvilinearity may be present at low fluxes such as those used in the one-ion HG experiments due to non-targeted (‘bystander’) effects, in which case a replacement for SEA synergy theory is needed. This paper comprises in silico modeling of published experimental data using a recently introduced, arguably optimal, replacement for SEA: incremental effect additivity (IEA). Customized open-source software is used. IEA is based on computer numerical integration of non-linear ordinary differential equations. To illustrate IEA synergy theory, possible rapidly-sequential-beam mixture experiments are discussed, including tight 95% confidence intervals calculated by Monte-Carlo sampling from variance–covariance matrices. The importance of having matched one-ion and mixed-beam experiments is emphasized. Arguments are presented against NASA over-emphasizing accelerator experiments with mixed beams whose dosing protocols are standardized rather than being adjustable to take biological variability into account. It is currently unknown whether mixed GCR beams sometimes have statistically significant synergy for the carcinogenesis endpoint. Synergy would increase risks for prolonged astronaut voyages in interplanetary space.


Introduction
The galactic cosmic ray (GCR) radiation field in interplanetary space includes ions that have high proton charge number Z and high kinetic energy (HZE ions). These often have high relative biological effectiveness (RBE) for various kinds of biological damage (reviewed, e.g., in Cucinotta and Cacao 2017;Goodhead 2018). In particular, murine Harderian gland (HG) tumourigenesis induced, in accelerator beam experiments, by individual ions in the GCR spectrum, with some of the HZE having high RBEs, has long been a NASA concern (Fry et al. 1985;Curtis et al. 1992;Alpen et al. 1993Alpen et al. , 1994Edwards 2001;Chang et al. 2016;Norbury et al. 2016). The Chang et al. paper concerns experiments at the Brookhaven (NY) NASA Space Radiation Laboratory (NSRL). The earlier papers concern experiments using the BEVALAC accelerator at what was then called Lawrence Berkeley Laboratory (LBL). The Norbury et al. paper mentions both. If substantial synergy is found in mixed-field experiments that correspond to these one-ion experiments, that would increase NASA concerns. The present paper describes synergy theory applicable to such mixed-field experiments.

3
The data considered in the present paper consist of already published results. In silico calculations are tailored to analyze the first relevant experimental mixed-beam results, which will become available soon. Additional calculations discuss hypothetical mixed beams that illustrate some key points about synergy theory and risk estimation. Some of the present analysis extends to radiobiology ideas about synergy (reviewed, for example, in Foucquier and Guedj 2015) developed in pharmacometrics, toxicology, evolutionary ecology and other fields.
Dose-effect relations (DERs) will play a central role. Typically, the main information on mixed-beam components comes from one-ion DERs. The present study uses new DERs for the one-ion data. These are more parsimonious (i.e., have fewer adjustable parameters) than other recent models (Cucinotta et al. 2013b;Chang et al. 2016;Cucinotta and Cacao 2017) for the same data. Because of the parsimony, the new models are of possible interest in their own right. However, they are used in the present paper mainly because their relative simplicity facilitates synergy analysis of mixed radiation fields whose components have these one-ion DERs.
Given one-ion DERs, synergy theory results can be calculated. The question to be answered is whether mixture data from rapidly sequential exposures that approximate simultaneous exposure manifest synergy, antagonism, or neither. Importantly, synergy theory is applicable only to mixtures where the experimental conditions closely match the oneion experiments that led to the one-ion DERs (Berenbaum 1989;Ham et al. 2018). For example, the one-ion experiments considered in the present paper did not intentionally add shielding between upstream beam entry and the target, so any mixtures which have shielding intentionally added cannot usefully be considered. Similarly, almost all the oneion experiments involved acute rather than protracted irradiation, so of necessity only experiments where the total mixture dose is applied as rapidly as possible, within less than 20 min even for the most complicated mixtures analyzed here, will be considered.
Synergy theory compares an experimentally observed mixture DER with a calculated baseline mixture DER defining the absence of synergy and absence of antagonism. Researchers in pharmacology and toxicology have known for a very long time (Fraser 1872;Loewe and Muischnek 1926) that the naive method of analyzing mixture effects with the simple effect additivity (SEA) approach to synergynamely just adding component effects, as is often taken for granted as obviously correct-is actually unreliable unless each mixture component one-agent DER is approximately linear-no-threshold (LNT). This problem is reviewed, e.g., in Zaider and Rossi (1980), Berenbaum (1989), Geary (2013), Foucquier and Guedj (2015), Piggott et al. (2015) and Tang et al. (2015).
As a simple example of this unreliability, consider two hypothetical HZE beams with respective one-ion DERs E 1 = βd 1 2 and E 2 = 2βd 2 2 , where β is a positive constant. These one-ion DERs, shown in Fig. 1, are curvilinear, since the second derivative (i.e., 2β or 4β, respectively) is positive. Suppose there is a 50-50 mixture of the two beams, so that d 1 = d/2 and d 2 = d/2, where d is the total mixture dose. Then, since (d/2) 2 = d 2 /4 the simple effect additivity (SEA) baseline no synergy and no antagonism (NSNA) effect height is only half the average of the two effect heights, instead of lying between the two. Any sensible synergy criterion would consider the dashed curve in Fig. 1 as specifying unexpectedly small mixture effects, i.e., antagonism, rather than specifying the baseline NSNA DER. Such discrepancies often arise when, as here, mixture component one-ion DERs are highly curvilinear. Consequently, many different replacements for SEA synergy theory are now in use to plan and interpret mixture experiments.
At sufficiently small radiation doses and high LETs, only a small fraction of all cell nuclei suffer a direct hit by a radiation track (Curtis et al. 1992;Hanin and Zaider 2014). Non-targeted effects (NTE) are then sometimes important (Cucinotta and Chappell 2010;Cucinotta et al. 2013a;Hada et al. 2014;Cacao et al. 2016;Chang et al. 2016;Cucinotta and Cacao 2017;Shuryak 2017), with cells directly hit by an ion influencing nearby cells through inter-cellular signaling (Hatzi et al. 2015). However, the question of whether NTE are significantly carcinogenic at very low HZE doses remains open (Piotrowski et al. 2017). Models of NTE action that are smooth (i.e., have continuous derivatives of all orders) use one-ion DERs that are very curvilinear, with negative second derivative, at low doses (Brenner et al. 2001 to investigate mixtures whose component one-ion DERs take NTE into account. The present paper will use a replacement for SEA theory called incremental effect additivity (IEA). IEA theory was introduced in two recent papers (Siranart et al. 2016;Ham et al. 2018). "Incremental" refers to the fact that one-ion DER slopes play an essential role in the theory. The underlying idea was suggested by Lam (1987). A one-ion DER slope defines a linear relation between a sufficiently small dose increment and the corresponding effect increment. Thus, by analyzing sufficiently small increments, one can circumvent the curvilinearities that plague SEA synergy theory. A systematic analysis of slopes requires using ordinary differential equations (ODEs). Implementing Lam's insight using ODEs has become practical because computers have become very adept at integrating non-linear ODEs. However, Lam (1994) did not use ODEs in his proposed replacement, called "independent action", for SEA.
The two recent references (Siranart et al. 2016;Ham et al. 2018) that introduced IEA, after reviewing many other synergy theories, gave evidence that IEA synergy theory is probably the optimal substitute for SEA synergy theory.
One important advantage IEA synergy theory has over SEA synergy theory and over many commonly used replacements for SEA synergy theory is that IEA synergy theory obeys a mixture of mixtures ("mixmix") principle (Ham et al. 2018). This principle is quite important for the purposes of the present paper because even nominally one-ion accelerator beams are often mixtures of different radiation qualities when they strike a target such as the HG, due to interaction of the primary radiation with animal tissue (resulting in self-shielding) or with other matter in the beam. Another advantage is that the IEA theory can be applied to mixtures whose one-ion component DERs have very heterogeneous shapes. The present paper concentrates on explaining in silico IEA methodology.
A potential source of confusion for all synergy calculations is that three conceptually different kinds of DERs must be considered. The three kinds are: (1) individual oneion DERs (solid lines in the example given in Fig. 1); (2) mixture baseline DERs that define absence of synergy and absence of antagonism (dashed line in the example given in Fig. 1); and (3) experimental mixture DERs, which may indicate synergy, or antagonism, or neither.
Because papers on synergy theory are sometimes overoptimistic as regards the usefulness of some specific synergy theory, drawbacks of IEA will also be discussed in the present paper. For example, many (but not all) synergy theories do not even try to predict whether mixed-agent synergy will occur; they merely try to define what synergy is (reviewed in Zaider and Rossi 1980;Berenbaum 1989;Geary 2013;Kim et al. 2015). Like SEA, IEA is among the synergy theories that have this drawback. In Online Resource 1, parts W3 and W4 include a systematic evaluation of IEA pros and cons.
There will be a number of acronyms in this paper. Table 1 lists the main ones, including less familiar ones often used here such as DER and IEA. The table also lists some of the most frequently used mathematical symbols.
To summarize, the main purpose of this paper is to explain how IEA synergy theory can be applied to accelerator experiments on murine HG tumourigenesis induced by mixed radiation fields some or all of whose beamline-entering components are one-ion HZE beams, with or without low-LET components present in the mixture.
The paper focusses on synergy theory techniques, emphasizing mathematical methods and customized computer programming more than biophysical insights. For example, all one-ion DERs for HZE ions used here will include terms that can model NTE in addition to including terms for TE. Competing HZE one-ion DERs that assume TE-only action are here omitted because the one-ion DERs that model NTE in addition to TE can illustrate synergy theory techniques adequately. No implication that joint TE-NTE one-ion DERs for this data are considered more plausible than TE-only one-ion DERs are intended.
The approach to one-ion DERs presented here for the data in Fry et al. (1985), Curtis et al. (1992), Alpen et al. (1993Alpen et al. ( , 1994 and Chang et al. (2016) has been influenced by the papers of Cucinotta and Chappell, including their seminal 2010 paper (Cucinotta and Chappell 2010). Their papers foreshadow a number of our techniques, including the use of one-ion DERs highly curvilinear at very low dose for analyzing the murine HG data.
Recent one-ion DERs for the data (Cucinotta and Chappell 2010;Cucinotta et al. 2013b;Chang et al. 2016;Cucinotta and Cacao 2017) are based on modifications of Katz' amorphous track structure approach (Katz 1988;Cucinotta et al. 1999;Goodhead 2006). The one-ion DERs used in the present paper for the same data are, as discussed above, more parsimonious than the one-ion DERs based on the amorphous track structure approach. They achieve extra parsimony by taking full advantage of a "hazard function" equation, favored by Cucinotta and coworkers and reviewed in Cucinotta and Cacao (2017). However, the present paper, due to its previously mentioned emphasis on mathematical synergy theory rather than biophysical insights, does not attempt a balanced comparison of parsimonious vs. more-parameter one-ion DERs that would take into account goodness of fit, not just parsimony.

Customized software
In the present study the open-source programming language R (Matloff 2011) is used; initially designed for statistical calculations, R has now gained wide acceptance among modelers (IEEE 2014). The customized programs developed in the course of the present study are available at https ://githu b.com/raine rsach s/mouse _tumou rs201 8. Readers can freely download, use, and modify them to evaluate the paper's conclusions critically. Detailed instructions for using the scripts are in Online Resources 1, part W1 "Guide to GitHub scripts".

One-ion DERs
Throughout this paper dose d (in cGy) is used as the main predictor variable. Fluence could have been used without changing the content or conclusions. Dose equivalents in Sv are never considered, mainly because they use photon data as a low-LET reference. That photon data and its modeling would be confounding factors, without being informative on any of the actual calculations the present paper carries out.
A mixed radiation field consists of N ≥ 2 components. Each component has a one-ion DER. DERs always model the radiogenic part of the effect-the background is subtracted out so a DER E(d) always obeys the equation E(0) = 0 by definition. Synergy theory starts with the one-ion DERs. Figure 2 illustrates the E(0) = 0 condition that holds for all DERs. The figure also illustrates "standard" properties that the one-ion DERs considered in this paper will be assumed to have unless explicitly stated to

Fig. 2
One-ion DERs-illustration of the terms "DER", "standard" one-ion DER, "strictly convex", "linear-no-threshold (LNT)" and "strictly concave" the contrary, and illustrates some convexity/concavity terminology. "Standard" one-ion DERs are defined as DERs with the following two properties.  Non-standard one-ion DERs are often useful but are not needed for this paper. An example of a DER that violates properties (1) and (2) is a linear DER with threshold. At the threshold there is a kink where the first derivative is discontinuous and the second derivative is in effect a Dirac delta function so that, roughly speaking, it is + ∞ there.
Standard one-ion DERs can have shapes more complicated than those shown in Fig. 2. Allowed, for example, are sigmoidal curves where the second derivative is positive at all small doses, is zero at one point of inflection, and is negative on up to the largest dose of interest (where, however, the slope is still required to be positive by property 1).

Scope of the data
All data used in this paper are included in open-source software freely downloadable from GitHub. The data includes GCR-simulating experimental observations on the fraction of female B6CF1/Anl mice that develop at least one HG tumour due to exposure at various doses to various one-ion beams. This fraction, the (radiogenic) tumour prevalence, is necessarily ≤ 1. It is 0 at dose = 0 because the background fraction, measured in control sham irradiation exposures, is subtracted out. Different ions differ in charge number Z, in atomic mass number u, and in ion speed β* relative to the speed of light; these parameters determine the LET L and specific kinetic energy SKE shown in Table 2. The data and data analyses are described in Fry et al. (1985), Curtis et al. (1992), Alpen et al. (1993Alpen et al. ( , 1994 and Chang et al. (2016). The Chang et al. paper concerns NSRL experiments; the other papers concern earlier experiments using the BEV-ALAC accelerator at LBL.
In Table 2, LETs for the LBL rows are for the HG in constrained mice; LETs for the NSRL rows are at the surface of the mice, which were allowed to move rather freely during irradiation. Rows 6 and 7 for 56 Fe at 600 MeV/u are an exception where data from both LBL and the NSRL were combined, as described in detail in Chang et al. (2016). Rows 1 and 2 are for low-LET protons and alpha particles. Rows 3-10 are for seven HZE ions, with 56 Fe at 600 MeV/u appearing in two different rows. The protons and alpha particles are special cases of swift light ions (SLI). In this paper "SLI" will refer to atoms with β* > 0.5 and charge number Z ≤ 3 that are almost fully ionized so that their charge can be taken to be Z proton charges. SLI with SKE ≤ 1 GeV/u, as will be the case for all SLI in this paper, are considered low LET (Norbury et al. 2016).

Omitted data
Relevant experiments at NSRL are continuing. They involve additional one-ion beams or involve, for the first time, corresponding mixed-beams obtained by rapid sequential exposure to one-ion beams. The HG tumourigenesis data from these ongoing experiments are under analysis, have not yet been published, and will not be used in this paper. However, some of the synergy theory examples presented here were selected with ongoing or planning-stage mixed-beam experiments in mind.

HZE one-ion DERs: preliminary remarks
As mentioned in "Introduction" all relevant one-ion DERs, in the literature and in the present paper, assume TE dominate at high doses. The HZE one-ion DERs used in this article assume in addition that NTE dominate at low doses. These one-ion DERs modify those tumour prevalence models reviewed in Chang et al. (2016) and Cucinotta and Cacao (2017) that assume NTE action in addition to TE action. The HZE one-ion DERs used in the present study will be standard as defined above in connection with Fig. 2. As a preview, Fig. 3 illustrates the shape of the HZE oneion DERs that result from the assumption that NTE dominate at very low doses. The calculations, including regression, that lead to these shapes will be described later in the paper. The ion is 56 Fe with Z = 26, L = 193 keV/µm, and therefore, as follows from Table 2, SKE = 600 MeV/u. In Fig. 3a it looks as if standard smoothness conditions might be violated. The slope at d = 0 looks like it might be infinite; in addition it appears as if at one point there is a kink where the first derivate is discontinuous. Figure 3b is magnified to show that actually there is no kink. Figure 3c is magnified again to show that the slope at the origin is finite. At higher doses there is concavity (despite the fact that in radiobiology high-LET TE are often linear) due merely to the fact that prevalence, because it is defined as fraction of animals at risk which have at least one HG tumour, can never go above 100%.
For the data used here all HZE experiments except one used doses less than 100 cGy. One NSRL experiment was performed using 56 Fe with SKE = 600 MeV/u at dose 160 cGy, and the theoretical curve shown in Fig. 3 provides a tolerable fit even at that high dose, as can be seen in Online Resources 1, part W2.2.
As detailed in Online Resources part W2.3, an HZE dose of 160 cGy is far larger than the total HZE ion dose, summed over all HZE components, that would be encountered on, e.g., a half-year voyage from Earth to Mars. Trying to bridge such differences, and the differences between murine HG tumourigenesis and human cancer as endpoints, by experimental design and theoretical modeling is of course a very formidable challenge. But other methods of trying to estimate astronaut cancer risk due to a prolonged voyage beyond low-earth orbit are also beset by difficulties (Cucinotta et al. 2013a). Here in the present study the focus is on the murine experiments, without trying to describe the RBE methods designed to set standards for astronaut protection.

HZE one-ion DERs: avoiding extra adjustable parameters
The starting point for the models developed in the present paper is an equation, reviewed in Cucinotta and Cacao (2017), suggested by Cucinotta and coworkers: Here E(d) is a one-ion DER. H(d) is a non-negative function, called a "hazard function" in Cucinotta and Cacao (2017), which can be chosen by biophysical modeling and then defines E(d) via Eq. (1). Short calculations show that if H(d) in Eq. (1) is chosen to have all the properties of a standard DER (as defined above in connection with Fig. 2), then (1) implies that E(d) has a positive first derivative, as required by the definition of a standard one-ion DER. The calculations that now follow will use, in Eq. (1), H(d) functions which have all the properties of standard DERs to get one-ion DERs which are standard. Equation (1) represents an important improvement over earlier models of the HG data. The equation implies that, no matter how large H(d) becomes, E(d) < 1. Thus, without the need to add any extra adjustable parameters, Eq. (1) incorporates the limitation that E(d) ≤ 1, which must hold due to the way prevalence is defined in this data set.
The HZE H(d) functions will use LET L in Table 2 as an auxiliary predictor variable. They will be taken to have additive NTE and TE contributions: Here the NTE and TE contributions are denoted by N(d) and T(d; L), respectively. The relation, Eq. (1) (3) η is an adjustable parameter, interpreted as approximately the prevalence at which NTE saturates, with saturation prevalence assumed to be the same for all HZE ions. For prevalences substantially larger than η, NTE, if any, are putatively small compared to TE, can be neglected compared to TE, or are already accounted for as cryptic modulations of observed TE.
There are no HZE data points in the region 0 < d < 3.3 cGy. Data at higher doses perhaps suggest NTE which lead to a large average positive slope, in the region 0 < d < 3.3 cGy, whose cumulative influence builds up a prevalence large enough to be detectable above background and noise at doses ≥ 3.3 cGy. To take into account NTE dose dependence in a way consistent with the strict concavity found in mechanistic models for NTE for other endpoints Brenner et al. (2001) (3) was used, with d 0 = 5 × 10 −4 cGy. Numerical explorations show that the final results of the present paper are insensitive to d 0 as long as d 0 ≪ 0.1 cGy. The factor 1 − exp( − d/d 0 ) was also used earlier when modeling NTE in this HG data (Cucinotta and Chappell 2010).
For the other term, T(d; L), in Eq.
(2) new equations were devised. After many attempts, an algebraic combination of two adjustable parameters in earlier models was found that is, for the dose range of main interest, nearly dose-independent in those earlier models. Choosing this combination as one adjustable parameter allowed the number of adjustable parameters in the one-ion HZE DER E(d; L) to be reduced from 4 to 3. "Parsimonious" models, with a minimal number of adjustable parameters in the spirit of Occam's razor, are often especially emphasized in radiobiology. Additional motivations for choosing a new HZE DER are detailed in Online Resources part W2.1.
Specifically, for H(d; L) a TE term LNT in dose was used with a coefficient involving a standard two-parametric L dependence that typically peaks at an LET of several hundred keV/µm: Combining Eqs. (1)-(4) gives comparatively simple equations for the new HZE model's DER: where The three adjustable parameters are a 1 , a 2 , and η.

HZE one-ion DERs: calibration methods and variance-covariance matrices
The background value, Y 0 , for sham-irradiated controls was taken from Chang et al. (2016) to be Y 0 = 2.7% prevalence, as estimated from all the zero-dose data, including older data not acquired at NSRL. In the present paper, Y 0 is regarded as an exact value. If one were attempting to compare a TE-only model with a model that allows for both TE and NTE, the value of Y 0 and its variance would be very important because at the very low doses where NTE putatively dominate TE, Y 0 and NTE have approximately the same magnitude. However, for reasons discussed in the last subsection of the Introduction, TE-only models were not considered in this paper.
Given Y 0 , the three adjustable parameters in the HZE one-ion DER of Eqs. (5) and (6) were obtained by inversevariance-weighted non-linear least squares regression, using all HZE non-zero-dose data globally with L as an auxiliary predictor variable. The variances for weighting were calculated by Ainsworth's formula p(1 − p)/n (Fry et al. 1985), where p is prevalence and n is the number of animals at risk. The nls() function of the programming language R determined the variance-covariance matrix during the regression calculation, and this matrix was used in subsequent 95% confidence interval (CI) estimates for IEA baseline NSNA mixture DERs. These calibration and CI estimation methods are based on assuming multivariate Gaussian distributions. Part W2.4 of Online Resources 1 gives evidence that this assumption is warranted for the data set used here.

3
After calibration, the one-ion HZE DERs were considered applicable to all ions in the Z, LET, and specific kinetic energy ranges covered by the data, i.e., applicable even to one-ion beams not in the data set. The relevant Z and L ranges were given above: 10 ≤ Z ≤ 57; approximate LET 25 ≤ L (keV/µm) ≤ 950. The specific kinetic energy range was 260 ≤ SKE (MeV/u) ≤ 1000.

One-ion DERs for low-LET proton and alpha particle beams
The one-ion DER for data on SLI, such as the HG tumourigenesis data for the low-LET proton and alpha particles was taken to be Here α is the only adjustable parameter, whose statistical distribution was determined by inverse-variance-weighted non-linear regression. Originally a linear-quadratic form for H(d) was assumed, but the quadratic term was found to be not significantly different from zero so it was omitted from the model. As was the case for HZE DERs, Eq. (1) thus facilitated use of parsimonious models.
When one-ion data for 2 < Z < 10 are added to this HG data set, it would become reasonable to treat L as having a continuous spectrum. Like other one-ion DERs used in the literature for this data (e.g., Chang et al. 2016), the DER used here in effect neglects the existence of biophysical similarities between the case L ≤ 1.6 keV/µm for Z ≤ 2 and the case L ≥ 25 keV/µm for Z ≥ 10.

Notation
Consider acute irradiation with a mixed beam of N ≥ 2 different radiation qualities. The jth radiation quality contributes some proportion r j of the total mixture dose d. The following equations, (8A)-(8D), hold In the subsequent calculations r j will always, for convenience, be independent of dose. Dose-independent proportions r j model one typical pattern for irradiation. The assumption of dose-independent proportions does not affect the final results. It implies that any one of the d j can be considered a control variable on essentially the same footing as the total mixture dose d since d j determines d, via d = d j /r j with r j > 0, and thereby determines each d i = r i d j /r j . However, a sharp distinction will be drawn between the dose control variables d and d j vs. total mixture effect considered as a control variable.
In the analyses effect size will sometimes be used to determine d and d j , instead of being determined by one of them.

SEA
Using the notations specified above, the baseline NSNA mixture DER of SEA synergy theory, denoted by S(d), is the naïve sum

Inverse functions
Inverse functions (sometimes called compositional inverse functions) play a prominent role in various synergy theories. Inverse functions are needed when using effect, rather than dose, as the independent variable. A familiar radiobiology example of inverse functions occurs when calculating the RBE of two different radiations.
The inverse of a monotonically increasing function undoes the action of the function. For example, for x > 0, √ x 2 = x so the positive square root function is the inverse of the squaring function; note that the inverse of x 2 is not x −2 . As another example exp[ln(x)] = x for x > 0, and ln[exp(y)] = y so the functions exp and ln are inverses of each other.
The IEA equation that defines absence of synergy and absence of antagonism.
When the SEA synergy theory is inappropriate, an IEA baseline NSNA mixture DER I(d) has a number of conceptual and practical advantages over other replacements for S(d) (Siranart et al. 2016). This subsection defines the elementary version of IEA synergy theory, which suffices for the present paper; an advanced version (Ham et al. 2018) is in general preferred; this advanced IEA version is described in Online Resources 1, parts W3 and W4.
Consider a mixture of N components, with each component one-ion DER "standard" as defined when discussing Fig. 2. Then each component one-ion DER increases monotonically and therefore has a compositional inverse function D j , defined for all doses of interest and thus defined for all sufficiently small non-negative effects. As discussed in subsection "Inverse Functions" above this means D j (e) = d when e = E(d). The baseline IEA no-synergy no-antagonism mixture DER I(d) is defined as the solution of the following autonomous initial value problem for a first-order, typically non-linear, ODE: In Eq. (10A), the square bracket with its subscript indicates the following calculations: first find the slope of the jth one-ion DER curve as a function of individual dose d j . Then evaluate d j using the inverse function D j with the argument of D j being the effect I already present due to the influence of all the components acting jointly. In this way different components, using the linearity of tangent lines of a curve to circumvent the curvilinearities that plague SEA synergy theory, appropriately track changes of slope both in their own DER and in the other DERs. Using d j = D j (I) in Eq. (10) instead of the seemingly more natural d j = D j (E j ) is the key assumption made. Using d j = D j (E j ) would merely lead back to the SEA baseline mixture DER S(d) (Ham et al. 2018). Parts W3 and W4 of Online Resource 1 discuss additional motivations for Eq. (10) and detail properties of its solution.
Equation (10) can be interpreted as follows. As the total mixture dose d increases slightly, every individual component dose d j has a slight proportional increase since dd j /dd = r j > 0. Therefore every mixture component contributes some incremental effect. The size of the incremental effect is determined by the state of the biological target, specifically by the total effect already contributed by all the components collectively-and not by the dose (or the effect) the individual component has already contributed. In this approach effect is considered as the control variable, and one-ion DER shapes at doses > 50 cGy may remain relevant even to a mixture with total mixture dose ≤ 50 cGy.

Calibrating background and radiogenic effects separately
Synergy is typically considered as due to interactions among agents. The present mathematical synergy analysis applies to radiogenic effects. Calibrating one-ion DERs E(d) from data always, as discussed above, used only data at non-zero doses, E(0) being 0 by definition. Background, designated by Y 0 , was based on the zero-dose data for sham-irradiated controls. Y 0 is needed when calibrating one-ion DERs from data or comparing baseline mixture DERs to data. However the main synergy calculations involve only one-ion DERs, not background plus radiogenic, effects.

Uncertainties in mixture effects
Synergy theory requires not only a way to calculate a baseline mixture DER defining NSNA but also a method of estimating uncertainties for that baseline mixture DER from mixture component one-ion DER uncertainties. Taken together these two elements constitute a default hypothesis useful for statistical significance tests on mixture observations. Without such tests, it is sometimes unclear if an unexpectedly large or small observed result does or does not call for a follow-up experiment. Here Monte-Carlo simulations (Binder 1995) were used to calculate 95% CI for I(d).
Because it is known that neglecting correlations between calibrated parameters tends to overestimate how large CI are (Hanin 2002;Ham et al. 2018) such correlations were here taken into account using variance-covariance matrices. Figure 4 summarizes the present synergy calculations in a flow chart that with minor rewording would be generally applicable in synergy theory. However, for brevity, SEA and IEA are the only synergy theories for which this paper gives any specific results.

Summary of in silico synergy theory methodology
The biophysical and statistical modeling of one-ion DERs (row a in Fig. 4) is not systematic. It uses educated guesswork trying to assemble many different kinds of biophysical, mathematical, and statistical information into equations. No results on mixed beams (row b in Fig. 4) are available in the data set of Table 2. Instead this paper emphasizes row c: mixed beams for which data will soon be available, mixed beams tentatively planned, mixed-beam experiments we believe should be carried out, and hypothetical mixed beams discussed because they illustrate some aspect of synergy theory, are analyzed.

Combinatorial complexity
The approach shown in Fig. 4 confirmed a surprising problem described in Ham et al. (2018): investigating synergy for multi-ion beams systematically is flatly impossible because there are too many possible mixtures that would need to be considered. For example, suppose a 50-50 mixture of two ions does not show synergy. A systematic approach would presumably call for two further experiments, one with, say, a 85-15% mixture of the same two radiation qualities, the other a 15-85% mixture-a total of three experiments instead of just one. The startling aspect is that if one starts with, say, a seven-ion mixture where each ion contributes 1/7 of the dose and applies the same line of reasoning one sees that a systematic approach would presumably require thousands of experiments rather than just one. Combinatorial complexity for possible dose proportions leads to an extraordinarily rapid increase for the number of potentially useful experiments as the number N of mixture components increases. Even merely trying to systematize in silico calculations becomes very difficult for N greater than about 8.

One-ion DERs
In this paper all one-ion beams for Z ≤ 2 are modeled by the same one-ion DER, namely the DER for swift light ions, Eq. (7). The equation contains one adjustable parameter, α. Calibrating α by inverse-variance-weighted non-linear regression on the non-zero-dose data gave α = 0.00153 ± 0.000124 (p < 10 −7 ). Figure 5 shows the DER, the non-zero-dose data points, and the 95% empirical CI for each data point. Since the protons and alpha particles have different LETs, namely 0.4 keV/µm and 1.6 keV/µm, respectively, plotting fluence-effect relationships would have given two separate curves.
The one-ion DERs emphasized in this paper are for high-LET ions. Calibration results for the three adjustable parameters of our HZE DERs given in Eqs. (5) and (6), obtained by inverse-variance-weighted non-linear regression, are shown in Table 3. It is seen that all three adjustable parameters differ significantly from 0, indicating parsimony.
The resulting HZE DERs are shown in Fig. W2.1 in Online Resource 1 together with the corresponding data. They are also contained in the scripts freely downloadable from GitHub. The adjustable parameter variance-covariance and correlation matrices obtained during the regression calibration are shown in Table 4.

Mixture baseline NSNA DERs
The preceding subsection describes the results obtained for one-ion DERs. It thus corresponds to row a of the synergy modeling flow chart (Fig. 4). The present subsection gives four figures describing mixture results corresponding to row c of the flow chart. The figures illustrate the IEA NSNA baseline I(d). The subsection also compares I(d) with the SEA synergy theory NSNA baseline S(d) that I(d) is designed to replace. Examples of 95% CI for I(d) are postponed until the next subsection. Figure 6 shows an example of results for a mixture of SLI and HZE beams. It is for a proton beam with specific kinetic energy SKE = 250 MeV/u and LET 0.4 keV/µm mixed with a 56 Fe ion beam with SKE = 600 MeV/u and LET 193 keV/ µm. A recent rapidly sequential NSRL irradiation whose outcome has not yet been scored used (in different proportions) the same two ions.
Comparing Fig. 6a, b, both the red and dashed black baseline NSNA mixture DER curves are lower in Fig. 6a would expect from the fact that the low-LET beam, which has the lower one-ion DER, contributes 80% of the dose in Fig. 6a but only 20% in Fig. 6b. The red IEA curve here always lies between the two one-ion DERs and thus below the green Fe curve, as one would expect from the qualitative idea that effects bigger than those estimated from analyzing one-ion DERs should be taken as indicating synergy rather than defining absence of synergy. However, it is (barely) seen that near d = 300 cGy in Fig. 6b the SEA theory's definition of no synergy is actually a little higher than the green curve.
In experiments one usually aims for a maximum HG tumour prevalence of around 25-35%, not the larger effects determined by S(d) or I(d) for doses above about 100 cGy in Fig. 6. A maximum prevalence of 25-35% is typically large enough to give a positive control for value and help give a positive control for trend without mouse mortality becoming a major confounding factor. The maximum dose of 300 cGy in Fig. 6 is substantially larger than the doses typically used in the one-ion accelerator experiments. Table W2.2 In Online Resources 1 gives, for comparisons, some details on the total doses from various radiation qualities that would be accumulated in a 6-month earth-Mars voyage which could initiate a Mars mission, assuming almost no shielding and assuming no solar energetic particle events. For example, as shown in Table W2.2, the maximum proton dose of 240 cGy in Fig. 6a is about 10 times as large as the total accumulated dose from all SLI during such a voyage. And the maximum Fe dose of 60 cGy is about four times as large as the accumulated total dose from all ionized atomic nuclei with charge number Z in the range 4 ≤ Z ≤ 28. Figure 6 also applies without any changes if the proton beam is replaced by any mixture of SLI having the same total dose as the proton beam. The reason has nothing to do with the mixmix principle. The reason is merely the fact that the low-LET DER used in the present paper is the same for all SLI.
The synergy formalism can handle mixtures with many components. Figure 7 gives an illustrative example somewhat similar to an upcoming experiment. The 20% HZE fraction of the total mixture dose in Fig. 6a is redistributed to four HZE ions. As in Fig. 6a, 80% of the total dose is contributed by low-LET protons.
It is not obvious from the figure but all five one-ion DERs are twice continuously differentiable and strictly concave. For example the apparent kink in each HZE DER at prevalence of about 5% is actually due to a very negative second derivative. 3.52 × 10 − 3 4.25 × 10 − 4 8.24 1.93 × 10 − 9 p ≪ 0.001 η (%) 5.12 × 10 − 2 1.01 × 10 − 2 5.08 1.58 × 10 − 5 p ≪ 0.001 Table 4 Adjustable parameter variance-covariance matrix and correlation matrix Parameter Variance-covariance matrix Correlations a 1 a 2 η a 1 a 2 η a 1 (µm keV −1 cGy −1 ) 1.36 × 10 −10 3.81 × 10 −9 − 5.84 × 10 −8 1 0.77 − 0.50 a 2 (µm keV −1 ) 3.81 × 10 −9 1.81 × 10 −7 − 5.59 × 10 − Mixtures of a proton beam and a 56 Fe ion beam. a Protons contribute 80% of the total mixture dose and 56 Fe ions contribute the remaining 20%; b protons contribute 20% of the total mixture dose and 56 Fe ions contribute 80%; red and dashed black lines-IEA and SEA baseline NSNA mixture DER curves, respectively; orange curves-the effect that would result if the proton beam contributed the entire mixture dose instead of only a proportion, cut off at the maximum dose the proton beam contributes; green curve-the effect that would result if the 56 Fe beam contributed the entire mixture dose instead of only a proportion, cut off at the maximum dose 56 Fe contributes; HG Harderian gland, IEA incremental effect additivity, SEA simple effect additivity, NSNA no synergy and no antagonism, DER dose-effect relation 1 3 In Fig. 7, the red baseline IEA curve I(d) is a reasonable definition of NSNA for the mixture, nested between the various component DERs. Intuitively speaking, it is here below all the HZE DERs because so much of the total dose is contributed by low-LET radiation (orange curve). However, for d < 30 cGy SEA theory makes the absurd claim that some mixture effects larger than any component radiation could produce on its own are not evidence for synergy but rather define absence of synergy. This inability of SEA to handle the curvilinearities shown is typical. S(d) gives unrealistic over-estimates when all components of a mixture have strictly concave DERs. Figure 7 also illustrates the comments in the Method section about allowing the auxiliary LET predictor variable to act as a numerical variable rather than a categorical variable. There is no data in our data set for the specific LETs assigned to the four HZE components of the Fig. 7 mixture. However, the HZE DER of Eqs. (5) and (6), Eq. (9) which defines S(d), and the IEA Eq. (10) all make sense if L is read as a numerical rather than a categorial variable, so it was possible to carry out all the calculations needed to get the SEA and IEA curves in Fig. 7.
Finally, Fig. 7, like Fig. 6, is applicable without any changes if the proton contribution is replaced by any SLI mixture contributing the same dose added up over all the SLI beams. As before, this flexibility is not related to the mixmix principle. Figure 8 shows synergy theory results for an NSRL experiment using a 50-50 mixture of two HZE components. Scoring the experimental results began in May 2018 and is ongoing. The two components (silicon and iron ions) are characterized in Table 2.
It is seen in Fig. 8 that the SEA NSNA function S(d) is, as in Fig. 7, unacceptably high: the mixture SEA NTE saturates at about the sum of the two saturation heights instead of saturating at about the height of each one.
There are many further examples which illustrate various aspects of IEA synergy theory. The final example in this subsection, Fig. 9, is a hypothetical mixture of all seven HZE ions in Table 2 (Ne, Si, Ti, Fe at two different energies, Nb, and La), each contributing 1/7 of the total mixture dose.
It is seen in Fig. 9 that SEA synergy theory gives a particularly unreasonable baseline here, essentially because there are many components all of whose DERs are strictly concave, with very strong concavity at doses just above zero.

Uncertainties in baseline mixture DERs: accounting for parameter correlations
The following two final figures focus on uncertainties in the IEA NSNA baseline DER of a mixture. Figure 10 adds uncertainty information to Fig. 8.
Because it used sampling from the variance-covariance matrix shown in Table 4, calculating the yellow ribbon in Fig. 10 properly took into account the adjustable parameter correlations shown in Table 4. Especially in view of the emphasis NASA risk limitation policies for astronauts place on 95% CI (reviewed, e.g., in Cucinotta and Chappell 2010;Cucinotta and Cacao 2017), taking parameter correlations into account properly is important. Figure 11 compares, for the hypothetical mixture of seven HZE used in Fig. 9, the appropriate 95% CI with the broader 95% CI that would be obtained by a failure to take adjustable parameter correlations into account. Figure 11 also illustrates the feasibility of the Monte-Carlo simulations for mixtures having many components.
In Fig. 11 it is seen that above 25 cGy the difference between top and bottom for the uncorrelated blue ribbon is more than twice as large as the difference for the yellow ribbon. For example, at 49 cGy the respective 95% CIs are 15.6% vs. 6.0% prevalence. The hypothetical point would indicate statistically significant synergy for this one particular cocktail of HZE, since there is no overlap of the error bars with the yellow ribbon. Neglecting correlations would have led to the incorrect interpretation that the point is suggestive of synergy but does not indicate statistically significant synergy.

Comments on the results and figures
To illustrate IEA synergy theory, results exemplifying the following were given: the importance of having high-quality one-ion DERs available; the use of IEA as an alternative to  Table 2; b magnified on the low-dose region of a; red and dashed black curve-IEA and SEA NSNA baseline DER curves, respectively; seven shorter curveseffect that one ion would have if it contributed the entire mixture dose d instead of contributing only d/7, cut off at d/7, with (solid black at the top, green, purple, almost hidden dark blue, orange, light blue, lilac at the bottom) curves corresponding to the ion in Table 2 which has LET L (in keV/µm) of (250,193,464,100,70,953,25) SEA when one-ion DER curvilinearity requires an alternative; calculation of baseline mixture DER 95% CI taking into account correlations between one-ion DER adjustable parameters so that much tighter CI are properly identified; and the extraordinarily rapid increase in the number of mixtures needed to look for possible synergy in N-component mixtures as N increases. Figures 7, 8 and 9 helped emphasize that taking SEA synergy theory for granted is a mistake. The absurdly large values of S(d) at very small doses violate not only intuitive ideas of what synergy means but also biophysical reasoning. Such large values correspond to the picture that during acute irradiation low-dose NTE in a mixture saturate at the sum of the components' NTE. But biophysical reasoning would suggest for acute doses something like the average instead of the sum. Consider two radiation qualities, say A and B. One-ion experiments strongly suggest that once radiation quality A causes cell #1 in a multi-cell interaction region corresponding to the NTE signal range to send out an NTE signal saturation occurs: a second hit on cell #2 in the same interaction region by radiation A causing cell #2 to send out additional bystander signals does not further enhance NTE noticeably. To suppose that if it is instead radiation quality B that hits cell #2, causes a bystander signal that would be almost ignored if the hit were by radiation quality A, and thereby doubles the NTE is far-fetched. How are the bystander cells, receiving signals from cells #1 and #2, supposed to know that the radiation quality of the radiation that directly affected cell # 2 was not the same as the radiation quality that directly affected cell # 1? On this argument, the SEA curves in Figs. 7, 8 and 9 should be rejected, and the figure illustrates, once again, that relying on SEA synergy theory can be misleading-a fact that has been known but too often been ignored in many subfields of biology during the last 150 years.
Using a biologically motivated microdosimetry modeling approach more mechanistic but less general than our mathematically motivated IEA approach Drs. D. Brenner and I. Shuryak (personal communication) obtained results on NTE saturation fully consistent with the present IEA results. They found that NTE saturation for a mixture each of whose components' NTE saturates at the same level should be expected to occur at about that level. They also found that then SEA overpredicts the mixture saturation level much as shown in Figs. 7, 8 and 9. IEA calculations generalize these results mathematically; for example the IEA baseline for a mixture of many NTE components which saturate at different levels, saturates at about the largest one-component NTE saturation level. Figures 8, 9, 10 and 11 concern HZE-only mixtures, with none of the total mixture dose entering the beam contributed by low-LET ions. Such HZE-only mixtures could help determine whether or not the extraordinarily high RBEs of HZEs in the data set used here are sometimes exacerbated by synergy for some specific HZE cocktail that happens to produce a particularly dangerous ionization pattern. Subsection "Choosing GCR-simulating mixtures and dosing protocols for accelerator experiments" below discusses implications of this possibility for deciding which mixed beams to use in NSRL experiments. Figures 10 and 11, on 95% CI for I(d), emphasize the importance of taking correlations into account. This is feasible and important not only for the approach here but also in most situations where adjustable parameter calibration via regression is used (Hanin 2002), so it should become routine.

Synergy theory in radiobiology
For the foreseeable future radiobiologists studying mixed radiation field effects will almost inevitably emphasize possible synergy and antagonism among the different radiation qualities in the mixture. Therefore, trying to find a systematic quantification of synergy, general enough to cover most cases of radiobiological interest and precise enough to enable credible estimates of statistical significance, is worthwhile. But it is not easy.
One main problem is the following: the common belief that synergy can always be defined as an effect greater than the SEA baseline mixture DER S(d) is wrong. In some important cases S(d) is clearly inappropriate. As detailed in Online Resources 1, parts W3 and W4 there are many published alternatives, among which the IEA synergy theory seems optimal. In the present paper IEA synergy theory was illustrated with detailed examples of mixtures of ions in the GCR spectrum, since more experimental information on such mixtures will soon become available. A caveat is that all known synergy theories, including IEA, have substantial limitations. Online Resources 1, subsection W3.2.1 lists many of these.

Choosing GCR-simulating mixtures and dosing protocols for accelerator experiments
Many recent papers, e.g. Norbury et al. (2016) and Straume et al. (2017) emphasize analyzing GCR-simulating mixtures where a large majority of the dose (in Gy) is contributed by a cocktail of low-LET protons and alpha particles. In contrast, Figs. 8, 9, 10 and 11 of the present paper concern mixtures that involve only HZE, no swift light ions at all. Moreover, the present paper discussed dosing protocols specifically tailored to the murine HG tumourigenesis endpoint by choosing doses high enough to produce prevalences that are well distinguishable from noise but low enough to avoid confounding mouse mortality effects. In contrast recent NSRL experiments are moving in the direction of standardizing, by choosing identical dosing protocols for many different biological endpoints. This subsection considers those two rather drastic discrepancies: all-HZE mixtures vs. mixtures with only small HZE dose proportions; standardized dose protocols vs. endpoint-dependent dose protocols.
As detailed in the present paper's "Introduction" above, the very high RBEs for HZE ions in murine HG tumourigenesis experiments have long been a concern for NASA. The HZE ions are less well understood than low-LET radiations and, for a number of complex biological endpoints, are very much more effective or even produce qualitatively different effects than low-LET radiation.
The Norbury et al. (2016) paper, which had a broad variety of co-authors including biologists and physicists, addressed the problem of what mixtures and dosing protocols should be used at NSRL. Emphasis was on seeking a specific mixed field whose composition would be representative of interplanetary space conditions. It was suggested that such a field, and a fixed dosing protocol, could be used as a standardized NSRL procedure for many different groups investigating various questions. However, it is important to note that the report also pointed out the need for flexibility. It argued that using mixtures and dosing protocols clearly not representative of conditions in outer space should sometimes be used, to allow full use of information from past one-ion experiments when studying synergy, and for a variety of other purposes.
That point about synergy seems well taken. The hypothesis that statistically significant synergy among HZE-only mixed-beam components does not occur is reasonable, since detailed theoretical analyses of track structure have not uncovered any clear mechanism whereby overlapping tracks from different ions would somehow cooperate to produce a particularly dangerous overall ionization-excitation pattern. But finding experimental evidence for or against the hypothesis should arguably be a high-priority project in view of the large RBEs found in the one-ion experiments. HZE-only mixed beams are the most appropriate beams for that task since in such mixed beams confounding effects of swift light ions are not present.
Speaking more generally, the comparatively flexible approach in Norbury et al. (2016) toward NSRL experiment guidelines seems to deserve more discussion than it has received. Perhaps the more flexibility the better. There may be no one mixed field truly representative of a majority of the many varied situations that will occur during interplanetary voyages. Different biological endpoints such as chromosome aberrations vs. murine tumourigenesis have different radio-sensitivities and different relaxation times, so that no one mixed-field dosing protocol can be appropriate for more than a minority of the different endpoints of interest. The protocols analyzed in the present paper included some, such as the protocol summarized in Fig. 11, whose experimental implementation at NSRL could be quite informative about synergy but which by no stretch of the imagination could be considered representative of any deep space mixed-field. Standardization neglects the factor of biological variability, typically dominant for any biological endpoint though not present for physics endpoints such as comparing responses of nominally identical semiconductor devices at risk for radiation damage. Flexibility would allow giving biological variability its due importance.

Conclusions
Synergy theory will continue to be used to plan experiments involving mixed radiation fields and interpret the results of such experiments. It can and should include calculations that give 95% confidence intervals based on sampling variance-covariance matrices. If NTE are important, SEA NSNA baseline mixture DERs should be ignored or used only cautiously. Probably IEA synergy theory is the optimal choice in radiobiology when a replacement for SEA synergy theory is needed. It may be optimal in other fields as well.
In any case, all synergy theories have more limitations than is generally realized.
Whether mixed GCR fields sometimes produce statistically significant synergy for human carcinogenesis or other deleterious effects is not known. Upcoming mixture experiments will hopefully provide clues to the answer.