Modeling cellular co-infection and reassortment of bluetongue virus in Culicoides midges

Abstract When related segmented RNA viruses co-infect a single cell, viral reassortment can occur, potentially leading to new strains with pandemic potential. One virus capable of reassortment is bluetongue virus (BTV), which causes substantial health impacts in ruminants and is transmitted via Culicoides midges. Because midges can become co-infected by feeding on multiple different host species and remain infected for their entire life span, there is a high potential for reassortment to occur. Once a midge is co-infected, additional barriers must be crossed for a reassortant virus to emerge, such as cellular co-infection and dissemination of reassortant viruses to the salivary glands. We developed three mathematical models of within-midge BTV dynamics of increasing complexity, allowing us to explore the conditions leading to the emergence of reassortant viruses. In confronting the simplest model with published data, we estimate that the average life span of a bluetongue virion in the midge midgut is about 6 h, a key determinant of establishing a successful infection. Examination of the full model, which permits cellular co-infection and reassortment, shows that small differences in fitness of the two infecting strains can have a large impact on the frequency with which reassortant virions are observed. This is consistent with experimental co-infection studies with BTV strains with different relative fitnesses that did not produce reassortant progeny. Our models also highlight several gaps in existing data that would allow us to elucidate these dynamics in more detail, in particular the times it takes the virus to disseminate to different tissues, and measurements of viral load and reassortant frequency at different temperatures.


Introduction
Some RNA viruses have segmented genomes, meaning that their genomes are divided into between two (e.g. Lassa virus) and twelve (e.g. Colorado tick fever virus) distinct RNA segments (McDonald et al. 2016). When two related segmented viruses co-infect a single host cell, virions produced by the cell can consist of segments taken from each of the infecting genotypes, potentially creating a new genotype that is a hybrid of the two infecting strains (McDonald et al. 2016;Lowen 2018). There are currently eleven families of viruses known to have segmented genomes, including the causative agents of medically important diseases of both humans and animals, such as influenza A virus, Rift Valley fever virus (RVFV), and bluetongue virus (BTV). Occasionally, reassortment can lead to the emergence of new strains of viruses that more effectively evade host immune responses, often by changing the surface antigens of an already successful strain while leaving other segments intact (McDonald et al. 2016). These new viruses may then have increased epidemic potential, most famously in the process known as antigenic shift for influenza virus, a major factor in precipitating flu pandemics (Kilbourne 2006;Lewis 2006;Nelson and Holmes 2007;Lycett, Duchatel, and Digard 2019). Similarly, reassortment can also allow virus strains to more effectively invade new regions or new hosts by generating novel combinations of traits (Ma et al. 2016).
For segmented viruses that are transmitted by insect vectors, there may be an increased potential for reassortment to occur. First, there is the potential for reassortment to occur in both the primary hosts and the vector. Secondly, as many vectors feed on different host species (Burkett-Cadena et al. 2008;Martínez-de la Puente, Figuerola, and Soriguer 2015;Ginsberg et al. 2021), there is an increased potential for co-infection with strains that typically infect different hosts. Moreover, once infected, insects typically remain infected for the remainder of their life span (Mellor 1990;Cheng et al. 2016), increasing the potential for co-infection during subsequent blood meals. Such co-infection can occur either by taking a blood meal on a co-infected host or by sequentially taking two or more blood meals on hosts infected with different strains.
Despite the high potential for reassortment, there are still within-vector obstacles to a reassortant virus strain emergingthe vector must be co-infected with two different strains, the two strains must co-infect a single cell, viable reassortant virions must be produced by that cell, and these virions must eventually be present in the salivary gland at sufficient concentration for transmission to occur during a blood meal. The vector immune response provides another obstacle to reassortment occurring. Vectors only have an innate immune system, and this is often characterized in terms of 'barriers' to infection, including a midgut infection barrier (MIB), midgut escape barrier (MEB), dissemination barrier (DB), and salivary gland infection barrier (Mellor 2000). Reassortant viruses or both of their parental strains must successfully navigate each of these barriers and make it to the salivary glands for transmission to occur. Understanding how these barriers interplay with the timing of blood meals and other parameters of virus replication is key to understanding how reassortant viruses could emerge in vectors.
In this paper, we seek to understand the interplay between blood meal timing, relative viral fitness, and the length of the incubation period in determining the amount of reassortment that occurs in a vector. We focus on BTV, a ten-segment virus that infects ruminant hosts, in Culicoides midges. To this end, we develop a series of three mathematical models of withinmidge viral dynamics of increasing complexity, each of which allows us to model a different type of data: infection prevalence, viral load dynamics, and reassortment frequency. The first model describes the MIB, the first barrier to infection a viral infection faces in a vector; the second describes the infection dynamics in a singly-infected midge; and the third describes the dynamics in a co-infected midge and allows for cellular co-infection and reassortment. By examining the behavior of these models and comparing them to available data, we highlight the potential effects that a variable incubation period (e.g. due to varying temperature), the length of the gap between infectious blood meals, and competition between co-infecting viruses can have on the frequency of reassortment. We also propose new experiments that would generate data that improve our understanding of viral dynamics and reassortment in the midge.

