Getting the most out of maths: How to coordinate mathematical modelling research to support a pandemic, lessons learnt from three initiatives that were part of the COVID-19 response in the UK.

In March 2020 mathematics became a key part of the scientific advice to the UK government on the pandemic response to COVID-19. Mathematical and statistical modelling provided critical information on the spread of the virus and the potential impact of different interventions. The unprecedented scale of the challenge led the epidemiological modelling community in the UK to be pushed to its limits. At the same time, mathematical modellers across the country were keen to use their knowledge and skills to support the COVID-19 modelling effort. However, this sudden great interest in epidemiological modelling needed to be coordinated to provide much-needed support, and to limit the burden on epidemiological modellers already very stretched for time. In this paper we describe three initiatives set up in the UK in spring 2020 to coordinate the mathematical sciences research community in supporting mathematical modelling of COVID-19. Each initiative had different primary aims and worked to maximise synergies between the various projects. We reflect on the lessons learnt, highlighting the key roles of pre-existing research collaborations and focal centres of coordination in contributing to the success of these initiatives. We conclude with recommendations about important ways in which the scientific research community could be better prepared for future pandemics. This manuscript was submitted as part of a theme issue on "Modelling COVID-19 and Preparedness for Future Pandemics".

conditions for Turing and wave instabilities in reaction-diffusion systems.

Necessary and sufficient conditions are provided for a diffusion-driven instability of a stable equilibrium of a reaction-diffusion system with n components and diagonal diffusion matrix. These can be either Turing or wave instabilities. Known necessary and sufficient conditions are reproduced for there to exist diffusion rates that cause a Turing bifurcation of a stable homogeneous state in the absence of diffusion. The method of proof here though, which is based on study of dispersion relations in the contrasting limits in which the wavenumber tends to zero and to [Formula: see text], gives a constructive method for choosing diffusion constants. The results are illustrated on a 3-component FitzHugh-Nagumo-like model proposed to study excitable wavetrains, and for two different coupled Brusselator systems with 4-components.

Unified framework for localized patterns in reaction-diffusion systems; the Gray-Scott and Gierer-Meinhardt cases.

A recent study of canonical activator-inhibitor Schnakenberg-like models posed on an infinite line is extended to include models, such as Gray-Scott, with bistability of homogeneous equilibria. A homotopy is studied that takes a Schnakenberg-like glycolysis model to the Gray-Scott model. Numerical continuation is used to understand the complete sequence of transitions to two-parameter bifurcation diagrams within the localized pattern parameter regime as the homotopy parameter varies. Several distinct codimension-two bifurcations are discovered including cusp and quadruple zero points for homogeneous steady states, a degenerate heteroclinic connection and a change in connectedness of the homoclinic snaking structure. The analysis is repeated for the Gierer-Meinhardt system, which lies outside the canonical framework. Similar transitions are found under homotopy between bifurcation diagrams for the case where there is a constant feed in the active field, to it being in the inactive field. Wider implications of the results are discussed for other pattern-formation systems arising as models of natural phenomena. This article is part of the theme issue 'Recent progress and open frontiers in Turing's theory of morphogenesis'.

Use of machine learning to analyse routinely collected intensive care unit data: a systematic review.

