Using thermodynamic parameters to calibrate a mechanistic dose-response for infection of a host by a virus.

Assessing the risk of infection from emerging viruses or of existing viruses jumping the species barrier into novel hosts is limited by the lack of dose response data. The initial stages of the infection of a host by a virus involve a series of specific contact interactions between molecules in the host and on the virus surface. The strength of the interaction is quantified in the literature by the dissociation constant (Kd) which is determined experimentally and is specific for a given virus molecule/host molecule combination. Here, two stages of the initial infection process of host intestinal cells are modelled, namely escape of the virus in the oral challenge dose from the innate host defenses (e.g. mucin proteins in mucus) and the subsequent binding of any surviving virus to receptor molecules on the surface of the host epithelial cells. The strength of virus binding to host cells and to mucins may be quantified by the association constants, Ka and Kmucin, respectively. Here, a mechanistic dose-response model for the probability of infection of a host by a given virus dose is constructed using Ka and Kmucin which may be derived from published Kd values taking into account the number of specific molecular interactions. It is shown that the effectiveness of the mucus barrier is determined not only by the amount of mucin but also by the magnitude of Kmucin. At very high Kmucin values, slight excesses of mucin over virus are sufficient to remove all the virus according to the model. At lower Kmucin values, high numbers of virus may escape even with large excesses of mucin. The output from the mechanistic model is the probability (p1) of infection by a single virion which is the parameter used in conventional dose-response models to predict the risk of infection of the host from the ingested dose. It is shown here how differences in Ka (due to molecular differences in an emerging virus strain or new host) affect p1, and how these differences in Ka may be quantified in terms of two thermodynamic parameters, namely enthalpy and entropy. This provides the theoretical link between sequencing data and risk of infection. Lack of data on entropy is a limitation at present and may also affect our interpretation of Kd in terms of infectivity. It is concluded that thermodynamic approaches have a major contribution to make in developing dose-response models for emerging viruses.

Introduction

Microbiological risk assessment (MRA) requires a dose-response relationship to translate the exposure (i.e. number of pathogen particles entering the host through a given route) into the probability of infection. Infection by an oral pathogen is defined as the multiplication of organisms within the host, followed by excretion (Haas et al., 1999) and, for the purpose of the work here does not include progression of disease or the host acquired immune response. Obtaining dose-response data for humans has generally relied on volunteer challenge experiments e.g. Cryptosporidium parvum in students (Okhuysen et al., 1998) or using outbreak data to back-calculate the relationship between measured exposures and infection rates (Teunis et al., 2004). There are limitations to both approaches particularly with emerging pathogens for which the exposure routes may not be fully elucidated, and for pathogens with serious clinical outcomes, e.g. Zaire ebolavirus (EBOV). Furthermore zoonotic viruses emerge through jumping the species barrier from an animal source to humans, e.g. Nipah virus (NiV) and EBOV, and in this respect the dose response would be for a one-off event that may be inefficient and difficult to reproduce without large numbers of animals. An additional complication is that the pathogen may adapt to the new host, such that its infectivity increases. This is well established for filoviruses in laboratory animals where the infectivity per plaque-forming unit (pfu) may change by several orders of magnitude with passaging (Gale et al., 2016), and has recently been demonstrated for EBOV Makona adapting to humans through an amino acid substitution in its glycoprotein during the recent catastrophic outbreak in West Africa (Diehl et al., 2016, Urbanowicz et al., 2016). That outbreak also raised many questions regarding the unknown potential for companion animals (cats and dogs) to serve either as a reservoir or vector for the virus and so be involved in transmission of EBOV to humans and other animals. The absence of dose-response data for EBOV in humans limits development of MRAs for the risk of infection of citizens in the EU for example from EBOV in illegally imported bushmeat. Indeed, it has been proposed that the infectivity to humans of an EBOV pfu may differ not only from bushmeat samples from different wildlife species (e.g. fruits bats and nonhuman primates) but also from different individuals of the same species depending on the degree of host adaptation (Gale et al., 2016). In effect no two pieces of bushmeat from EBOV-infected wildlife may be the same in terms of infectivity to humans, although this remains to be proved. There is clearly a need for novel approaches to calibrate dose-response relationships for the purposes of MRA for emerging pathogens.

The infection process of a host cell can be broken down into the component steps and modelled mathematically (Handel et al., 2014) and the probability of infection can be expressed as a function of the combined probabilities of each step (Gale et al., 2014). These steps include overcoming the initial host defenses, binding of the virion to its host cell receptor, entry to the host cell (i.e. internalisation and uncoating of the virion), and replication, capsid assembly and budding (Gale et al., 2014). Previously it was demonstrated that a dose-response model could, in part, be parameterized using thermodynamic data for some of the key molecular interactions in the infection process (Gale, 2017). The beauty of thermodynamic data is that they can be measured experimentally by biochemists (in some cases just using molecular components e.g. cloned virus protein and host receptor protein (Wang et al., 2016)) and do not involve live animal or human volunteer studies, which is a major advantage for dangerous pathogens. Furthermore the effect of amino acid substitutions in the host receptors on binding affinity can be measured directly (Yuan et al., 2015). The possibility of applying thermodynamics is further developed here for two of the key steps in the infection process of a host by a virus. The first step modelled is the probability of the virus overcoming the innate host defenses posed by mucin protein molecules and the pathogen recognition receptors (PRRs) produced by the host. Mucins have sugar units on their surface which bind to components on the surface of the virus, for example the haemagglutinnin (HA) glycoprotein molecules of influenza virus (de Graaf and Fouchier, 2014) or the VP1 of norovirus (NoV) (de Rougemont et al., 2011). Mucus present in the respiratory tract hampers influenza virus infection and in the case of humans predominantly contains  α2,3-sialic acid receptors. Indeed influenza viruses with  α2,3 specificity were inhibited by human mucins (de Graaf and Fouchier, 2014). The PRRs include the mannose binding protein (MBP) which has carbohydrate-recognition domains (CRD) which bind to regularly repeating sugar units on pathogen surface (Taylor and Drickamer, 2006). The second step modelled here is the binding of the virus to its specific receptors on the host cell surface. The approach here is developed for a generic faecal/oral virus such as NoV and rotavirus which infects epithelial cells lining the intestine (Boshuizen et al., 2005, De Rougemont et al., 2011), but could be applied to influenza A viruses which are inhaled and infect cells of the trachea and lung (de Graaf and Fouchier, 2014).

This paper first gives an overview of a mechanistic dose response model to introduce two probability parameters, namely the fraction, Fv, of virus escaping the mucin defense barrier and the fraction, Fc, of host cells with bound virus. The Methods section sets out a difference equation method to model Fv and Fc as a function of the mucin: virus ratio and virus dose in the intestine, respectively. Central to determining Fv and Fc is the strength of binding of the virus to the mucin and host cell as defined by the equilibrium constants Kmucin and Ka respectively. In the Theory section, the application of published data on the binding of the virus surface envelop glycoprotein (GP) to host cell receptor (Cr) molecules or to mucin molecules is reviewed in terms of determining Kmucin and Ka in order to parameterize the dose-response. Particular reference is made to using the dissociation constant Kd which is routinely determined experimentally for virus GPs binding to Cr molecules (Gambaryan et al., 2005, Raman et al., 2014; Yuan et al., 2016). It is then shown how the strength of virus/host cell binding (i.e. the magnitude of Ka) may be predicted from changes in two thermodynamic parameters, namely enthalpy (H) and entropy (S). The effects of amino acid changes at the contact surfaces of the virus GP and Cr on the enthalpy are considered with a view to the future parameterization of dose-response models based on genetic sequencing data. Entropy changes are also considered both in terms of virus binding and also in the interpretation of Kd data.

Methods

Overview of the development of a mechanistic dose-response model for infection in the intestine

The model parameters and variables are summarised in Table 1 . On ingestion of the initial virus challenge dose, Vinitial, by the host there are a number of immediate host defences in the mouth and gastrointestinal (GI) tract including the mucus barrier, decoy receptors and the innate immune system that selectively bind and hence remove the virus (McGuckin et al., 2011). For example the histoblood group antigens (HBGAs) are genetically determined glycans to which NoV selectively binds and are present on both decoy receptors in the saliva and on mucin, the main protein component of mucus (Shanker et al., 2011). The total number of viruses surviving the mucin barrier, and getting through to the intestine is given bywhere Fv is the fraction of free virus, i.e. that not bound to mucin. As shown in Fig. 1 , Fv can be modelled by two parameters, namely the total number, Muctotal, of mucin molecules in the mucus in the saliva and GI tract and an association constant, Kmucin (defined below) that quantifies the strength of binding of the virus to a mucin molecule. Thus by inserting Fv from Fig. 1 for a given mucin concentration into Eq. (1), the total number of free virus particles in the intestine and available to initiate infection of the epithelium may be modelled. On reaching the intestinal epithelium, a free virus particle binds to the surface of a host cell. The probability of infection of the host, phost, equals the probability of successful infection of at least one cell and is related to the number of cells (C.V) with bound virus by:-where pcell is the probability of successful infection of a host cell given a virus has bound to its surface. Thus the more cells with bound virus then the greater the chance that infection will be successful in at least one of them. The probability pcell depends on ability of the bound virus to enter the cell, replicate and bud (Gale et al., 2014, Gale, 2017) and is not discussed further here. Nowwhere Fc is the fraction of cells with bound virus, and Ctotal is the number of cells in the host intestinal epithelium. As shown in Fig. 2 , Fc is directly proportional to the total number of virus particles in the intestine, Vintestine, and is also dependent on the strength of the binding interaction between the virus and the epithelial cell surface as quantified by the association constant, Ka, which is now defined.

