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.

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.

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.

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.

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.

Author Correction: Amazonian forest-savanna bistability and human impact.

This corrects the article DOI: 10.1038/ncomms15519.

Amazonian forest-savanna bistability and human impact.

A bimodal distribution of tropical tree cover at intermediate precipitation levels has been presented as evidence of fire-induced bistability. Here we subdivide satellite vegetation data into those from human-unaffected areas and those from regions close to human-cultivated zones. Bimodality is found to be almost absent in the unaffected regions, whereas it is significantly enhanced close to cultivated zones. Assuming higher logging rates closer to cultivated zones and spatial diffusion of fire, our spatiotemporal mathematical model reproduces these patterns. Given a gradient of climatic and edaphic factors, rather than bistability there is a predictable spatial boundary, a Maxwell point, that separates regions where forest and savanna states are naturally selected. While bimodality can hence be explained by anthropogenic edge effects and natural spatial heterogeneity, a narrow range of bimodality remaining in the human-unaffected data indicates that there is still bistability, although on smaller scales than claimed previously.

The impacting cantilever: modal non-convergence and the importance of stiffness matching.

The problem of an Euler-Bernoulli cantilever beam whose free end impacts with a point constraint is revisited from the point of view of modal analysis. It is shown that there is non-uniqueness of consistent impact laws for a given modal truncation. Moreover, taking an N-mode compliant, bilinear formulation and passing to the rigid limit leads to a sequence of impact models that does not converge as N--> ∞. The dynamics of such truncated models are studied numerically and found to give rise to quite different dynamics depending on the number of degrees of freedom taken. The simulations are compared with results from simple experiments that show a propensity for multiple-tap dynamics, in which higher-order modes lead to rapidly cycling intermittent contact. The conclusion reached is that, to derive an accurate model, one needs to avoid the impact limit altogether, and take sufficiently many modes in the formulation to match the actual stiffness of the constraining stop.