BACKGROUND: Intensive care units (ICUs) face financial, bed management, and staffing constraints. Detailed data covering all aspects of patients' journeys into and through intensive care are now collected and stored in electronic health records: machine learning has been used to analyse such data in order to provide decision support to clinicians. METHODS: Systematic review of the applications of machine learning to routinely collected ICU data. Web of Science and MEDLINE databases were searched to identify candidate articles: those on image processing were excluded. The study aim, the type of machine learning used, the size of dataset analysed, whether and how the model was validated, and measures of predictive accuracy were extracted. RESULTS: Of 2450 papers identified, 258 fulfilled eligibility criteria. The most common study aims were predicting complications (77 papers [29.8% of studies]), predicting mortality (70 [27.1%]), improving prognostic models (43 [16.7%]), and classifying sub-populations (29 [11.2%]). Median sample size was 488 (IQR 108-4099): 41 studies analysed data on > 10,000 patients. Analyses focused on 169 (65.5%) papers that used machine learning to predict complications, mortality, length of stay, or improvement of health. Predictions were validated in 161 (95.2%) of these studies: the area under the ROC curve (AUC) was reported by 97 (60.2%) but only 10 (6.2%) validated predictions using independent data. The median AUC was 0.83 in studies of 1000-10,000 patients, rising to 0.94 in studies of > 100,000 patients. The most common machine learning methods were neural networks (72 studies [42.6%]), support vector machines (40 [23.7%]), and classification/decision trees (34 [20.1%]). Since 2015 (125 studies [48.4%]), the most common methods were support vector machines (37 studies [29.6%]) and random forests (29 [23.2%]). CONCLUSIONS: The rate of publication of studies using machine learning to analyse routinely collected ICU data is increasing rapidly. The sample sizes used in many published studies are too small to exploit the potential of these methods. Methodological and reporting guidelines are needed, particularly with regard to the choice of method and validation of predictions, to increase confidence in reported findings and aid in translating findings towards routine use in clinical practice.

Revealing instabilities in a generalized triadic supply network: A bifurcation analysis.

Supply networks are exposed to instabilities and thus a high level of risk. To mitigate this risk, it is necessary to understand how instabilities are formed in supply networks. In this paper, we focus on instabilities in inventory dynamics that develop due to the topology of the supply network. To be able to capture these topology-induced instabilities, we use a method called generalized modeling, a minimally specified modeling approach adopted from ecology. This method maps the functional dependencies of production rates on the inventory levels of different parts and products, which are imposed by the network topology, to a set of elasticity parameters. We perform a bifurcation analysis to investigate how these elasticities affect the stability. First, we show that dyads and serial supply chains are immune to topology-induced instabilities. In contrast, in a simple triadic network, where a supplier acts as both a first and a second tier supplier, we can identify instabilities that emerge from saddle-node, Hopf, and global homoclinic bifurcations. These bifurcations lead to different types of dynamical behavior, including exponential convergence to and divergence from a steady state, temporary oscillations around a steady state, and co-existence of different types of dynamics, depending on initial conditions. Finally, we discuss managerial implications of the results.

Respiratory modulated sympathetic activity: a putative mechanism for developing vascular resistance?

KEY POINTS: Sympathetic activity exhibits respiratory modulation that is amplified in hypertensive rats. Respiratory modulated sympathetic activity produces greater changes in vascular resistance than tonic stimulation of the same stimulus magnitude in normotensive but not hypertensive rats. Mathematical modelling demonstrates that respiratory modulated sympathetic activity may fail to produce greater vascular resistance changes in hypertensive rats because the system is saturated as a consequence of a dysfunctional noradrenaline reuptake mechanism. Respiratory modulated sympathetic activity is an efficient mechanism to raise vascular resistance promptly, corroborating its involvement in the ontogenesis of hypertension. ABSTRACT: Sympathetic nerve activity (SNA) exhibits respiratory modulation. This component of SNA is important - being recruited under cardiorespiratory reflex conditions and elevated in the spontaneously hypertensive (SH) rat - and yet the exact influence of this modulation on vascular tone is not understood, even in normotensive conditions. We constructed a mathematical model of the sympathetic innervation of an arteriole, and used it to test the hypothesis that respiratory modulation of SNA preferentially increases vasoconstriction compared to a frequency-matched tonic pattern. Simulations supported the hypothesis, where respiratory modulated increases in vasoconstriction were mediated by a noradrenergic mechanism. These predictions were tested in vivo in adult Wistar rats. Stimulation of the sympathetic chain (L3) with respiratory modulated bursting patterns, revealed that bursting increases vascular resistance (VR) more than tonic stimulation (57.8 ± 3.3% vs. 44.8 ± 4.2%; P < 0.001; n = 8). The onset of the VR response was also quicker for bursting stimulation (rise time constant = 1.98 ± 0.09 s vs. 2.35 ± 0.20 s; P < 0.01). In adult SH rats (n = 8), the VR response to bursting (44.6 ± 3.9%) was not different to tonic (37.4 ± 3.5%; P = 0.57). Using both mathematical modelling and in vivo techniques, we have shown that VR depends critically on respiratory modulation and revealed that this pattern dependency in Wistar rats is due to a noradrenergic mechanism. This respiratory component may therefore contribute to the ontogenesis of hypertension in the pre-hypertensive SH rat - raising VR and driving vascular remodelling. Why adult SH rats do not exhibit a pattern-dependent response is not known, but further modelling revealed that this may be due to dysfunctional noradrenaline reuptake.