Table 1

Summary of model parameters.

Parameter	Description	Comments
phost	Probability host organism is infected	Overall objective of dose-response
p1	Probability of infection from ingestion of a single virion by the host	Parameter to be obtained from mechanistic dose response approach developed here for direct use in conventional dose response model (Eq. (16))
pcell	Probability that a host cell becomes infected given virus has bound to its surface	Not discussed further here
Vinitial	Challenge dose of virus to the host organism	Oral exposure in MRA. [Vinitial] is the concentration of total virus in the simulated intestine at 2.19 × 10−15 M
Vintestine	Total number of viruses not bound to mucin, and getting through to the intestine to initiate infection of host cell.	Within intestine, Vintestine includes both virus bound to host cells and not bound to host cells, i.e Vintestine = C.V + Vfree.
V.Muc	Number of viruses bound to mucin
Mucfree	Number of mucin molecules with no bound virus
Muctotal	Total number of mucin molecules in the host mucus in the host saliva and GI tract	This is varied relative to the fixed virus challenge dose (Vinitial) in Fig. 1
Fv	Fraction of virus that is not bound to mucin, i.e. survives to infect host cells	Probability virus breaks through mucin barrier into intestine
Fc	Fraction of host cells (Ctotal) with bound virus	Probability that cell has virus bound to it
Ctotal	Number of cells in the host intestinal epithelium	In the simulations, Ctotal is constant at 4.15 × 108. Concentration of total host cells [Ctotal] is 2.19 × 10−15 M.
C.V	Number of host cells with bound virus	[C.V] is concentration (M) of cells with bound virus.
Cfree	Number of host cells without bound virus, i.e. free	[Cfree] is concentration (M) of free cells (i.e. cell with no virus attached)
Vfree	Number of viruses in the intestine not bound to cells, i.e. free	[Vfree] is concentration (M) of free virus (i.e. not bound to cells)
Ka/Kmucin	Association constants for binding of virus to host cell and mucin respectively	Units M−1. Related to reciprocal of Kd (Eq. (17))
Kd	Dissociation constant measured experimentally between individual virus molecules and host molecules.	Units M. The smaller Kd in magnitude, the stronger the binding.

Fig. 1

Fraction, Fv, of virus not bound to mucin plotted as a function of the total mucin: total virus ratio. The virus challenge dose in the 0.314 dm3 volume of intestine was fixed at 4.15 × 108 virus particles and the number of mucin molecules was increased from 103 to 4.15 × 1012 molecules as represented by the symbols. Binding of all the virus is achieved below the horizontal dotted line which represents 1 unbound virus remaining in the intestine. Solid lines with symbols represent points calculated with difference equation approach (see text) with Kmucin values of 1022 (x), 1020 (●), 1018 (▲), 1015 (■), 1013 (Δ), 1011.7 (♦) and 109 (□) (M−1). Dashed lines represent Eq. 11 assuming [Muctotal] ∼ [Mucfree] with Kmucin values from left to right of 1020, 1018, 1015, 1013, 1011.7, 109, 107, and 105 (M−1). Eq. (11) fails for high Kmucin values at mucin: virus ratios of < 1:1 (arrow).

Fig. 2

Fraction of host cells with bound virus plotted as a function of virus dose in intestine for pathogen/host cell interaction of increasing binding affinity as represented by Ka values of 105 (o), 107 (□), 109 (Δ), 1011.7 (♦), 1013 (■), 1015 (▲), 1018 and 1020 (●) (M−1). Note the lines for 1018 and 1020 are superimposed because the affinity is so high that either all the virus is bound at low virus doses or all the host cells are saturated at high virus doses. Vertical dotted line represents virus: host cell ratio of 1:1.

Expressing the fraction, F

Each binding process is represented by a dynamic equilibrium. Thus, within the given volume of the intestine, free virus (Vfree) and host cells with no bound virus (Cfree) are in dynamic equilibrium with cells with bound virus (C.V) as represented by:-

The association constant, Ka, is expressed in terms of the concentrations (Gale, 2017) as:-where the square brackets, [], represent the concentration in moles dm−3 (M). The term “in dynamic equilibrium” means the process is reversible (Handel et al., 2014) such that free virus binds to free host cells to form C.V with a rate, kon, and the C.V complexes then dissociate into Cfree and Vfree at a slower rate, koff, depending on the strength of binding. Thus the association constant, Ka (and similarly Kmucin) may also be written in terms of the association/dissociation rates as:-

Visualising Ka (or Kmucin) in terms of on/off rates may be easier conceptually and has implications when host or virus factors selectively change kon or koff (see below) thus affecting Ka (or Kmucin) according to Eq. (6). As shown by Gale (2017), the fraction, Fc, of host cells with bound virus is given by:-

Substituting [C.V] with Eq. (5) and rearranging gives Fc in terms of the free virus concentration [Vfree] and Ka:-

Similarly, when the virus enters the host there will be a dynamic equilibrium between virus that is not bound to mucin and hence reaches the intestinal epithelium (Vintestine), free mucin (Mucfree) and virus bound to mucin (V.Muc).

The strength of the binding is reflected by the association constant, Kmucin, between virus and mucin and is expressed as:-

The fraction of free virus (Fv) is given by:-

Replacing [V.Muc] with Eq. (9) and rearranging gives Fv in terms of the free mucin concentration and Kmucin:-

Modelling how Fc varies with dose of virus in the intestine

The objective is to construct a plot of Fc as a function of Vintestine as shown in Fig. 2. Eq. (7) may be used with [Vfree] representing [Vintestine] when the total number of virus particles, Vintestine, greatly exceeds the number of host cells (Ctotal) such that [Vfree] is relatively unaffected by virus binding and therefore [Vfree] approximates [Vintestine]. This is also acceptable for low binding affinity viruses as in the species barrier model of Gale (2017) such that very little virus is bound even at high virus doses. However, in many natural infection processes (e.g. through drinking water), the host may be challenged by very low numbers of pathogen which bind with high affinity to host cells such that [Vfree] is greatly diminished compared to [Vintestine] and tends to zero as all the virus is bound. The problem is that Eq. (7) cannot be expressed mathematically in terms of [Vintestine]. The solution adopted here to model Fc in terms of [Vintestine] comprises three steps. The first step involves setting up a model for a host intestine so that the concentrations of virus and host cells may be defined in order to calculate Ka using Eq. (5). In the second step, a difference equation approach is used to produce a range of [Vfree], [Cfree] and [C.V] combinations at six different total virus (Vintestine) doses. From these concentrations, Ka values and Fc values are calculated and are plotted in Fig. 3 . In the third step, for each of the six Vintestine doses, the Fc is read off Fig. 3 for given Ka values. Fc is then plotted against Vintestine for each Ka in Fig. 2.

Fig. 3

Fraction of host cells with bound virus plotted as a function of binding affinity as represented by the association constant Ka. Doses of 1.0 × 103 (♦), 1.0 × 106 (■), 4.15 × 108 (o), 4.15 × 109 (Δ), 4.15 × 1010 (▲) and 4.15 × 1011 (●) viruses in a simulated intestine with 4.15 × 108 susceptible host cells.

The concentration of host cells in the model intestine

Developing a dose-response model based on Ka and Kmucin requires an estimation of the volume in dm3 (i.e. litres) within which the challenge dose of pathogen has access to susceptible cells in the host. This is needed to convert particle numbers in moles into concentrations for use in Eq. (5). The volume into which the pathogen enters within the host is for example the lumen of the intestine or even a drop of body fluid (e.g. blood) on a mucosal surface in the case of EBOV. For the purpose of the model, a 100 cm length of host intestine with a radius of 1 cm is simulated. The total surface area is therefore 628.3 cm2 and the volume inside the intestinal lumen is 0.314 dm3 (1 dm3 = 1000 cm3). According to Rosen and Misfeldt (1980) the density of cells in dog kidney epithelia is 6.6 × 105 cells per cm2. The total number, Ctotal, of epithelial cells over the 628.3 cm2 surface of the model intestine is thus 4.15 × 108, and the density of cells by volume is 1.32 × 109 per dm3. Dividing by the Avogadro number, L, of 6.02 × 1023 particles per mole (Price and Dwek, 1979) gives a mole concentration by volume [Ctotal] of total host cells of 2.19 × 10−15 M. It is reassuring to note that this is similar to the cell concentration of 4.1 × 10−15 M reported for the canine kidney cells in the avian influenza virus binding experiments of Nunes-Correia et al. (1999). Handel et al. (2014) in their simulation for influenza virus infection use a concentration which is ∼100-fold higher although they simply assumed a packed collection of cells each with a volume of 20 µm3. The model here therefore takes into account the large void volume of the intestine within which a faecal-oral pathogen is diluted.

Virus doses used in the simulation