Data
For single-blood-meal infection studies, we largely relied on data from Fu et al. (1999). While there are six studies in total that report viral dynamics following a single BTV infection (Foster and Jones 1979;Chandler et al. 1985;Fu et al. 1999;Wittmann, Mellor, and Baylis 2002;Veronesi et al. 2013;Mills et al. 2017;Federici et al. 2019), we chose to focus on Fu et al. due to the fact that they reported viral titers, infection prevalence, and dissemination in detail and that they undertook both oral and intrathoracic infections in parallel. We focused on one study as the viral dynamics are often qualitatively different between studies, for reasons that are unclear (e.g. Foster and Jones observed two distinct proliferation phases, while others only observed one (Foster and Jones 1979;Fu et al. 1999;Mills et al. 2017)). Fu et al. infected Culicoides variipennis midges with BTV-1, both orally and intrathoracically, and kept them at 24 ∘ ± 1 ∘ C. Viral load and infection rate measurements were then made at various intervals up to 2 weeks. In another arm of the experiment, this study used an additional colony of midges that were refractory to oral infection, though we did not use data from that arm of the experiment.
We used data on reassortment during co-infection from a pair of studies: Samal et al. (1987) and el Hussein et al. (1989) ( Supplementary Fig. S1). These studies were conducted by the same group. Both studies blood-fed C. variipennis with BTV Serotypes 10 and 17, keeping the midges at 25 ∘ C. Blood meal viral loads in all experiments were ∼6 log 10 TCID 50 /ml, similar to Fu et al. (1999). They then used electrophoretic analysis of progeny virus from co-infected flies to assess whether isolated plaques were of reassortant viruses or one of the parental strains. The main difference between the studies was that Samal et al only blood-fed midges simultaneously (i.e. with one coinfected blood meal), whereas el Hussein et al. did this but also looked at sequential infection (i.e. with two separate blood meals separated by a period of time). The periods of time tested by el Hussein were 1, 3, 7, 9, and 11 days following the initial blood meal.

MIB model
The literature on viral dynamics in insects generally frames whether or not an infection is established in terms of a sequence of barriers to infection. These can be either physical barriers, such as the barrier between the midgut cavity and the hemocoel, or immunological barriers, such as when cells utilize RNA interference to inhibit virion production. Culicoides midges are typically thought to have an MIB (i.e. a barrier preventing virus particles in the blood meal from infecting midgut cells), a DB (i.e. a barrier preventing an infected midge from gaining a fully disseminated infection), and possibly an MEB (i.e. a barrier to virions escaping the midgut into the hemocoel and the secondary tissues). Other barriers exist in other insects, such as a salivary gland infection barrier and a transovarial infection barrier, but there is no evidence of these in Culicoides infected with BTV. Due to the difficulty of distinguishing an MEB from a DB, we model these two barriers as one. It is worth noting that when midges that were refractory to oral infection were instead intrathoracically infected in the Fu et al. (1999) study, 100 per cent of them were able to transmit through their saliva, compared to 0 per cent when infected orally, highlighting the importance of the midgut infection and escape barriers.
We first focus on the MIB to gain an understanding of the factors governing this barrier. First, we fitted a logistic function to the prevalence data from Fu et al. to estimate the proportion of midges that have an MIB and the timing of when midges that have such a barrier cease to have detectable levels of virus in their midgut. We fitted a functional relationship, between prevalence, P, time since inoculation, t, and three parameters, ∞ , 1/2 , and k (Fig. 1). The parameters include ∞ , the proportion of midges that achieve a disseminated infection, 1/2 , the time when the prevalence is equal to ∞ +1 2 , and k, which controls the steepness of the curve. Fitting this relationship to data from Fu et al., on the basis of least squares using the optim function in R yielded ∞ = 0.346 infected midges per blood-fed midge, 1/2 = 28.9 h, and = 0.835 h −1 .
Next, we estimated the rate of clearance of virions from the midgut, assuming that they are cleared at a constant rate (Fig. 2). We can model this process as follows: where c b is the rate virion clearance and V 0 is the initial viral load. We then assume that when = 1/2 , V b is at the limit of detection, or Λ, because this was the time when half of those midges that did not develop a detectable infection ceased to have detectable  There are several different reasons that a viremic blood meal could fail to establish an infection in the midge: a physical barrier to infection; host defenses decrease viral production (e.g. through RNA interference) or increase viral clearance (e.g. via antimicrobial peptides) such that the reproduction number is less than 1; or all viral particles in the blood meal could be cleared before they infect a cell. In the third case, the probability that all virions are cleared before one manages to infect a cell is given by the following equation: where p MIB is the probability that the MIB is passed given a sufficiently viremic blood meal, c b is the rate of viral clearance in the midgut, m is the per cell rate of viral entry, T m is the number of target cells in the midgut, and V 0 is the initial number of virions in the blood meal. This expression allows us to explore the impact of varying blood meal size on the probability of infection and to distinguish the classical conception of an MIB from the stochastic elimination of the viral population.