Modelling the vascular response to sympathetic postganglionic nerve activity.

This paper explores the influence of burst properties of the sympathetic nervous system on arterial contractility. Specifically, a mathematical model is constructed of the pathway from action potential generation in a sympathetic postganglionic neurone to contraction of an arterial smooth muscle cell. The differential equation model is a synthesis of models of the individual physiological processes, and is shown to be consistent with physiological data. The model is found to be unresponsive to tonic (regular) stimulation at typical frequencies recorded in sympathetic efferents. However, when stimulated at the same average frequency, but with repetitive respiratory-modulated burst patterns, it produces marked contractions. Moreover, the contractile force produced is found to be highly dependent on the number of spikes in each burst. In particular, when the model is driven by preganglionic spike trains recorded from wild-type and spontaneously hypertensive rats (which have increased spiking during each burst) the contractile force was found to be 10-fold greater in the hypertensive case. An explanation is provided in terms of the summative increased release of noradrenaline. Furthermore, the results suggest the marked effect that hypertensive spike trains had on smooth muscle cell tone can provide a significant contribution to the pathology of hypertension.

The use of normal forms for analysing nonlinear mechanical vibrations.

A historical introduction is given of the theory of normal forms for simplifying nonlinear dynamical systems close to resonances or bifurcation points. The specific focus is on mechanical vibration problems, described by finite degree-of-freedom second-order-in-time differential equations. A recent variant of the normal form method, that respects the specific structure of such models, is recalled. It is shown how this method can be placed within the context of the general theory of normal forms provided the damping and forcing terms are treated as unfolding parameters. The approach is contrasted to the alternative theory of nonlinear normal modes (NNMs) which is argued to be problematic in the presence of damping. The efficacy of the normal form method is illustrated on a model of the vibration of a taut cable, which is geometrically nonlinear. It is shown how the method is able to accurately predict NNM shapes and their bifurcations.

Increased intrinsic excitability of muscle vasoconstrictor preganglionic neurons may contribute to the elevated sympathetic activity in hypertensive rats.

Hypertension is associated with pathologically increased sympathetic drive to the vasculature. This has been attributed to increased excitatory drive to sympathetic preganglionic neurons (SPN) from brainstem cardiovascular control centers. However, there is also evidence supporting increased intrinsic excitability of SPN. To test this hypothesis, we made whole cell recordings of muscle vasoconstrictor-like (MVClike) SPN in the working-heart brainstem preparation of spontaneously hypertensive (SH) and normotensive Wistar-Kyoto (WKY) rats. The MVClike SPN have a higher spontaneous firing frequency in the SH rat (3.85 ± 0.4 vs. 2.44 ± 0.4 Hz in WKY; P = 0.011) with greater respiratory modulation of their activity. The action potentials of SH SPN had smaller, shorter afterhyperpolarizations (AHPs) and showed diminished transient rectification indicating suppression of an A-type potassium conductance (IA). We developed mathematical models of the SPN to establish if changes in their intrinsic properties in SH rats could account for their altered firing. Reduction of the maximal conductance density of IA by 15-30% changed the excitability and output of the model from the WKY to a SH profile, with increased firing frequency, amplified respiratory modulation, and smaller AHPs. This change in output is predominantly a consequence of altered synaptic integration. Consistent with these in silico predictions, we found that intrathecal 4-aminopyridine (4-AP) increased sympathetic nerve activity, elevated perfusion pressure, and augmented Traube-Hering waves. Our findings indicate that IA acts as a powerful filter on incoming synaptic drive to SPN and that its diminution in the SH rat is potentially sufficient to account for the increased sympathetic output underlying hypertension.