Simulations were undertaken for six virus doses (Vintestine) from 1.0 × 103 to 4.15 × 1011 virions in the model intestine (Fig. 2), representing virus: host cell ratios from 2.4 × 10−6:1 to 1000:1 and concentrations from 5.3 × 10−21 M to 2.2 × 10−12 M.

Difference equation approach to model Fc against Ka for each virus dose (Fig. 3)

For low virus doses (i.e. Vintestine < Ctotal), Kas were calculated using [C.V], [Vfree] and [Cfree] in Eq. (5) over the full range of virus binding i.e. from one bound virus to all Vintestine viruses being bound. According to Eq. (4), each cell can only bind one virus, and the number of cells with bound virus therefore equals the number of viruses bound to cells which was calculated as:-over the range of Vfree from 0 (all viruses bound) to Vintestine (no viruses bound). For high virus doses (i.e. Vintestine > Ctotal) the simulation was run for Vfree from 0 to Ctotal (such that C.V is always positive) because the model in Eq. (4) assumes each cell can only bind one virus. For each C.V, the number of free cells, Cfree, was calculated as:-

Values for [Cfree], [Vfree] and [C.V] were calculated as Cfree, Vfree and C.V respectively, divided by the volume of the intestine (in dm3) and L and then used to calculate Ka in Eq. (5) over the range of C.Vs. For each virus challenge dose, values of Fc are plotted as a function of Ka in Fig. 3 where Fc is calculated as:-and Ctotal, is constant at 4.15 × 108 cells (Table 1).

Modelling how FV varies with ratio of mucin to virus: surviving the mucin defence

The fraction of free virus, Fv, represents the probability that virus is not bound to the mucin and in effect escapes the mucin barrier. The objective is to model Fv as a function of the mucin to virus ratio over a range of Kmucin values as shown in Fig. 1. The simple approach to model Fv is to use Eq. (11) with [Mucfree] representing the total concentration of mucin in the intestine, [Muctotal], which can be measured in a host experimentally. As discussed below, this fails at mucin: virus ratios of <1:1. Therefore a difference equation approach was used. For this, the challenge dose, Vinitial, in the host intestine was fixed at 4.15 × 108 virions and the fraction of virus not bound to mucin, Fv was then calculated for ten values of Muctotal ranging from 103 mucin molecules to 4.15 × 1012 mucin molecules, representing mucin: virus ratios ranging from 2.4 × 10−6:1 to 10,000:1. This was done using the same stages as described above for virus binding to host cells for seven Kmucin values ranging from 109 to 1022 M−1. Thus for the number of free mucin molecules, Mucfree, ranging from 0 (i.e. all mucin molecules bound to virus) to Muctotal (i.e. all mucin molecules free of virus), values of V.Muc and Vintestine were calculated as:

By dividing by the volume of the intestine and L, the values of V.Muc, Vintestine and Mucfree were converted to corresponding concentrations, namely [V.Muc], [Vintestine] and [Mucfree], from which Kmucin values were calculated using Eq. (9). For each Kmucin value, Fv was calculated as:and a plot (not shown) of Fv versus Kmucin constructed in the same was as for virus binding to host cells in Fig. 3. From that plot, values of Fv for each of the ten Muctotals were read off for a given Kmucin and plotted against the mucin molecule:virus ratio (i.e. Muctotal:Vinitial) in Fig. 1.

Effects of stochasticity

Stochasticity in the challenge dose (Vinitial) would be addressed in the exposure calculation and is outside the scope of this work. The mechanistic dose-response model developed here using the difference equation approach is not affected by stochasticity because the values of C.V and Vintestine calculated in Eqs. (12) and (14) respectively are “given” integers, and thus Fc and Fv calculated by Eqs. (13) and (15) respectively are exact for each integer C.V and Vintestine (Ctotal and Vinitial being constant integers in the simulation). The probability of infection of the host according to Eq. (2) is calculated for each integer value of C.V from 0, 1, 2, … to Ctotal or Vintestine (depending on which is lower) and is not affected by stochasticity. The effect of stochasticity is considered here for the use of Eq. (7) to calculate Fc and where C.V is a fraction.

Spatial heterogeneity in the model of the intestine

The model assumes that the free pathogens and mucins are homogeneously distributed within the lumen of the intestine, such that the concentrations at equilibrium are constant along the 1 m length of the simulated intestine. The spatial heterogeneity of the pathogen in the lumen is not known and could vary depending on the pathogen distribution in the food, water or faecal/vomit contamination ingested by the host. For example, ingestion of a small amount of faeces laden with virus could give much higher virus concentrations at certain parts of the lumen, depending on the degree of mixing within the lumen. In contrast ingestion of a 0.314 dm3 volume of water contaminated with pathogen could give a homogeneous (Poisson) distribution along the 1 m length. This is not considered further here. Spatial heterogeneity will exist in the mucin concentration because the mucus forms a layer lining the intestine wall (McGuckin et al., 2011). This needs consideration in the further development of this mechanistic approach.

Results

Here an intestine is simulated with an internal volume of 0.314 dm3 and a total of 4.15 × 108 susceptible cells in the intestinal epithelium. In Fig. 1, an oral challenge dose of 4.15 × 108 viruses is administered, and the fraction of these escaping the mucin barrier and reaching the epithelium is modelled. Fig. 3 shows the fraction of epithelium cells with bound virus as a function of Ka for six doses of virus that have got through the mucin barrier ranging from 1000 virions to 4.15 × 1011 virions. In Fig. 2, the points from Fig. 3 are replotted in the form of a dose-response which relates the fraction of host cells with bound virus to the dose of virus that has got through the mucin barrier. These dose-response relationships are presented for a range of Ka values from 105 to 1020 M−1 in Fig. 2.

Assessing the magnitude of Kmucin needed for an effective mucus barrier

The fraction, Fv, of virus escaping the mucus barrier is plotted as function of the mucin: virus ratio for a range of Kmucin values from 109 to 1022 (M−1) in Fig. 1. The horizontal dotted line in Fig. 1 represents just one free virus remaining in the 0.314 dm3 volume of the simulated intestine. At values of Fv below this line there is < 1 free virus in the simulated intestine and in effect all of the 4.15 × 108 virions in the initial challenge dose (Vinitial) are bound to mucin. Thus, the effectiveness of the mucin barrier can be assessed simply in terms of Kmucin and the mucin: virus ratio required to bring Fv to below this line in Fig. 1. At the very high Kmucin value of 1022 M−1 all of the virus is bound as the mucin: virus ratio exceeds 1:1. However, at progressively lower Kmucin values, less and less of the virus is bound at a given mucin: virus ratio. Thus even at Kmucin values as high as 1018 and 1020 M−1 large numbers of viruses (20,330 and 180 respectively) are still free at mucin: virus ratios of 10:1. According to the model in Fig. 1, in the 0.314 dm3 volume of simulated intestine, large excesses of mucin over virus are needed to make a significant impact on the fraction of free virus at the lower Kmucin values. For example, at Kmucin values of 1013 M−1 and 1011.7 M−1 4% and 50% of the virus is free, respectively, at mucin: virus ratios of 1000:1. Thus mopping up virus is relatively inefficient and, except at very high Kmucin values, large numbers of virus may break through. This reflects the dilution in the large volume of the simulated intestine. The main conclusion from Fig. 1 is that as the mucin: virus ratio increases from < 1:1 through 1:1 to > 1:1 (i.e. going from left to right along the x-axis) for high Kmucin values ( > 1015 M−1), then the fraction of free virus falls dramatically and non-linearly, with all of the virus being mopped up at Kmucin of ∼1022 M−1.

The virus binding capacity of the mucin is finite

As expected at mucin: virus ratios < 1:1, the fraction of free virus approaches 100%, even with very high Kmucin values (Fig. 1). This simply reflects the fact that there is not enough mucin capacity to mop up all the virus such that very high virus challenge doses overwhelm the mucin barrier. With very high Kmucin values (1022 M−1), a very slight excess of mucin over virus is sufficient to mop up all the virus. At intermediate Kmucin values (1018–1020 M−1) mopping up the last remaining viruses is less efficient and proportionately higher mucin: virus ratios are required at lower Kmucin values.

Simplifying the approach: using Eq. (11) to model Fv in MRA

The total mucin: total virus ratio is also expressed as the concentration of total mucin, [Muctotal], on the x-axis scale of Fig. 1. This is relative to the virus concentration in the intestine which is fixed at 2.2 × 10−15 M. Thus 4.15 × 108 virus particles in a volume of 0.314 dm3 on dividing by L represent a concentration of 2.2 × 10−15 M. The dashed lines in Fig. 1 show Fv as calculated by Eq. (11) for each Kmucin value using [Muctotal] as an approximation for [Mucfree]. This is an appropriate approximation at high mucin to virus ratios. However, at mucin: virus ratios of < 1:1, [Mucfree] becomes much less than [Muctotal] as all the mucin is bound, particularly at high Kmucin values, and Eq. (11) is not applicable. As an example, the arrow in Fig. 1 shows that Eq. (11) predicts that 4.4% of virus is free (i.e. 95.6% of virus is bound) when there is 10,000-fold more virus particles than mucin molecules. This is clearly not possible. At lower Kmucin values, Fv decreases linearly with increasing mucin: virus ratio and may be modelled by Eq. (11) in good agreement with Fv calculated by the difference equation approach (symbols in Fig. 1). Thus Eq. (11) holds at low Kmucin values or at high mucin: virus ratios (irrespective of Kmucin). This is important because high mucin: virus ratios might be expected in a natural infection situation even with the high virus challenge dose used in the simulation here. From the practical point of developing MRA methodology, applying Eq. (11) is easier than the complicated difference equation approach developed for the symbols in Fig. 1. Thus understanding the limitations of Eq. (11) is important. In summary Eq. (11) fails at mucin: virus ratios of < 1:1 at the higher Kmucin values, as represented by the arrow, but is generally applicable at lower Kmucin values particularly at higher mucin: virus ratios. The results show that the value of Kmucin is critical not only in determining whether Eq. (11) can be used in MRA but also in assessing how effective the mucin barrier is in mopping up viruses.