Single infection model
We model the infection process and viral dynamics as a system of ordinary differential equations, where variables with a b subscript pertain to the blood meal, those with an m to the midgut, and those with an s to the secondary tissues (Fig. 3). State variables with a double subscript (mm or ss) refer to cells that are co-infected. Allowing for this is necessary to ensure consistency with the co-infection model. In this system, V b virions are initially ingested in the blood meal. To estimate the initial number of virions in the blood meal, we multiply the concentration in TCID 50 /midge by ln(2) to convert this into approximate plaque-forming units (PFU)/midge. This assumes that the number of PFUs per positive tube follows a Poisson distribution and is derived from the expression ( = 0) = − = 0.5, where is the mean number of PFUs per positive tube. We then assume that the number of PFUs is equal to the initial number of virions. Free virions then infect the uninfected target cells in the midgut (T) at rate . Once infected, a cell enters its eclipse phase, which lasts on average, during which time it does not produce new virions. During the eclipse phase, a cell can become co-infected, in this case only with virions of the same genotype. Following the eclipse phase, cells become productively infected, producing new virions at rate p, a fraction d m of which pass the MEB into the secondary tissues. Similarly, for infected secondary tissues, a fraction d s pass back into the midgut. All cells experience mortality at rate and virions at rate c. New uninfected cells are created at rate 0 to keep a fixed population of cells of 0 . Infected cells can die at a different rate than uninfected cells, though in this study we set these rates to be the same as increased mortality in BTV-infected midges has not been documented.
To model the barriers to infection, we simulated the model under three scenarios: (1) with 0, = 0 to represent midges with an MIB; (2) with 0, > 0 and 0, = 0 to represent midges that do not have an MIB but do have a DB; and (3) with both 0, > 0 and 0, > 0 to represent those midges that have neither barrier, resulting in a fully disseminated infection. We then define  (2) or (3) (Fu et al. 1999). The average viral load near the start of the experiment, when all midges are infected, will then be given as follows: where p MIB is the proportion of midges that will be infected, and p DB is the proportion of midges that develop a fully disseminated infection. Similarly, the average viral load near the end of the experiment, when only midges without an MIB are infected, will be given as follows: From these, we can calculate the average viral load among positive midges at time t as the weighted sum where P is the prevalence over time. V pos is then compared to the average measured viral load in Fu et al.
By modeling the barriers in this way, we are implicitly assuming that some subset of midges have an MIB or DB, and the remainder do not, reflecting the fact that the high or low frequencies of these barriers can be obtained through selective breeding (Mellor, Boorman, and Baylis 2000). In practice, the presence of the barrier reflects many underlying biological mechanisms, including the host genotype, the virus genotype, their interaction, the age of the host, and the temperature at which the host is housed (Mills et al. 2017). While it is still unclear the specific mechanisms underlying each of these barriers, but some research has highlighted the importance of cells inhibiting viral production via RNA interference (Mills, Nayduch, and Michel 2015). Our approach can be thought of as being agnostic to the actual underlying mechanism, but treating the presence/absence of each barrier as a binary outcome determined by the specific conditions of the experimental set-up. This binary outcome reflects the way in which infection barriers are often described in the literature (Mellor, Boorman, and Baylis 2000).