Time-frequency composition of mosquito flight tones obtained using Hilbert spectral analysis.

Techniques for estimating temporal variation in the frequency content of acoustic tones based on short-time fast Fourier transforms are fundamentally limited by an inherent time-frequency trade-off. This paper presents an alternative methodology, based on Hilbert spectral analysis, which is not affected by this weakness, and applies it to the accurate estimation of mosquito wing beat frequencies. Mosquitoes are known to communicate with one another via the sounds generated by their flapping wings. Active frequency modulation between pairs of mosquitoes is thought to take place as a precursor to courtship. Studying the acoustically-based interactions of mosquitoes therefore relies on an accurate representation of flight frequency as a time-evolving property, yet conventional Fourier spectrograms are unable to capture the rapid modulations in frequency that mosquito flight tones exhibit. The algorithms introduced in this paper are able to automatically detect and extract fully temporally resolved frequency information from audio recordings. Application of the technique to experimental recordings of single tethered mosquitoes in flight reveals corroboration with previous reported findings. The advantages of the method for animal communication studies are discussed, with particular attention given to its potential utility for studying pairwise mosquito interactions.

Comparison of nonlinear mammalian cochlear-partition models.

Various simple mathematical models of the dynamics of the organ of Corti in the mammalian cochlea are analyzed and their dynamics compared. The specific models considered are phenomenological Hopf and cusp normal forms, a recently proposed description combining active hair-bundle motility and somatic motility, a reduction thereof, and finally a model highlighting the importance of the coupling between the nonlinear transduction current and somatic motility. It is found that for certain models precise tuning to any bifurcation is not necessary and that a compressively nonlinear response over a range similar to experimental observations and that the normal form of the Hopf bifurcation is not the only description that reproduces compression and tuning similar to experiment.

A new phase model of the spatiotemporal relationships between three circadian oscillators in the brainstem.

Analysis of ex vivo Per2 bioluminescent rhythm previously recorded in the mouse dorsal vagal complex reveals a characteristic phase relationship between three distinct circadian oscillators. These signals represent core clock gene expression in the area postrema (AP), the nucleus of the solitary tract (NTS) and the ependymal cells surrounding the 4th ventricle (4Vep). Initially, the data suggests a consistent phasing in which the AP peaks first, followed shortly by the NTS, with the 4Vep peaking 8-9 h later. Wavelet analysis reveals that this pattern is not consistently maintained throughout a recording, however, the phase dynamics strongly imply that oscillator interactions are present. A simple phase model of the three oscillators is developed and it suggests that realistic phase dynamics occur between three model oscillators with coupling close to a synchronisation transition. The coupling topology suggests that the AP bidirectionally communicates phase information to the NTS and the 4Vep to synchronise the three structures. A comparison of the model with previous experimental manipulations demonstrates its feasibility to explain DVC circadian phasing. Finally, we show that simulating steadily decaying coupling improves the model's ability to capture experimental phase dynamics.

Modelling a dynamic magneto-agglutination bioassay.