Fc is related to both Ka and the virus dose in the intestine

The fraction of host cells with bound virus increases linearly with the dose of virus in the simulated intestine for all Ka values until at high virus doses the host cells become saturated (Fc → 1) with virus (Fig. 2). At low Ka values (<1013 M−1) unrealistically high virus doses (>1012) are required for saturation of the host cells. From Fig. 3, Fc increases linearly with Ka at a given virus dose but, not surprisingly, is limited by the virus dose at Vintestine: Ctotal ratios of < 1:1. Thus, there are 4.15 × 108 cells in the simulated intestine, and for the virus doses of 103 and 106 virions, the maximum values achievable for Fc are 2.4 × 10−6 and 2.4 × 10−3 respectively when all those virions are bound. At Ka values of > 1015 M−1, the binding is so strong that all the virions in the dose are bound (Fig. 3) and Fc reaches its plateau. At Vintestine: Ctotal ratios > 1 (i.e. the 4.15 × 1010 and 4.15 × 1011 doses in Fig. 3 representing Vintestine: Ctotal ratios of 100:1 and 1000:1, respectively) saturation of the host cells occurs at high Ka values with progressively lower Ka values required for saturation of the host cells as the virus dose increases. Thus increasing Ka above certain values has no effect on Fc either because all the viruses are bound at low virus doses (i.e. Vintestine: Ctotal ratios < 1) or because all the host cells are saturated at high virus doses (i.e. Vintestine: Ctotal ratios > 1). This is borne out in Fig. 2 which shows that for Ka values > 1015 M−1, the dose-response curves are superimposed such that Fc is limited by available dose and not by Ka. The main conclusion of Figs 2 and 3 is that the value of Ka over the range 105 to 1015 M−1 is critical in affecting Fc.

Simplifying the approach: using Eq. (7) to model Fc in MRA

The total virus dose in the intestine (Vintestine) may be approximated to Vfree, converted to a concentration and used as [Vfree] in Eq. (7) to calculate Fc as an alternative to the difference equation approach used to construct Fig. 2. However, this simplification fails at Ka values greater than 1015 M−1 (results not shown) for which Eq. (7) overestimates Fc by orders of magnitude at low virus doses. The results show that Eq. (7) is appropriate to calculate Fc up to Ka values of at least 1013 M−1 irrespective of dose.

Demonstration of the potential application of the model

The application of the model is demonstrated in Table 2 for four scenarios representing different combinations of initial challenge dose and affinity of the virus for the cell (as determined by Ka). Thus the Ka is 1010 M−1 for the low affinity binding virus and 1013 M−1 for the high affinity binding virus. These broadly reflect those measured experimentally (Nunes-Correia et al., 1999) for low binding and high affinity binding sites for influenza A virus H5N1 binding to canine kidney cells (Table 3 ). Also as discussed below, this difference of 1000-fold in Ka could reflect the result of a single amino acid change in the GP or Cr molecule affecting a salt bridge at the GP/Cr interface. For the purpose of the demonstration, pcell is assumed to be 0.1, i.e. a cell with bound virus has a 10% chance of successful infection. For these scenarios, the predicted probability of the host being infected, phost, ranges from 4.4 × 10−5 to 0.988. Fv is constant at 0.19 (Table 2) and is calculated using Eq. (11) with a Kmucin of 109 M−1 (at which Eq. (11) is appropriate) and a [Mucfree] of 4.25 × 10−9 M. This value of [Mucfree] is calculated using a mucin mass concentration of 1.7 mg/ml for saliva (Kejriwal et al., 2014), a dilution factor of 100-fold in the food or water matrix and a mucin protein molecular weight of 4000,000 Daltons (Kesimer and Sheehan 2012). It is solely calculated for demonstration purposes in Table 2 in the absence of published data. Fc in Table 2 varies for each scenario and is calculated using Eq. (7) with [Vfree] calculated from Vintestine (in the 0.314 dm3 volume of the simulated intestine) and L. The values of Ka used in these scenarios are low enough such that Eq. (7) is applicable (see above) as can be seen for the Ka = 1013 M−1 line from the difference equation approach in Fig. 2 where Fc ∼ 1.0 × 10−7 for a Vintestine of 1904 virions in agreement with Table 2.

Table 2

Demonstration of the application of the mechanistic dose-response model: Predicting the probability of infection of the host (phost) from the initial challenge dose (Vinitial) for a low and a high affinity binding virus.

Scenario	A: High dose, low affinity	B: High dose, high affinity	C: Low dose, low affinity	D: Low dose, high affinity.
Vinitial	10,000	10,000	100	100
Kmucin (M−1)	1.0 × 109	1.0 × 109	1.0 × 109	1.0 × 109
[Mucfree] (M)	4.25 × 10−9	4.25 × 10−9	4.25 × 10−9	4.25 × 10−9
Fv (Eq. (11))	0.1904	0.1904	0.1904	0.1904
Vintestine (Eq. (1))	1904	1904	19.04	19.04
Ka (M−1)	1.0 × 1010	1.0 × 1013	1.0 × 1010	1.0 × 1013
Fc (Eq. (7))	1.0 × 10−10	1.1 × 10−7	1.0 × 10−12	1.0 × 10−9
Ctotal	4.15 × 108	4.15 × 108	4.15 × 108	4.15 × 108
C.V (Eq. (3))	0.042	41.8	0.00042	0.42
phost (Eq. (2))	4.4 × 10−3	0.988	4.4 × 10−5	4.3 × 10−2
Calculation of phost using conventional dose response model (Eq. (16)) with p1 calculated from the ID50 obtained using the mechanistic method
Vinitial	10,000	10,000	100	100
p1	4.4 × 10−7	4.4 × 10−4	4.4 × 10−7	4.4 × 10−4
phost	4.4 × 10−3	0.988	4.4 × 10−5	4.3 × 10−2

Table 3

Summary of three approaches to parameterise Ka.

Term		Comments	Examples
Approach 1: Experimental measurement
Ka		Determined experimentally, e.g. for influenza A H5N1 virus binding to canine kidney cells	2.7 × 1012 M−1 and 2.0 × 1010 M−1 at 37 °C for high and low affinity binding sites, respectively (Nunes-Correia et al., 1999)
Approach 2: Calculation of Ka from published Kd data
Kd		Lot of experimental data for GP/Cr interactions	Range from 10−4 M to 10−12 M (see text)
n		Number of GP/Cr molecular contacts made on virus attachment to cell	Ranges from 1 to multiple contacts as suggested for EBOV. Carneiro et al. (2002) estimated 7 for VSV binding to PL bilayers.
ΔSa		Large and negative (see below) hence lowering magnitude of Ka
Ka		Calculated from Kd, n and ΔSa using Eq. (17)	*For n = 3 and Kd = 10−4 M, Ka = 1012 M−1. *For n = 1 and Kd = 10−12 M, Ka = 1012 M−1.*assuming ΔSa = 0 J/mol/K
Approach 3: Calculation of Ka from enthalpy and entropy terms
Enthalpy term, ΔHa		Favourable interactions (e.g. formation of salt bridges) between amino acid residues at GP/Cr contact interface and good spatial fits give a large negative value driving virus binding. Limited by lack of specific data for GP/Cr contacts or for virus/host cell binding	−56.5 kJ/mol for antibody binding to its antigen (Bostrom et al., 2011). −4958 kJ/mol for VSV binding to PL bilayers (Carneiro et al., 2002)
Entropy term, ΔSa	Component entropy terms
	ΔSsolvent	Likely to be positive as disordering and hence entropy of water solvent molecules displaced from GP/Cr contact surfaces during binding may increase	+401 J/mol/K for antibody binding to its antigen (Bostrom et al., 2011)
	ΔSconf	Negative due to reduction in conformational mobility. Estimated to be −6.1 J/mol/K per amino acid for ordering of disordered regions of proteins (Rajasekaran et al., 2016)	−301 J/mol/K for antibody binding to its antigen (Bostrom et al., 2011)
	ΔSrt	Negative for all interactions when a particle is immobilised on a surface.	−209 J/mol/K for a pentapeptide binding to a lipid membrane (Ben-Tal et al., 2000). −33.5 J/mol/K for antibody binding to its antigen (Bostrom et al., 2011)
	ΔSmem	Negative giving a repulsive force pushing virus away from host cell (Sharma, 2013)	No data
Overall ΔSa calculated as sum of component entropy terms		Large and negative reflecting general immobilisation of components during binding.	−16,062 J/mol/K for VSV binding to PL bilayers (Carneiro et al., 2002)
ΔGa		Calculated from ΔHa, ΔSa and T using Eq. (18)
Ka		Calculated from ΔGa using Eq. (19)