Co-infection model
First, we analyzed how the proportion of reassortant viruses depended on the time since the first blood meal and the gap between the blood meals in the studies of Samal et al. and el Hussein et al. (Samal et al. 1987;el Hussein et al. 1989). To do this, we fitted a generalized additive model (GAM) to the combined data set with a logistic link, with the proportion of plaques that were reassortant viruses as the dependent variables, and two independent variables: the gap between blood meals and the time since the first blood meal at which a measurement was taken. As the data on the proportion of plaques that were reassortant were quite noisy, using a GAM allowed us to smooth the relationship between blood-meal timing and reassortment and isolate the most important features. We used four and five knots, respectively, for the gap between blood meals and the time since the first blood meal, based on visual inspection of results obtained with varying numbers of knots.
Next, we extended our compartmental model to accommodate co-infection (Fig. 4, see Supplementary material for equations). We now have three different types of virions, labeled i, j for the two infecting strains of BTV and r for reassortant viruses. We do not distinguish between different segment combinations, grouping all reassortants together for simplicity. In this model, cells can become co-infected with one or two different types of virions, but not three. The second infection can occur at any point during the eclipse phase of a cell, but we assume that once a cell is productively infected, it becomes refractory to further infection due to superinfection exclusion (Zou et al. 2009). Note that cells infected with a particular type of virion can also be reinfected with the same strain during their eclipse phase, but that cells co-infected with the same virion have no difference in virion production than cells singly infected with that virus. Allowing re-infection in this way ensures that the model achieves ecological neutrality in the sense described by Lipsitch et al. (2009); i.e. the dynamics of the ecological variables (the total number of cells infected with one type, the total number infected with two types, etc.) only depends on the ecological variables, not the specific strains involved. The model also achieves population genetic neutrality, in the sense that there is no a priori stable strain equilibrium of strain frequencies and that it instead depends on the initial frequencies ( Supplementary Fig. S2). These two types of neutrality help ensure the model does not have a hidden assumption about the equilibrium level of virus types and, hence, the level of re-assortment.
When a cell is co-infected with two different types of virions, it can produce new virions of either parental type or reassortant virions. We assume that the proportion of produced virions that are of a particular parental type is as follows: where is the proportion of produced virions that are with-like; i.e. for which reassortment does not occur. This term can be thought of as describing the combined effect of viruses occupying different parts of the cell, e.g. due to inclusions within the cell (Desmet, Anguish and Parker 2014;Lowen 2018), and incompatible segment combinations. The term 1− 2 10 describes those produced virions for which reassortment does occur, but each of the ten segments is from the same parent. Similarly, the proportion of produced virions that are reassortants is as follows: This formulation assumes that all segment combinations are equally likely and that there is no preference in production for either of the co-infecting virus types. We also assume that the overall rate of viral production is unaffected by the multiplicity of infections. In the absence of reassortment (i.e. when = 1), the co-infection model shows identical viral load dynamics as the single infection model, irrespective of the initial mix of virus types, provided that both virus types have the same parameters ( Supplementary Fig. S2).
The results of Samal et al. (1987) and El Hussein et al. (1989) (Fig. 5) suggest that when there is a small gap between infectious blood meals, more reassortment is observed. We hypothesize that this might be due to a decreased MEB when multiple blood meals are taken, consistent with the work of Brackney, LaReau, and Smith (2021) and Armstrong et al. (2020), which suggests that micro-perforations in the basal lamina caused by multiple blood meals aid viral dissemination. We model this as an increase in the probability of passing the midgut escape/DB on the logistic scale to ensure that the probability remains less than 1, resulting in the following equation: Hence, if we assume that the parameters governing the infection dynamics do not differ between the two infecting viruses and reassortant viruses, then the co-infection model has only two additional parameters compared to the single-infection model: and inc db .

Parameter estimation
We  predicted proportion of virions that are reassortant, at time t when blood meals were separated by g. We fitted to the GAM prediction as the data are somewhat noisy, and we felt this approach might help the dynamical model match the important features of the relationship between blood meal timing and reassortment, particularly the peak in reassortment when the gap between blood meals is around 3 days. We then summed the log-likelihoods to obtain the overall log-likelihood, where T is the set of time-points used by Fu et al., and G is the set of combinations of time-points and gaps between blood meals used in Samal et al. and el Hussein et al. As there is little prior knowledge of the values of any parameters, we used uniform priors on all parameters, with ranges given in Supplementary Table S1. We used the differential evolution Markov chain Monte Carlo (DEzs MCMC) sampler as implemented in the BayesianTools 0.1.7 package in R 4.0.0 (R Core Team 2018; Hartig, Minunno, and Paul 2019).