The process of developing an end-to-end model of a magneto-immunoassay is described, simulating the agglutination effect due to the specific binding of bacteria to paramagnetic particles. After establishing the properties of the dose-specific agglutination through direct imaging, a microfluidic assay was used to demonstrate changes in the magnetophoretic transport dynamics of agglutinated clusters via transient inductive magentometer measurements. End-to-end mathematical modelling is used to establish the physical processes underlying the assay. First, a modified form of Becker-Döring nucleation kinetic equations is used to establish a relationship between analyte dose and average cluster size. Next, Stokes flow equations are used to establish a relationship between cluster size and speed of motion within the fluid chamber. This predicts a cluster-size dynamic profile of concentration of PMPs versus time when the magnetic field is switched between the two actuated magnets. Finally, inductive modelling is carried out to predict the response of the magnetometer circuit in response to the dynamics of magnetic clusters. The predictions of this model are shown to agree well with the results of experiments, and to predict the shape of the dose-response curve.

Process modelling of NHS cardiovascular waiting lists in response to the COVID-19 pandemic.

OBJECTIVE: To model the referral, diagnostic and treatment pathway for cardiovascular disease (CVD) in the English National Health Service (NHS) to provide commissioners and managers with a methodology to optimise patient flow and reduce waiting lists. STUDY DESIGN: A systems dynamics approach modelling the CVD healthcare system in England. The model is designed to capture current and predict future states of waiting lists. SETTING: Routinely collected, publicly available data streams of primary and secondary care, sourced from NHS Digital, NHS England, the Office of National Statistics and StatsWales. DATA COLLECTION AND EXTRACTION METHODS: The data used to train and validate the model were routinely collected and publicly available data. It was extracted and implemented in the model using the PySD package in python. RESULTS: NHS cardiovascular waiting lists in England have increased by over 40% compared with pre- COVID-19 levels. The rise in waiting lists was primarily due to restrictions in referrals from primary care, creating a bottleneck postpandemic. Predictive models show increasing point capacities within the system may paradoxically worsen downstream flow. While there is no simple rate-limiting step, the intervention that would most improve patient flow would be to increase consultant outpatient appointments. CONCLUSIONS: The increase in NHS CVD waiting lists in England can be captured using a systems dynamics approach, as can the future state of waiting lists in the presence of further shocks/interventions. It is important for those planning services to use such a systems-oriented approach because the feed-forward and feedback nature of patient flow through referral, diagnostics and treatment leads to counterintuitive effects of interventions designed to reduce waiting lists.

The benefits of peer transparency in safe workplace operation post pandemic lockdown.

The benefits of different levels of engagement with test, trace and isolate procedures are investigated for a pandemic in which there is little population immunity, in terms of productivity and public health. Simple mathematical modelling is used in the context of a single, relatively closed workplace such as a factory or back-office where, in normal operation, each worker has lengthy interactions with a fixed set of colleagues. A discrete-time SEIR model on a fixed interaction graph is simulated with parameters that are motivated by the recent COVID-19 pandemic in the UK during a post-peak phase, including a small risk of viral infection from outside the working environment. Two kinds of worker are assumed, transparents who regularly test, share their results with colleagues and isolate as soon as a contact tests positive for the disease, and opaques who do none of these. Moreover, the simulations are constructed as a 'playable model' in which the transparency level, disease parameters and mean interaction degree can be varied by the user. The model is also analysed in the continuum limit. All simulations point to the double benefit of transparency in both maximizing productivity and minimizing overall infection rates. Based on these findings, public policy implications are discussed for how to incentivise this mutually beneficial behaviour in different kinds of workplace, and simple recommendations are made.

Conveyance of texture signals along a rat whisker.