Effect of stochasticity

The approach using Fc from Eq. (7) predicts fractions of a C.V (Table 2) for both low dose scenarios and also for the high dose, low affinity scenario. In reality CV must be an integer as in the difference equation approach. The issue of stochasticity arises in the model application in Table 2 when C.V is a fraction. To investigate the effect of stochasticity, the mechanistic model developed in Table 2 is used to determine the infectious dose 50% (ID50), i.e. the value of Vinitial for which phost = 0.5. It was found that for the high affinity virus (Ka = 1013 M−1), the ID50 is 1574 viruses and for the low affinity virus (Ka = 1010 M−1) the ID50 is 1574,000 viruses (Table 4 ). The C.V value for the ID50 is 6.6 and, being much greater than 1, is in the range where stochasticity is less of an issue. The classic dose response model for infection of a host is written as:-where p1 is the risk of infection from ingestion of a single virion by the host. By rearranging Eq. (16) and setting Vinitial = ID50 and phost = 0.5, p1 is calculated as 4.4 × 10−4 and 4.4 × 10−7 for the high affinity and low affinity viruses respectively. These values may be used directly in the conventional dose-response model in the form of Eq. (16). Furthermore, setting Vinitial = 1 in the mechanistic model in Table 4 predicts the same probability values of 4.4 × 10−4 and 4.4 × 10−7 for the high affinity and low affinity viruses respectively suggesting that stochasticity is not a problem as expected for dose-response models which are linear at low doses.

Table 4

Obtaining parameters from the mechanistic dose-response model for use in conventional dose-response models in the form of Eq. (16): The infectious dose 50% (ID50) and the probability of infection from ingestion of a single virion by the host (p1) for a high and low affinity virus.

Scenario	ID50, low affinity	ID50, high affinity	Dose of one low affinity virus	Dose of one high affinity virus
Vinitial	1574,000	1574	1	1
Kmucin (M−1)	1.0 × 109	1.0 × 109	1.0 × 109	1.0 × 109
[Mucfree] (M)	4.25 × 10−9	4.25 × 10−9	4.25 × 10−9	4.25 × 10−9
Fv (Eq. (11))	0.1904	0.1904	0.1904	0.1904
Vintestine (Eq. (1))	299,809	299.8	0.19	0.19
Ka (M−1)	1.0 × 1010	1.0 × 1013	1.0 × 1010	1.0 × 1013
Fc (Eq. (7))	1.6 × 10−8	1.6 × 10−8	1.0 × 10−14	1.0 × 10−11
Ctotal	4.15 × 108	4.15 × 108	4.15 × 108	4.15 × 108
C.V (Eq. (3))	6.6	6.6	4.2 × 10−6	4.2 × 10−3
phost (Eq. (2))	0.5	0.5	4.4 × 10−7	4.4 × 10−4

Applying the outputs of the mechanistic model to the MRA

The values of phost are also calculated in Table 2 by using p1 (obtained from the mechanistic model as described above) directly in Eq. (16) to estimate the risk from the initial challenge dose (Vinitial). The results are identical with those predicted by the mechanistic method. This not only provides further evidence that stochasticity is not an issue but also demonstrates how the output from the mechanistic model can be applied directly to conventional dose-response models.

Theory

Parameterising the model: determining the virus/host cell association constant, Ka

Three approaches for parameterising Ka are set out in Table 3. In a few cases, the Ka is determined experimentally by directly measuring virus binding to cells (Table 3). Approach 2 developed here is to estimate Ka from published Kd values which are routinely measured experimentally. Thus, the magnitude of Ka for binding of a virus to its host cell is determined by the strength of the interaction(s) between viral surface components, typically the viral glycoprotein (GP) in the case of enveloped viruses such as EBOV, and the receptors (Cr), e.g. T-cell immunoglobulin and mucin domain protein 1 (TIM-1) for EBOV, on the human host cell surface (Yuan et al., 2015). Biochemists quantify the binding affinity between molecules in terms of the dissociation constant (Price and Dwek, 1979):for the reversible dissociation process given by:-

The magnitude of Kd in units of mol/dm3 (M) may be determined experimentally between purified parts of the virus GP and the host Cr using surface plasmon resonance (SPR) in which one protein “partner” is immobilized to a chip surface and changes in refractive index are used to measure binding of the other “partner” which is free in solution. Thus, Yuan et al. (2015) immobilized EBOV-GP to a chip and used SPR to measure a Kd of 2.67 × 10−5 (M) for binding of soluble human TIM-1 from solution. EBOV entry involves a second GP/Cr binding step within the endosome in which the GP binds to the host Niemann-Pick C1 (NPC1) protein. Using SPR, Wang et al. (2016) reported a Kd of 1.58 × 10−4 (M) for EBOV GP binding to human NPC1. These binding affinities for EBOV GP are much weaker than that for MERS-CoV binding to its cellular receptor CD26 for which the Kd was 1.67 × 10−8 (M) as measured by SPR (Lu et al., 2013). Raman et al. (2014) summarise Kds for influenza virus HA binding to host cell glycans on human tracheal and alveolar sections. Some HAs bind to human receptors with Kds in the 10−12 M range, while others show Kds of ∼10−9 M. Gambaryan et al. (2005) present Kds for binding of H1 and H3 influenza viruses to glycans, typically in the range of 10 – 100 × 10−9 M.

Taking into account multiple interactions between the virus and host cell

For a single GP/Cr interaction, Ka is the reciprocal of Kd. However, the reciprocal Kd values measured by SPR for the GP/Cr interaction cannot directly be used as the Ka values in Figs. 2 and 3. This is because Kd as measured by SPR does not take into account the number of GP/Cr contacts per virus/cell interaction and also does not allow for any changes in the entropy (S) on the binding and immobilization of a whole virus particle to a cell (discussed below). For example, Beniac et al. (2012) calculated that an EBOV virion filament of 982 nm in length would have about 1888 copies of the GP spike protein and could therefore make multiple contacts with Cr molecules on the host cell surface. The change in Gibbs free energy (ΔGd) for dissociation of one mole of GP/Cr contacts is related to Kd (Price and Dwek, 1979, Kastrits and Bonvin, 2013) by:-where T is the temperature and R is the ideal gas constant (8.31 J/mol/K (Price and Dwek 1979)). Since the Gibbs free energy (G) is a state function in thermodynamics, the overall change in G for dissociation of n moles of GP/Cr contacts, is given by

Since the change in Gibbs free energy (ΔGa) for association of one mole of virus with one mole is cells is related to Ka by (Gale, 2017) and the change in Gibbs free energy for association of n mole of GP/Cr contacts, (Kastrits and Bonvin, 2013) then Taking into account the change in entropy on whole virus binding, ΔSa, (discussed below) Ka may be expressed as:-

The value of G is affected by temperature and pressure, and for this reason G is used by biochemists. Thus providing the temperature and pressure are constant, as assumed here for virus in the host intestine, then any changes in G relate only to those changes in energy within the system of interest, namely from the molecular interactions between virus and host during infection within the simulated intestine.

Estimating Ka through the thermodynamic parameters, enthalpy and entropy

Approach 3 in Table 3 is to calculate the association constant Ka through the changes in two thermodynamic parameters, namely enthalpy (H) and entropy (S) on binding. Thus Ka, is related to ΔGa for association of the virus with the cell (Gale, 2017, Kastrits and Bonvin, 2013) according to:which on rearranging giveswhere ΔH and ΔS are the changes in enthalpy and entropy, respectively, of the virus/host cell system on association (Carneiro et al., 2002, Gale, 2017). Therefore in principle knowing ΔH and ΔS would enable calculation of Ka and hence Fc for a given dose of virus from Fig. 2. Bostrom et al. (2011) demonstrated how the binding affinity of an antibody could be broken down into the enthalpy and entropy terms. How this approach could be applied to a virus during infection with the aim of calibrating dose-response models is summarized as Approach 3 in Table 3. To the author's knowledge there are no data for ΔH and ΔS for viruses binding to host cells. However, Carneiro et al. (2002) measured the forces between vesicular stomatitis virus (VSV) and an artificial phospholipid (PL) bilayer (as opposed to a host cell) and obtained values for ΔHa and ΔSa of −4958 kJ/mol and = −16,062 J/mol/K respectively. Carneiro et al. (2002) explain the huge value for ΔHa as the result of the cumulative ΔHas from multiple interactions due to the dense packing of G protein on the VSV virus surface. Indeed Carneiro et al. (2002) estimate seven G proteins’ binding per VSV and calculate the ΔHa and ΔSa values for binding of a single G protein at −694 kJ/mol and −2196 J/mol/K respectively.

The importance of understanding entropy changes during virus binding to a cell for calculating Ka

Approaches 2 and 3 in Table 3 are dependent on estimating ΔSa. As shown by Carneiro et al. (2002) for VSV binding to PL bilayers, the decrease in entropy is huge (Table 3), demonstrating the importance of understanding ΔSa in estimating Ka. Unlike ΔHa which is based on molecular contacts at the binding faces, ΔSa involves changes in order and mobility and is difficult to visualise, but may be broken down into a number of components which act additively (Bostrom et al., 2011). Thuswhere ΔSsolvent is the change in entropy of the water solvent molecules, ΔSconf is the change in the internal conformational freedom of the proteins, ΔSrt is the change in rotational and translational freedom of the virus, and ΔSmem is an entropic pressure associated with bringing two membranes close together, as for example, when the virus envelope approaches the host cell membrane prior to fusion (Sharma, 2013). The likely contributions from the four entropy component terms are summarized in Table 3.