The roles of temperature and viral fitness
Wittman et al. showed that temperature has a large impact on the extrinsic incubation period (EIP) of BTV in Culicoides, finding that 1/EIP is linearly related to temperature (Wittmann, Mellor, and Baylis 2002). The same study also showed that temperature had a negligible effect on vector competence. Additionally, there is evidence that the viral dynamics differ between serotypes, and between strains within a serotype (Mills et al. 2017;Kopanke et al. 2021). Both of these forms of heterogeneity are likely to have important impacts on the proportion of virions that are reassortants, so we explored their effects on reassortment. Both EIP and the overall dynamics are a complex combination of all parameters, so we made the following simplifications: (a) To assess the role of EIP, we varied the length of the eclipse phase in both the midgut and secondary tissues. This changes the EIP, without affecting equilibrium viral loads, in line with Wittmann et al., who defined the EIP as the first time following the eclipse phase that viral loads reach 10 2.5 TCID 50 /midge.
(b) To assess the role of differing dynamics between serotypes, we varied the ratio of the viral production rates, p, between the two viral types involved in the co-infection while keeping the geometric mean of these parameters fixed at its fitted value. We kept the production rate of reassortant viruses at the baseline fitted level (so it will be between the two production rates of the parental types) but tried an additional sensitivity analysis where it was set to the maximum of the two parental types. Increasing p increases both the equilibrium viral load of that virus type and its growth rate.
We assessed three targets in this analysis: (1) the proportion of virions that are reassortants on Day 10 when midges were co-infected on Day 0; (2) the proportion of virions that are reassortants on Day 10 when midges were infected with the second virus type 3 days after the first; and (3) the ratio of Target (1) to Target (2), which gives an insight into whether the proportion of virions that are reassortants increases when there is a short gap between infections, as observed in el Hussein et al. (1989). As a complement to (b), we also undertook this analysis for different values of instead of different values of p.
In addition, we also undertook a systematic exploration of parameter space using a Sobol sweep before calculating Saltelli indices for the same three targets. These indices describe the amount of variation in the output that is explained by a particular input: which parameters were varied and their ranges are shown in Supplementary Table S1. We used the sensobol 1.

Midgut infection barrier
A midge feeding on an infectious blood meal could fail to become infected either because it has an MIB or because the infecting virions fail to infect a cell before they are all cleared. At initial viral loads similar to those observed in Fu et al. and a number of other studies (∼10 3 TCID 50 /midge), we found that the overall probability of infection does not substantially differ from the probability that the midge has an MIB, unless the rate of viral entry into cells is very low (Fig. 6A). If the initial viral load is much lower than The relationship between the initial viral load in the blood meal, the rate at which virions enter cells, and the probability of a midgut infection being established (color-scale), assuming that two-thirds of midges have an MIB. In both panels, vertical lines indicate the baseline value of that parameter, and circles indicate the best fitting parameters when fitting the model. that observed in Fu et al., then the probability of infection could be lower than that suggested by the MIB alone at more moderate rates of viral entry ( ≤ 1 new infected cells per hour per virion) (Fig. 6B). However, at values of the rate of viral entry similar to those that we found when we calibrated the model ( ≃ 5new infected cells per hour per virion, i.e. it takes 0.2 h on average for a virion to infect a cell), we found that even very low infecting doses would lead to an infection in the absence of an MIB (Fig. 6B). In other words, there is little dose-response effect on the prevalence of infection.

Single-infection model
The single-infection model was able to reproduce several key features of within-midge viral dynamics following both oral and intrathoracic infection. First, the eclipse phase following an infectious blood meal (i.e. the time until the average viral load of positive midges has reached a level higher than the infecting dose) lasted on average about 100 h, and the viral load reached its minimum value after around 30 h (Fig. 7A). Our model highlights two differing mechanisms causing this eclipse phase: (1) when a blood meal is imbibed, it takes time before cells are infected and until those cells produce new virions (see black lines in Fig. 7) and (2) early in the experiment, the average viral load of positive midges consists of both infected midges and midges that are not infected but still have virus remaining from the blood meal, so averaging these decreases the measured viral load (compare the black, red, and green lines in Fig. 7A). Second, our model is able to capture the observation of Fu et al. that some midges become infected but do not develop a fully disseminated infection, thereby remaining positive for the duration of the experiment but with a steadily declining viral load (gray line in Fig. 7A). Our model is also able to reproduce the viral dynamics in those midges that do not have any barriers to infection, although in this case the equilibrium viral load is slightly underestimated (Fig. 7B). This suggests that this population of midges can be distinguished by more than just an absence of midgut and DBs; for instance they might also have higher rates of viral production in their secondary tissues.