Neuronal activities underlying a percept are constrained by the physics of sensory signals. In the tactile sense such constraints are frictional stick-slip events, occurring, amongst other vibrotactile features, when tactile sensors are in contact with objects. We reveal new biomechanical phenomena about the transmission of these microNewton forces at the tip of a rat's whisker, where they occur, to the base where they engage primary afferents. Using high resolution videography and accurate measurement of axial and normal forces at the follicle, we show that the conical and curved rat whisker acts as a sign-converting amplification filter for moment to robustly engage primary afferents. Furthermore, we present a model based on geometrically nonlinear Cosserat rod theory and a friction model that recreates the observed whole-beam whisker dynamics. The model quantifies the relation between kinematics (positions and velocities) and dynamic variables (forces and moments). Thus, only videographic assessment of acceleration is required to estimate forces and moments measured by the primary afferents. Our study highlights how sensory systems deal with complex physical constraints of perceptual targets and sensors.

Emergence of the London Millennium Bridge instability without synchronisation.

The pedestrian-induced instability of the London Millennium Bridge is a widely used example of Kuramoto synchronisation. Yet, reviewing observational, experimental, and modelling evidence, we argue that increased coherence of pedestrians' foot placement is a consequence of, not a cause of the instability. Instead, uncorrelated pedestrians produce positive feedback, through negative damping on average, that can initiate significant lateral bridge vibration over a wide range of natural frequencies. We present a simple general formula that quantifies this effect, and illustrate it through simulation of three mathematical models, including one with strong propensity for synchronisation. Despite subtle effects of gait strategies in determining precise instability thresholds, our results show that average negative damping is always the trigger. More broadly, we describe an alternative to Kuramoto theory for emergence of coherent oscillations in nature; collective contributions from incoherent agents need not cancel, but can provide positive feedback on average, leading to global limit-cycle motion.

Explanation of the onset of bouncing cycles in isotropic rotor dynamics; a grazing bifurcation analysis.

The dynamics associated with bouncing-type partial contact cycles are considered for a 2 degree-of-freedom unbalanced rotor in the rigid-stator limit. Specifically, analytical explanation is provided for a previously proposed criterion for the onset upon increasing the rotor speed Ω of single-bounce-per-period periodic motion, namely internal resonance between forward and backward whirling modes. Focusing on the cases of 2 : 1 and 3 : 2 resonances, detailed numerical results for small rotor damping reveal that stable bouncing periodic orbits, which coexist with non-contacting motion, arise just beyond the resonance speed Ω p:q . The theory of discontinuity maps is used to analyse the problem as a codimension-two degenerate grazing bifurcation in the limit of zero rotor damping and Ω = Ω p:q . An analytic unfolding of the map explains all the features of the bouncing orbits locally. In particular, for non-zero damping ζ, stable bouncing motion bifurcates in the direction of increasing Ω speed in a smooth fold bifurcation point that is at rotor speed O ( ζ ) beyond Ω p:q . The results provide the first analytic explanation of partial-contact bouncing orbits and has implications for prediction and avoidance of unwanted machine vibrations in a number of different industrial settings.

Tropical tree cover in a heterogeneous environment: A reaction-diffusion model.

Observed bimodal tree cover distributions at particular environmental conditions and theoretical models indicate that some areas in the tropics can be in either of the alternative stable vegetation states forest or savanna. However, when including spatial interaction in nonspatial differential equation models of a bistable quantity, only the state with the lowest potential energy remains stable. Our recent reaction-diffusion model of Amazonian tree cover confirmed this and was able to reproduce the observed spatial distribution of forest versus savanna satisfactorily when forced by heterogeneous environmental and anthropogenic variables, even though bistability was underestimated. These conclusions were solely based on simulation results for one set of parameters. Here, we perform an analytical and numerical analysis of the model. We derive the Maxwell point (MP) of the homogeneous reaction-diffusion equation without savanna trees as a function of rainfall and human impact and show that the front between forest and nonforest settles at this point as long as savanna tree cover near the front remains sufficiently low. For parameters resulting in higher savanna tree cover near the front, we also find irregular forest-savanna cycles and woodland-savanna bistability, which can both explain the remaining observed bimodality.