The change in conformational entropy (ΔSconf): intrinsically disordered proteins in viruses and their role in ligand binding

Many ligand binding proteins have intrinsically disordered regions (IDRs) which undergo a disorder-to-order transition on or near the interface on binding the ligand (Fong et al., 2009) decreasing the entropy such that ΔSconf is negative (Rajasekaran et al., 2016). Viral proteins are rich in intrinsic disorder (Dolan et al., 2015) and intrinsically disordered proteins serve as host cell receptors for viruses. An example is the ephrin receptor ligand binding domain (Fong et al., 2009) which also serves as receptor for Hendra virus (HeV) and NiV (Xu et al., 2012). Thus the unbound ephrin receptor contains partially disordered loops. In the complex (bound to ephrin) these loops are ordered to form the ligand-binding channel (Fong et al., 2009). This raises the question of whether the ephrin receptor undergoes a disorder-to-order transition on HeV/NiV binding. Interestingly Xu et al. (2012) show from the X-ray crystal structure that the binding of HeV G protein involves the movement of a tryptophan “latch” on the ephrin-B2 receptor. There would be a change in the conformational entropy ΔSconf associated with this during formation of the HeV G protein/ephrin-B2 receptor complex.

According to Eq. (17), increasing the number of GP/Cr contacts during virus binding increases Ka and hence the infectivity according to Fig. 3. Counterintuitively, in the case of EBOV GP pseudoviruses, increasing the surface density of GP actually decreases the infectivity by about 10-fold, perhaps reflecting steric hindrance from tightly packed GP proteins (Mohan et al., 2015) blocking a conformational change that gives an increase in ΔSconf required for binding of GP to TIM-1 (Gale, 2017). This is consistent with the observation that EBOV GPs are well separated in space on the EBOV filament surface (Beniac et al., 2012) thus allowing plenty of space for conformational changes without steric hindrance. A demonstration of how changes in ΔSconf could explain the 10-fold increase in infectivity reported by Mohan et al. (2015) on decreasing the GP density may be explored with Eq. (19). Thus increasing ΔSa by 20 J/mol/K (irrespective of the values of ΔHa and T) increases Ka by 11-fold. A decrease in conformational entropy (ΔSconf) due to steric hindrance of 20 J/mol/K represents the ordering of just four amino acid residues on the basis of the −6.1 J/mol/K per residue proposed by Rajasekaran et al. (2016) for disordered proteins. Therefore blocking the disordering of just four amino acids in the EBOV GP could explain the 10-fold decrease in infectivity of EBOV GP pseudoviruses with high GP densities.

Booth et al. (2013) show that filovirus filaments are very flexible, more flexible than filamentous paramyxoviruses and much more flexible than rhabdoviruses. Thus a flexible filament, making multiple GP/Cr contacts, can adopt many more conformations in 2D-space when bound to a host cell surface than a rigid rod of the same length. A bound filament will therefore have a higher entropy than a stiff rod, and will have a higher Ka. A flexible filament can also test out multiple GP/Cr contacts, thus maximising contacts and making ΔHa more negative in magnitude which increases Ka according to Eq. (19). Indeed Booth et al. (2013) suggest that the extreme aspect ratio of filaments may be an adaptation that enhances cellular attachment.

The change in rotational and translational entropy (ΔSrt)

The association of two species to form a bound complex, e.g. the binding of a ligand to a protein or the adsorption of a peptide on a lipid membrane, always involves an entropy loss (Ben-Tal et al., 2000). This no doubt applies to the binding of viruses to cells. Consideration of ΔSrt may be of particular importance for binding of filoviruses (Gale, 2017) which comprise repeating modular units linked together into single very long filaments (Beniac et al., 2012, Booth et al., 2013). Thus Gale (2017) argued that the prior immobilisation of multiple virus units through linking together into a polyploid filament enhanced cell binding compared to that of single virus units in the case of filoviruses through a less unfavourable ΔSrt term. Furthermore for this reason, natural polyploid filovirus filaments may bind more strongly to cells than suggested by results from experiments to determine filovirus infectivity using EBOV GP-expressing pseudoviruses, which are spherical and single.

Using enthalpy to relate changes in amino acid sequence at the binding faces of the virus GP and host receptor protein to a change in Ka

While lack of data on ΔHa limits the use of Approach 3 in Table 3 to determine Ka, Eq. (19) may be applied to determine the impact of amino acid changes in GP and Cr on Ka and hence link genetic sequencing data to the risk of infection. The ultimate aim would be to predict how such sequence changes affect infectivity, host range and jumping the species barrier from the perspective of the dose-response. This could include assessing the relative risk of novel strains of virus with mutations in the GP and also assessing the risk of a “standard” virus jumping the species barrier into a novel host which may have amino acid differences in its Cr relative to the “standard” host. The starting point is having information on the infectivity of the “standard” virus in the “standard” host and also crystal structure data of the “standard” virus GP bound to the “standard” host Cr. Crystal structure data are available for a number of viral GPs docked to their human host cell Cr at atomic detail including EBOV and MERS-CoV (Zhao et al., 2016, Wang et al., 2016, Lu et al., 2013). Knowing the crystal structure of the GP/Cr protein complex assists in understanding how changes in amino acids at the binding interface of GP and Cr affect ΔHa and hence Ka through Eq. (19). Thus although ΔHa itself may not be known, it may be possible to predict the change in ΔHa (i.e. ΔΔHa) for an amino acid substitution from basic biophysics knowledge of the energies of the intermolecular forces that hold proteins together. These are typically salt bridges and hydrogen bonds as shown for MERS-CoV bound to its Cr (Lu et al., 2013). An example of how this could be applied is demonstrated in Table 5 for changing amino acids involved in salt bridges. In a salt bridge a negatively charged amino acid residue on one protein interacts electrostatically with a positively charged residue on the other protein resulting in a strong attraction. Salt bridges, e.g. between histidine 30 (positively charged) and aspartate (Asp) 70 (negatively charged) in bacteriophage T4 lysozyme contribute −11 to −21 kJ/mol to ΔG (Anderson et al., 1990, Dong and Zhou, 2002). For the purpose of the demonstration in Table 5 it is assumed that the ΔHa for a salt bridge formation is −17.8 kJ/mol as this gives a 1000-fold change in Ka (Price and Dwek 1979) for each salt bridge removed or added through mutation. Proximity of two residues of the same charge results in electrostatic repulsion as in Scenarios A and B which represent viruses with very weak affinities for the host. Scenarios D and E represent viruses which have adapted to the host with hugely increased affinities. Scenarios A and E are only two mutations away from each other. Thus changing two positively charged amino acids (Scenario A) into negatively charged amino acids (Scenario E) could increase Ka by 1012-fold.

Table 5

Theoretical consideration of how changes in ΔHa (ΔΔHa) through mutations affecting two amino acid residues involved in salt bridges at the contact interface between GP and Cr could affect the binding affinity Ka at 310 K (37 °C).

Virus/host scenario	Effect of mutation on the number of charged amino acid residues opposite each other at the contact surface of GP and Cr	Representation (+, positive amino acid residue; −, negative amino acid residue; 0, electrostatically neutral amino acid residue)	Comments	ΔΔHa (kJ/mol) relative to Scenario C. Calculated assuming one salt bridge contributes 17.8 kJ/mol (Anderson et al., 1990)	Effect on Ka from Eq. (19) relative to Scenario C	Increase in Ka from Eq. (19) relative to Scenario A
A	Presence of 2 similarly charged residues on each surface in close proximity	GP….Cr   + ….. +  + …. .+	Very strong electrostatic repulsion, e.g. as proposed for an insect virus jumping the species barrier to a mammal (Gale 2017)	2 x +17.8 = +35.6	106 fold decrease	1
B	Presence of 1 similarly charged residue on each surface in close proximity	GP….Cr 0 ….. +  + …. .+	Strong electrostatic repulsion	1 x + 17.8 = + 17.8	103-fold decrease	103
C	0 ionic interactions	GP….Cr 0 ….. + 0…. .+	No effect	0	1	106
D	Generation of 1 salt bridge	GP….Cr  -….. + 0…. .+	Strong electrostatic attraction in a virus partially adapted to its new host	1 x −17.8 = −17.8	103-fold increase	109
E	Generation of a second salt bridge	GP….Cr -….. + -…. .+	Very strong electrostatic attraction in virus fully adapted to new host	2 x −17.8 = −35.6	106-fold increase	1012

EBOV adaptation to humans through changes in GP