Co-infection model
We jointly calibrated the model to single-infection dynamics and the proportion of virions that are reassortant during co-infection, assuming identical parameters for the two infecting strains in the latter case. Following this, we estimated that 17.5 per cent (95 per cent credible interval (CrI): 14.4-18.2 per cent) of virions produced by a co-infected cell are reassortant viruses, and the remainder is one of the parental virus types. While there are few empirical estimates of the frequency of reassortment from single cells, it is worth noting that the overall frequency of reassortant viruses in co-infected hosts can be quite high, such as those seen in Samal et al. (1987) andel Hussein et al. (1989), who observed that 42 per cent and 48 per cent, respectively, of viruses isolated from simultaneously co-infected midges were reassortant. High reassortant frequencies have also been observed in other systems, such as influenza in guinea pigs, where between 46 per cent and 86 per cent of isolated viruses were reassortant (Marshall et al. 2013). High levels of reassortment have also been observed in cell culture for influenza; Hockman et al. observed that even when only 20 per cent of cells were co-infected, around 70 per cent of isolated viruses were reassortant (Hockman et al. 2020). We also estimated that when the midge took two blood meals rather than just one (or in other words, when the two infections were separated in time), the probability that the midge had a DB substantially decreased from 88 per cent to 2.4 per cent [95 per cent CrI: 2.0-2.4 per cent].
The dynamics of the total viral load when we modeled coinfection with a single blood meal on Day 0 was identical to the dynamics of a single infection (Figs 7 and 8, Supplementary  Fig. S3). This is to be expected, as both viruses had the same parameters and the total number of virions per blood meal was fixed, but further verifies that the model behaves in a neutral manner with respect to co-infection. When the second infection occurred later than the first, we observed higher equilibrium loads and an increase in the growth rate around 24 h after the second blood meal (Fig. 8, particularly for 120-h gap or more). These effects are caused by the decrease in the probability of a DB following the second blood meal, perhaps due to micro-perforations in the midgut basal lamina (Brackney, LaReau, and Smith 2021). The reassortant viral load reached 100 TCID 50 /midge around 24 h after the second blood meal, though this length of time increased as the gap between the blood meal increased (Fig. 8). Notably, it takes much longer (around 48 h) to reach 100 TCID 50 /midge when midges are fed a single co-infecting blood meal, in part due to the higher DB, and in part because there are no pre-existing eclipsing cells as there are when the infections are separated in time.
In the empirical data, the proportion of virions that were found to be reassortants varied with the gap between blood meals, showing an initial increase to around 70 per cent reassortants as the gap increased from 0 to ∼3 days, before declining to zero at longer gaps (Fig. 5B). There was not a strong relationship between the proportion reassortant and the time at which the measurement was taken, with the GAM fit suggesting the possibility of a slight decrease for later measurements (Fig. 5A). Our dynamical model was able to capture the approximate magnitude of the proportion reassortant for midges that were co-infected with a single blood meal on Day 0 (Fig. 9), and also predicted little or no reassortment for longer gaps. However, the model could not capture the increase in the proportion reassortant as the gap between infections increased to ∼3 days. Although this type of pattern is within the range of the model's behavior, parameterizations that allowed for this increase were incommensurate with the overall viral load dynamics found by Fu et al. Additionally, our model always predicts increasing proportions of reassortant viruses with later measurements, something which was not found by el Hussein et al. and Samal et al.

The role of temperature and viral fitness
To explore the effect of differing viral fitness and a variable incubation period, we varied the eclipse phase and the ratio of the viral production rates between the two viruses, keeping all other parameters fixed. When the two virus types involved with the co-infection had different production rates, we only observed an appreciable amount of reassortment when the production rates were similar (Fig. 10A, Supplementary S4). When production rates differed by a factor of two or more, we found that less than 5 per cent of virions were reassortants (Fig. 10C). This was the case whether it was the first infection or the second infection that had a higher production rate. Most reassortment occurred when the second virus had a slightly higher production rate and the eclipse phase was short (Fig. 10A, C-see the dark red region above the dashed line). When we reduced the eclipse phase, we saw a larger proportion of reassortant viruses by Day 10 than with the calibrated value of the eclipse phase (Fig. 10A, B, dotted line). The effect of the eclipse phase on the proportion of reassortment observed depends on the relative viral production rates. For example, if the second virus has a lower production rate, then a longer eclipse phase may lead to more reassortment (Fig. 10, bottom right  region). When the ratio in the viral production rates was greater than ∼2 (or less than ∼0.5), we were able to recover the increase in the proportion of reassortants observed when the two infections were separated by 3 days compared to when they both occurred on Day 0 (Fig. 11, blue region), as indicated by the empirical data (Fig. 5B). However, at most of these parameter values, we also saw very little reassortment in general (Fig. 10), incommensurate with the data from el Hussein et al. The only region of parameter space where there were both substantial reassortment and more reassortment when infections were separated by 3 days was when the eclipse phase was long and the first virus had slightly higher production rates (Figs 10 and 11, near the top of the lower right section). Varying the relative rate of effective contact between virions and target cells ( ), instead of the relative production rate, yielded almost identical patterns of reassortment ( Supplementary Figs S5-S7).
In the co-infection model with the calibrated value of the eclipse phase and identical growth parameters for each virus, a variance-based sensitivity analysis suggested that the parameter describing the proportion of produced virions that are with-like ( ) is most important in determining the proportion of virions that are reassortants overall (Supplementary Figs S8 and S9). This parameter explains 79 per cent of the variance in the proportion of virions that are reassortants when both infections occur on Day 0, and 39 per cent when the second infection occurs on Day 3. Other important parameters are the rate at which cells die ( ) and the rate at which virions are cleared in the secondary tissues (c s ). For co-infections, m explains 5.7 per cent of the variance, s explains 1.2 per cent, and c s explains 0.13 per cent, while for infections separated by 3 days, m explains 15.6 per cent of the variance, Figure 9. The proportion of virions that are reassortants at different gaps between blood meals (x-axis) and different times since the first infection (hue of lines). The left hand plot show predictions from a GAM fitted to data from Samal et al. and el Hussein et al, and the right hand the predictions from the dynamical model. Figure 10. The proportion of reassortant viruses observed 10 days after the initial infection (color-scale) when the second blood meal occurs 3 days after the initial infection, for different lengths of eclipse phase (x-axis) and for different viral production rates between the two infections (y-axis). When the fold-change in p is less than 1, the first infecting virus has a higher production rate than the second, and when it is greater than 1, the second virus has a higher production rate. Dashed lines represent baseline values. Note that the eclipse phase in the secondary tissues was also varied proportionally by the same amount from baseline as the eclipse phase in the midgut (not shown). Note that the left-hand plot has an x-axis that is on a log scale-reassortment decreases very quickly as the difference in viral fitness increases. account for 13 per cent of the variance for co-infections and 32 per cent for infections separated by 3 days. The interaction between m and in particular accounts for 6.7 per cent of the overall variance for co-infections and 7.8 per cent when infections are separated by 3 days.