Amino acid changes in the EBOV GP affect its host specificity. Thus, EBOV Makona in the 2013 West Africa outbreak adapted to humans (and in doing so became less infectious to bats) through substituting alanine (A) with valine (V) at residue 82 of the GP (Diehl et al., 2016, Urbanowicz et al., 2016). This is referred to as the A82V substitution. However, the changes in infectivity to humans (as measured by EBOV GP pseudoviruses infecting human cell lines) appear to be relatively small and in the range of 2-fold to 4-fold. A 4-fold increase in infectivity due to a 4-fold increase in Fc in Fig. 3 would reflect a 4-fold increase in Ka, and a change in ΔGa according to Eq. (18) of −3.6 kJ/mol. This would probably reflect changes in ΔHa (and not the T.ΔSa term) since residue 82 on the EBOV GP is between residues 80 and 83 which make multiple atom contacts directly with six surface residues on the NPC1 protein (Wang et al., 2016). Thus, the A82V substitution in the GP may change ΔHa by ∼−3.6 kJ/mol on the basis of measured changes in infectivity. This relatively small change (compared to the −11 to −21 kJ/mol for a salt bridge) is consistent with the chemical similarity of alanine and valine.

Changes in the host receptor (NPC1) that protect Eidolon helvum fruit bats from EBOV infection

In the human NPC1 protein (Cr), an Asp at residue 502 protrudes from the surface and contacts the EBOV GP in the docked complex (Zhao et al., 2016). The natural presence of phenylalanine (Phe) at residue 502 (instead of Asp) in the NPC1 of the African Straw-coloured fruit bat (Eidolon helvum) appears to protect this species from EBOV (Ng et al., 2015). Indeed, replacing Phe at residue 502 of E. helvum NPC1 with Asp completely restored its binding to EBOV GP (Ng et al., 2015). From the X-ray structure, the presence of Phe at residue 502 of the NPC1 in E. helvum would cause severe clashes on binding EBOV GP (Zhao et al., 2016) such that ΔHa in the E. helvum NPC1/EBOV GP complex would be significantly less negative than that for binding to human NPC1 and perhaps even positive (reflecting a repulsion). Indeed, determining ΔΔHa would enable the effect of the Asp to Phe substitution at residue 502 on Ka and hence on the infectivity of EBOV to E. helvum to be assessed through Eq. (19) as in Table 5.

Hoffmann et al. (2016) developed VSV pseudoviruses with filovirus GPs expressed on their surfaces to test the efficiency of GP-mediated cell entry into cell lines of different bat species. Hoffmann et al. (2016) demonstrated that the efficiency of entry of EBOV-GP VSV pseudovirus into E. helvum cell lines was markedly reduced compared to that with GPs from other filoviruses and at standard virus levels no infection was detected. However, at high levels of (i.e. undiluted) virus, EBOV-GP VSV pseudovirus did show significantly increased infection compared to the negative controls. Thus EBOV GP is capable of mediating entry into E. helvum cells, albeit with low efficiency. This is entirely consistent with the dose-response model simulated in Fig. 2 which predicts for low Kas (105–107 M−1) that very small fractions (10−7 to 10−5) of the host cells have bound virus at very high virus doses, i.e. 1011 to 1012. Thus some infection would be expected in E. helvum on the basis of Eqs. (2) and (3) and the probability of infection which, although low at low Ka values, would never be zero because there is no threshold Ka below which binding does not occur (Fig. 3). Supporting this, Hayman et al. (2010) reported EBOV antibodies in one of 256 E. helvum bats tested in Ghana and suggested this resulted from EBOV infection.

Parameterising the model: determining values for Kmucin

Values for Kmucin may be derived from published data on the strength of interactions between lectins and glycans. Lectins are proteins which bind selective glycan groups (Taylor and Drickamer, 2006). Viruses such as influenza virus, NoV and rotavirus behave as lectins by binding selectively to glycans. Indeed, Kmucin is related to the number of lectin/glycan interactions and their respective Kds in the same way as Ka according to Eq. (17). Dam and Brewer (2010) report Kds of 2.0 × 10−10 M for the lectin soy bean agglutinin binding to porcine submaxillary mucin which has multiple GalNAc sugars to bind. Similarly Vataira macrocarpa lectin has a Kd of 1.0 × 10−10 M with porcine submaxillary mucin. Acting singly (i.e. using n = 1 in Eq. (17)), these Kds would translate into Kmucin values of ∼1010 M−1 at which very high excesses of mucin would be needed to achieve removal of the virus according to Fig. 1. There are two mechanisms by which the magnitude of Kmucin could be increased:-

Through multiple binding interactions; and

Through irreversible binding.

The effect of multiple contacts on the strength of binding has been shown for binding of the mannose binding protein (MBP) to carbohydrate structures on the surfaces of pathogens. Each carbohydrate-recognition domain (CRD) of the MBP interacts only with the terminal sugar residue in an oligosaccharide chain (Taylor and Drickamer, 2006). The Kd for interaction with a high mannose oligosaccharide is ∼10−3 M, i.e. low affinity. High-affinity binding of MBP requires interaction of multiple CRDs with multiple terminal mannose residues. The trimeric structure of MBP presents a cluster of CRDs for interaction with appropriately spaced terminal mannose residues on the pathogen surface. Thus the arrays of terminal sugar residues on the surfaces of microorganisms can interact with multiple sites simultaneously. As shown in Eq. (17), a three way interaction involving three binding sites with Kds of 1.0 × 10−3 M can result in an overall Ka of up to 1.0 × 109 M−1 for a multivalent ligand (i.e. pathogen surface) with appropriately spaced terminal sugar residues (Taylor and Drickamer, 2006). In the case of influenza virus A, the binding site of each HA polypeptide is relatively shallow and interacts primarily with the terminal sialic acid residues linked to galactose (Taylor and Drickamer, 2006). The affinity for monomeric sialosides is weak (Kd ∼10−3 M), but binding to cell surfaces is enhanced by the simultaneous interaction of multiple HA-binding sites with multiple sialic acid residues on the target host cell (Taylor and Drickamer, 2006). The cumulative effect of n binding sites each with the same Kd is given by the n th power (Eq. (17)). Thus seven interactions with Kd ∼10−3 M, for example, would give a Ka ∼1021 M−1. It is concluded that mucins and MBP PRRs of the innate immune system could bind virus with sufficient affinity such that slight mucin excesses could remove all the virus (Fig. 1). Furthermore, Kmucin could become higher if the virus were irreversibly bound for example by being “folded in” to the inside of mucous droplets such that it can no longer exchange with free virus on the outside. In effect koff would tend to 0 in Eq. (6) and hence Kmucin would tend to infinity. It should be noted that the magnitude of Kd (and hence Ka and Kmucin) is constant for a given molecular interaction at a given temperature. However, the magnitude of Kmucin may change during the infection process through both host mechanisms and virus mechanisms as is now discussed.

Changes in Kmucin through host effects

Mucin glycosylation changes during infection such that mucins isolated from rotavirus-infected mice at 4 days post infection were more potent at inhibiting rotavirus infection than mucins from control mice (Boshuizen et al., 2005). Interestingly there are also age-dependent differences in mucin quantities, composition, and/or structure which alter the antiviral capabilities of the small intestine mucins (Boshuizen et al., 2005). Furthermore, mucin production by the host increases during infection by rotavirus with mucin-coding mRNA levels peaking at 1 day post infection (Boshuizen et al., 2005). Thus both the mucin: virus ratio and the magnitude of Kmucin could change during progression of the infection, hence affecting Fv in Fig. 1.

Changes in Kmucin through virus effects

There are viral mechanisms that remove the virus from the mucin so increasing the concentration of free virus. For example in the case of influenza virus, the neuraminidase removes sialic acids from glycans, which enables virus particles to be removed from the cell surface after assembly and from decoy receptors e.g. in mucus (Guo et al., 2017). This has two effects on virus binding in Fig. 1. First it increases the rate of removal of the virus (koff in Eq. (6)) and hence decreases Kmucin. Second it decreases the concentration of mucin thus decreasing the mucin: virus ratio. From Fig. 1 both these act synergistically to increase the fraction of free virus.

In the case of influenza viruses, the HA protein binds to sialoside receptors on the host cell surface, preferentially binding to sialic acids linked to a penultimate galactose (Guo et al., 2017). The balance between activities of HA and neuraminidase proteins has a critical role in optimal viral fitness, tropism and transmission (Guo et al., 2017). The opposing effects of HA and neuraminidase in influenza A virus infectivity could be modelled through their effect on Kmucin (through kon and koff in Eq. (6)) and the mucin: virus ratio in Fig. 1. Handel et al. (2014) use differential equations to model how sticky an influenza virus should be to maximize fitness, with stickiness representing the balance between attachment and detachment. Shanker et al. (2011) proposed that local flexibility in part of the NoV surface protein could allow the virus to disassociate from salivary mucin-linked HBGA in the changing microenvironment (pH for example) during its passage through the GI tract to subsequently associate with HBGAs linked to intestinal epithelial cells. In effect the koff rates increase in Eq. (6) such that Kmucin varies in different parts of the host.

Discussion

This is a concept paper which presents a generic approach to use data from biochemistry and molecular biology to parameterize dose-response models for MRA through thermodynamic equations. In the model developed here a host intestine is simulated with 4.15 × 108 cells in a volume of 0.314 dm3. The initial infection process is broken down into two stages, namely escape of the virus from the innate host defense barriers (Fig. 1) and the subsequent binding of any remaining virus to its specific receptor on the host cell (Fig. 2). The two parameters, namely Fv and Fc, which quantify these stages represent probabilities which together with pcell can be used with Eqs. (1)–(3) to give a complete dose-response model which calculates the risk of infection of the host from a given challenge dose (Vinitial). (Although Vintestine from Eq. (1) does not appear in Eqs. (2) or (3), it is used in Fig. 2 to estimate Fc for Eq. (3)). The difference equation approach developed here is tedious to apply to MRA and applying Eq. (7) and Eq. (11) would greatly simplify the calculation of Fc and Fv, respectively for MRA purposes. However, using Eq. (7) for Fc fails at low doses of a virus which has a high binding affinity for the host cells (Ka > 1015 M−1) and the difference equation approach has to be used as in Fig. 2. Using Eq. (11) for Fv is appropriate for MRA if the mucin concentration exceeds the virus concentration in the host by a factor of > 10 or if Kmucin < 1013 M−1 (Fig. 1).

The mechanistic approach developed here may be used to determine the ID50 (see Table 4) from which it is easy to calculate p1, the risk to the host from ingestion of a dose of one pathogen. The parameter p1 may then be used directly in conventional dose-response models (Eq. (16)) thus linking the application of the mechanistic model developed here to MRA. The work here appears to demonstrate that stochasticity in the mechanistic model is not an issue. If stochasticity issues were to arise for low doses of very high affinity pathogens (i.e. those with very high Ka), the approach should be to use the difference equation approach to determine the ID50 from which p1 can then be calculated and used directly in Eq. 16. Thus, depending of course on the value of pcell, C.V is much greater than 1 for the ID50 (Table 4) and stochasticity issues are minimised.

Although more information on the concentrations of mucin proteins in the saliva and intestine is needed, lack of data on the thermodynamic parameters is the main limiting factor. Approach 3 (Table 3) which uses Eq. (19) to calculate Ka is limited by the current lack of data for ΔHa and ΔSa for the virus/host cell interaction and Approach 2 using published Kd data should therefore be pursued in the absence of experimental Ka data from Approach 1. While data are available for Kd for GP/Cr binding for several viruses, Approach 2 using Eq. (17) to estimate Ka is also limited by lack of information on ΔSa for binding of the whole virus to cells. This is important because the magnitude of ΔSa for immobilization of a whole virus is likely to be huge as demonstrated for VSV binding to artificial PL bilayers (Table 3). Lack of information on ΔSa is therefore a major data gap, and furthermore may also have wider implications in affecting our ability to interpret measured Kd values in terms of the actual infectivity of the virus. Thus it is the Ka as determined by Eq. (17) taking into account ΔSa for whole virus binding which defines infectivity (see results from Table 2) and not the Kd alone. This is because ΔSa for binding of the whole virus to cells will be much greater than the ΔSa for binding of a soluble protein fragment as in Kd determinations by SPR. Indeed the ΔSa for binding of VSV is a huge −16,062 J/mol/K compared to the ΔSrt and ΔSconf terms for binding of an antibody to its antigen (Table 3) which is comparable to GP binding to its Cr in terms of molecular sizes.

The model developed here based on Eq. (4) does not allow for multiple virus particles to bind to a single host cell. However, the cell membrane of the host cell typically contains many copies of Cr and in theory each host cell could therefore bind multiple copies of the virus. Furthermore, the virus particles themselves could be ingested as a clump. Indeed in the case of EBOV, many genome copies may be linked together into long polyploid filaments (Beniac et al., 2012) which could affect the thermodynamics of binding (Gale 2017). Eq. (4) could be modified to accommodate multiple viruses binding to a single host cell (as described for ligands binding to a protein by Price and Dwek (1979)) if the number of Crs per host cell is known. It should be stressed that homogeneity in both the mucin and pathogen is assumed here for simplicity. In further developing this mechanistic model, data would be needed not only on the heterogeneity of the mucin concentration but also on the statistical distribution for the number of viruses bound per cell. For example the distribution could be Poisson such that for a ratio of 1 virus bound per host cell (representing a high Vintestine), some host cells have zero, most have one while a few have two, three or even four viruses bound. Alternatively the distribution could be over-dispersed, such that a few host cells have thousands of bound virus while most have none. Clearly, over-dispersion would decrease C.V compared to the Poisson distribution. In the model presented here for demonstrating proof of principle, the modelling of Fc as function of Vintestine assumes each host cell binds just one virus through Eq. (4). This is appropriate for two reasons. First, exposure to faecal/oral pathogens through routes such as water may involve very low numbers of pathogen per person (see Gale 2017), such that it would be unlikely that more than one pathogen would bind to the same host cell. Second, it maximises C.V and hence phost. Thus, the binding of more than one virus to a host cell may be waste of virus resource because C.V in Eq. 2 is not maximised. Binding of two viruses to the same cell will not affect pcell, unless there is co-operation between the two viruses in some way, such that one facilitates infection by the other. In this sense, the value of phost predicted through the mechanistic dose-response here is worst case. It should be noted that Gale (2017) proposed that formation of polyploid filaments in the case of EBOV did indeed enhance cell binding by optimising ΔSa, in effect representing a co-operation between viral particles in the infection process. This demonstrates the importance of not only considering the statistical distribution aspects of virus clumps but also the thermodynamic implications on their binding.

As shown in Fig. 1, the probability of the virus escaping the mucin barrier (as represented by Fv) is controlled by the mucin: virus ratio and the mucin/virus binding affinity as represented by Kmucin. The range for which Kmucin is biologically significant is < 1022 (M−1); 1022 being the Kmucin value at which there is no free virus at mucin: virus ratios > 1:1 according to the simulation in Fig. 1. Values for Kmucin could be high due to multiple contacts with repeating units and “folding in” of bound virus into the interior of the mucus. Reading the x-axis from right to left in Fig. 1 suggests a threshold effect for virus dose, as the binding capacity of the mucin is exceeded at high virus doses. This is not considered further here, but is of interest for the development of MRA which generally assumes there is no threshold dose. The approach developed here in Fig. 1 for mucins could also be applied to modelling the probability of the virus being inactivated by other components of the host innate immune system such as the PRRs, which include the MBP, which recognize repeating units on the pathogen surface. Once bound to the MBP, the virus would be taken up (phagocytosis) through the complement system by macrophages and destroyed (Taylor and Drickamer, 2006).

The probability of a host cell having bound virus (as represented by Fc) increases linearly with dose (Fig. 2). The biologically significant range for Ka is 〈1015 (M−1) within which Fc increases linearly with increasing Ka (Fig. 3). Kd values measured by SPR are in the 10−9 M–10−12 M range for GP/Cr interactions for MERS-CoV and influenza virus suggesting very strong binding (Lu et al., 2013, Raman et al., 2014). Even for EBOV GP with Kds in the 10−4 to 10−5 M range (Yuan et al., 2015, Wang et al., 2016), high Ka values could be achieved through multiple contacts (Eq. (17)) although this may be offset by the ΔSa term (as discussed above). Increasing the Ka above 1015 M−1 does not further increase the fraction (Fc) of host cells with bound virus (Fig. 3) such that Fc versus virus dose curves are superimposed in Fig. 2 at Ka 〉 1015 M−1. At such high Kas, receptor binding is not likely to be rate-limiting in the infection process. Thus Xu et al. (2012) found that the stability of the HeV-G/ephrin-B2 association does not strongly correlate with the efficiency of viral entry, suggesting that, for ephrin-B2-expressing cells, viral attachment is not the rate limiting step in the viral entry process. This is consistent with highly efficient binding (i.e. large Ka) such that further increasing Ka does not increase the fraction of cells with bound virus, Fc, according to Fig. 3. Hence C.V (through Eq. (3)) is constant and not rate-limiting, while a subsequent part of the infection process accommodated in pcell is less efficient than Fc, such that pcell is very small and controls the overall value of phost according to Eq. (2). Indeed, the magnitude of Ka is much more important when receptor binding is inefficient such as in the jumping of the species barrier into a novel host by an emerging virus as suggested by Gale (2017). This is confirmed by Fig. 3 with Fc increasing linearly with Ka at lower Ka values. Gale (2017) demonstrated that virus binding to its receptor is only important in controlling the probability of infection of a host cell if it is a major barrier, i.e. relatively inefficient. The Ka for a virus jumping the species barrier into a novel host species would be expected to be very low perhaps due to electrostatic repulsion as for Scenarios A and B in Table 5. The subsequent adaptation of the virus to its new host through mutation would increase Ka as suggested here in Table 5 and reported for EBOV in humans (Diehl et al., 2016, Urbanowicz et al., 2016).

The magnitude of ΔHa is largely controlled by how good the molecular fit is between the virus GP surface and the Cr surface and depends on hydrogen bonds and salt bridges. As shown in Table 5, changes in ΔHa (ΔΔHa) may be used to predict changes in Ka providing the molecular basis of the GP/Cr interaction is available at atomic detail to enable some assessment of ΔΔHa. The infectivity of the new virus/host combination relative to the “standard” virus/host would be directly proportional to the change in Ka (as demonstrated from results in Table 2). Thus, thermodynamics provides the link between changes in the virus/host identified through sequencing data and the risk of infection. This in theory allows predictions of the infectivity of emerging virus strains or the susceptibility of novel hosts to be assessed on the basis of the effects of changes in sequence data. It is concluded that thermodynamic approaches have a major contribution to make in developing dose-response models for emerging viruses.

Conflict of interest

None declared.

Disclaimer

The views expressed in this paper are those of the author and not necessarily those of any organisations.