Discussion
In this paper, we used three models of increasing complexity to elucidate the viral dynamics and emergence of reassortment in an insect vector. The first model allowed us to estimate that virions in the midgut have a life span of around 6 h on average. The second Figure 11. The ratio of the proportion of reassortant viruses observed 10 days after the initial infection when both blood meals occur on Day 0 to the same proportion when the second blood meal occurs on Day 3 (color-scale), for different lengths of eclipse phase (x-axis) and for different viral production rates between the two infections (y-axis). When this ratio is below 1, it indicates an increase in the proportion of viruses that are reassortant, as observed in el Hussein et al. When the fold-change in p is less than 1, the first infecting virus has a higher production rate than the second, and when it is greater than 1 the second virus has a higher production rate. Dashed lines represent baseline values. Note that the eclipse phase in the secondary tissues was also varied proportionally by the same amount from baseline as the eclipse phase in the midgut (not shown). Note also that the left-hand plot has an x-axis that is on a log scale. model is able to capture the temporal pattern of viral load within the midge but is poorly identifiable, suggesting that more data are needed to constrain elements of the model. Experiments to obtain such data are proposed below. Our final model captures the broad patterns in the frequency of reassortment at different times since the first infection and gaps between the two infections, but when virus parameters are the same across virus genotypes, it cannot capture the specific pattern of increased reassortant frequency at small gaps between blood meals. We were able to recover this pattern when viral production rates of the infecting strains differed, however. The final model also showed that even small relative differences in viral production rates can cause large differences in the frequency of reassortment.
Our results highlight the importance of the relative fitness of the infecting strains and blood meal timing. As previously mentioned, we could only reconstruct the increase in reassortment when blood meals were separated by 3 days by allowing the strains to have different production rates. We initially tried to reconstruct this increase by allowing second blood meals to have a higher probability of passing the MIB, something previously observed in Aedes aegypti mosquitoes (Armstrong et al. 2020), but surprisingly this could not recreate the observed increase in reassortment when viral production rates were identical. Prior studies of BTV co-infection have found conflicting results regarding reassortment frequency; while both Samal et al. and el Hussein et al. found substantial levels of reassortment when coinfecting C. variipennis with BTV-10 and BTV-17, a more recent study by Kopanke et al. observed no reassortment when infecting Culicoides sonorensis with BTV-2 and BTV-10 (Samal et al. 1987;El Hussein et al. 1989;Kopanke et al. 2021). A possible explanation for this discrepancy is suggested by our observation that when viral production rates varied by more than a modest amount, very little reassortment was observed. These experimental results and our modeling results suggest that whether reassortment is observed, and how much, may be a consequence of the delicate balance between relative viral fitness and the timing of infections. Relative viral fitness will in turn depend on the combination of vector and virus. It would be interesting to experimentally infect midges with a less fit virus followed by a fitter virus later to observe whether these effects are observed in practice. Our modeling results also corroborate the implications of studies by Samal et al. and el Hussein et al. that when infections are separated by more than around parasites in Anopheles mosquitoes (Teboh-Ewungkem and Yuster 2010; Childs andProsper 2017, 2020;Stopard, Churcher, and Lambert 2021). All four studies modeled transitions between the life stages of malaria parasites in the mosquito, or a subset of them, and did not explicitly model infection of mosquito cells. Prior models of viral dynamics in vectors have focused on Zika virus dynamics in Aedes mosquitoes (Lequime et al. 2020;Tuncer and Martcheva 2021). Both of these models modeled virus dynamics phenomenologically as logistic growth, though Tuncer et al. included a semi-mechanistic component describing both midgut and secondary tissues and transport between them. Like the malaria models but unlike our model, neither study directly modeled the infection of mosquito cells. These modeling approaches, though suitable for the questions asked in those studies, are less conducive to modeling cellular co-infection and reassortment, as was our aim here. At least two non-vector models of within-host reassortment have been developed for influenza (Leonard et al. 2017;Koelle et al. 2019). Koelle et al. proposed a parsimonious compartmental model of cellular co-infection and reassortment. In their model, all co-infected cells produce only reassortant viruses, meaning that eventually all virions are reassortant, which is incommensurate with the data on reassortment of BTV in midges.
Our study has at least four limitations. First, our model did not account for stochastic effects, which might be important if there is a bottleneck following one or more of the barriers that lead to low viral population sizes. Second, we did not allow for inter-midge variability in viral dynamics, which would allow us to capture the distribution in viral loads and reassortant frequencies. Third, we made the simplifying assumption that reassortants had a production rate equal to the geometric mean of their parental viruses. If all reassortant viruses are much less fit than their parental strains, we would obviously expect much lower reassortant frequency, and this might provide an alternative explanation to the findings of Kopanke et al. (2021). Fourth, we implicitly assume that all segment combinations are equally likely, with the exception of the parental combinations whose frequency is determined by the with-like parameter, . However, as we do not explicitly track different segment combinations, our reassortment compartment could be thought of as a frequency-weighted sum of all reassortants. The latter two limitations are both partially addressed by : both reduced reassortant fitness and a reduced set of possible segment combinations could be interpreted as an increase in .
Our modeling suggests several in vivo experiments that could help improve future modeling efforts and in turn our understanding of BTV dynamics in Culicoides. First, it would be useful to know the distribution of the timing of dissemination to different body parts, especially the midgut cells, the fat body cells, the ganglia, and the salivary glands (cell types that BTV is commonly observed in), including the time virus is first observed or crosses some threshold in that body part or the viral load in different tissue types over time. The much smaller size of the midge has made this more difficult for BTV in Culicoides, though Fu et al. did record this for a subset of midges without measuring the full distribution of timings (Fu et al. 1999). Second, the dose-response effect on viral load, infection prevalence, and reassortant frequency could help inform the model's structure. For example, if we were to compare the dose-response effect in the prevalence of infection in midges, then we could constrain the rate of cellular entry by virions in the midgut (Fig. 6B). Third, elucidation of the mechanisms governing the barriers to infection (i.e. the vector immune response) would allow a more mechanistic treatment of the barriers. Our model currently treats the barriers phenomenologically, without directly incorporating them into the model structure. Fourth, to account for inter-midge variability and include stochasticity, it would be useful to know the distribution of viral loads across midges over time, and how this differs between those midges that develop a fully disseminated infection and those that do not. This has been characterized for dengue virus in A. aegypti (Novelo et al. 2019). Finally, while we explored one possible effect of different temperatures via examination of a varying eclipse phase, a thorough description of viral load trajectories and reassortant frequencies at different temperatures would enable our more systematic inclusion of temperature dependence in viral dynamics.
Our approach and findings are most immediately relevant to other vector-borne viruses with segmented genomes. One such example is RVFV (Sall et al. 1999), for which reassortment in the mosquito vector has been shown to both occur frequently and be readily transmitted in laboratory conditions (Turell et al. 1990). When translating our approach to other vector-virus combinations, the specific barriers will certainly differ, and this would need to be incorporated into the model. Our findings around relative viral fitness and infection timing also have implications for other segmented virus systems, such as influenza viruses and rotaviruses in mammals. While such hosts, unlike vectors, typically have acute infections, this may only serve to increase the importance of relative fitness and timing. When translating our approach to these systems, the inclusion of the barriers to infection would no longer be needed, but it would become necessary to model cellular death and perhaps to more explicitly model the immunological response of the host.
Our study adds a mechanistic description of cellular coinfection and reassortment in an insect vector. Our results give a quantitative understanding of the typical time a virion is resident in the midgut and highlight how the amount of BTV reassortment observed during midge co-infection will depend on the interplay between the relative fitness of the two infecting strains and the timing of the respective blood meals. The relatively simple approach we employ would be straightforward to adapt to other vector-virus combinations with different barriers to infection or with barriers of different strengths. It would also be straightforward to incorporate new data types, and we suggest a range of experiments to generate such data. Such experimental work would allow us to better constrain the model's parameters and relate them to different temperature regimes. In addition, future modeling studies could focus on coupling this model with an epidemiological model to try and understand the epidemic potential of reassortment in the midge or applying this work to other vector-borne segment viruses.

Data availability
All data and code are available at https://github.com/scavany/ bluetongue_within_midge.

Supplementary data
Supplementary data are available at Virus Evolution online.
Conflict of interest: None declared.