Past Junior Applied Mathematics Seminar

We present models of a vapour bubble produced during ureteroscopy and laser lithotripsy treatment of kidney stones. This common treatment for kidney stones involves passing a flexible ureteroscope containing a laser fibre via the ureter and bladder into the kidney, where the fibre is placed in contact with the stone. Laser pulses are fired to fragment the stone into pieces small enough to pass through an outflow channel. Laser energy is also transferred to the surrounding fluid, resulting in vapourisation and the production of a cavitation bubble.

While in some cases, bubbles have undesirable effects – for example, causing retropulsion of the kidney stone – it is possible to exploit bubbles to make stone fragmentation more efficient. One laser manufacturer employs a method of firing laser pulses in quick succession; the latter pulses pass through the bubble created by the first pulse, which, due to the low absorption rate of vapour in comparison to liquid, increases the laser energy reaching the stone.

As is common in bubble dynamics, we couple the Rayleigh-Plesset equation to an energy conservation equation at the vapour-liquid boundary, and an advection-diffusion equation for the surrounding liquid temperature.1 However, this present work is novel in considering the laser, not only as the cause of nucleation, but as a spatiotemporal source of heat energy during the expansion and collapse of a vapour bubble.
 

Numerical and analytical methods are employed alongside experimental work to understand the effect of laser power, pulse duration and pulse pattern. Mathematically predicting the size, shape and duration of a bubble reduces the necessary experimental work and widens the possible parameter space to inform the design and usage of lasers clinically.

Optimal execution of large positions over a given trading period is a fundamental decision-making problem for financial services. In this talk we explore reinforcement learning methods, in particular policy gradient methods, for finding the optimal policy in the optimal liquidation problem. We show results for the case where we assume a linear quadratic regulator (LQR) model for the underlying dynamics and where we apply the method to the data directly. The empirical evidence suggests that the policy gradient method can learn the global optimal solution for a larger class of stochastic systems containing the LQR framework, and that it is more robust with respect to model misspecification when compared to a model-based approach.

Michael Negus

Modelling high-speed droplet impact onto an elastic membrane

The impact of a high-speed droplet onto an elastic membrane is a highly nonlinear process and poses a formidable modelling challenge due to both the multi-scale nature of the flow and the fluid-structure interaction between the droplet and the membrane. We present two modelling approaches for droplet impact onto elastic membranes: matched asymptotics and direct numerical simulations (DNS). Inviscid Wagner theory is used in the former to derive analytical expressions which approximate the behaviour of the droplet during the early stages of impact, while the DNS builds on the open-source volume-of-fluid code Basilisk. We demonstrate the strong influence that the thickness, tension and stiffness of the membrane have on the dynamics of the droplet and the membrane. We also quantitatively show that the speed the droplet spreads across the substrate is notably decreased when the membrane is more compliant, which is consistent with experimental findings that splashing can be inhibited by impacting onto a soft substrate. We conclude by showing how these methods are complementary, as a combination of both can lead to a thorough understanding of the droplet impact across timescales.

Helen Saville

Lubrication model of a valve-controlled, gravity-driven bioreactor

Hospitals sometimes experience shortages of donor blood platelet supplies, motivating research into in vitro production of platelets. We model a novel platelet bioreactor described in Shepherd et al. [1]. The bioreactor consists of an upper channel, a lower channel, and a cell-seeded porous collagen scaffold situated between the two. Flow is driven by gravity, and controlled by valves on the four inlets and outlets. The bioreactor is long relative to its width, a feature which we exploit to derive a lubrication reduction of Navier-Stokes flow coupled to Darcy. Models for two cases are considered: small amplitude valve oscillations, and order one amplitude valve oscillations. The former model is a systematic reduction; the latter incorporates a phenomenological approximation for the cross-sectional flow profile. As the shear stress experienced by cells influences platelet production, we use our model to quantify the effect of valve dynamics on shear stress.

1: Shepherd, J.H., Howard, D., Waller, A.K., Foster, H.R., Mueller, A., Moreau, T., Evans, A.L., Arumugam, M., Chalon, G.B., Vriend, E. and Davidenko, N., Biomaterials, 182, pp.135-144. (2018)

Several problems of great importance in the study of glaciers and ice sheets involve processes of attachment and reattachment of the ice from the bedrock. Consider, for example, an ice sheet sliding from the continent into the ocean, where it goes afloat. Another example is that of subglacial cavitation, a fundamental mechanism in glacial sliding where the ice detaches from the bedrock along the downstream area of an obstacle. Such problems are generally modelled as a viscous Stokes flow with a free boundary and contact boundary conditions. In this talk, I will present a framework for solving such problems numerically. I will start by introducing the mathematical formulation of these viscous contact problems and the challenges that arise when trying to approximate them numerically. I will then show how, given a stable scheme for the free boundary equation, one can build a penalty formulation for the viscous contact problem in such a way that the resulting algorithm remains stable and robust.

Organoids are three–dimensional multicellular tissue constructs. When cultured in vitro, they recapitulate the structure, heterogeneity, and function of their in vivo counterparts. As awareness of the multiple uses of organoids has grown, e.g. in drug discovery and personalised medicine, demand has increased for low–cost and efficient methods of producing them in a reproducible manner and at scale. We are working in collaboration with the biotechnology company Cellesce, who develop bioprocessing systems for the expansion of organoids at scale. Part of their technology includes a bioreactor, which utilises flow of culture media to enhance nutrient delivery to the organoids and facilitate the removal of waste metabolites. A key priority is ensuring uniformity in organoid size and reproducibility; qualities that depends on the bioreactor design and operating conditions. A complete understanding of the system requires knowledge of the spatial and temporal information regarding flow and the resulting oxygen and metabolite concentrations throughout the bioreactor. However, it is impractical to obtain this data empirically, due to the highly–controlled environment of the bioreactor posing difficulties for online real–time monitoring of the system. Thus, we exploit a mathematical modelling approach, to provide spatial as well as temporal information.

In the bioreactor, organoids are seeded as single cells in a layer of hydrogel. We present a general model for the nutrient and waste metabolite concentrations in the hydrogel and organoid regions of the bioreactor. Resolving for the millions of organoids within the hydrogel is computationally expensive and infeasible. Hence, we take a mathematical homogenisation approach to understand how the behaviour of the organoids on the microscale influences the macroscale behaviour in the hydrogel layer. We consider the case of growing organoids, with a temporally and spatially dependent radii, and exploit the separation of scales to systematically derive an effective macroscale model for metabolite transport. We explore some canonical problems to understand our homogenised system.

Many problems in engineering can be understood as controlling the bifurcation structure of a given device. For example, one may wish to delay the onset of instability, or bring forward a bifurcation to enable rapid switching between states. In this talk, we will describe a numerical technique for controlling the bifurcation diagram of a nonlinear partial differential equation by varying the shape of the domain. Our aim is to delay or advance a given branch point to a target parameter value. The algorithm consists of solving a shape optimization problem constrained by an augmented system of equations, called the Moore–Spence system, that characterize the location of the branch points. We will demonstrate the effectiveness of this technique on several numerical experiments on the Allen–Cahn, Navier–Stokes, and hyperelasticity equations.

Tendon tissue engineering aims to grow functional tissue in the lab. Tissue is grown inside a bioreactor which controls both the mechanical and biochemical environment. As tendon cells alter their behaviour in response to mechanical stresses, designing suitable bioreactor loading regimes forms a key component in ensuring healthy tissue growth.  

Linking the forces imposed by the bioreactor to the shear stress experienced by individual cell is achieved by homogenisation using multiscale asymptotics. We will present a continuum model capturing fluid-structure interaction between the nutrient media and the fibrous scaffold where cells grow. Solutions reflecting different experimental conditions will be discussed in view of the implications for shear stress distribution experienced by cells across the bioreactor.  

Neural controlled differential equations (Neural CDEs) are a continuous-time extension of recurrent neural networks (RNNs). They are considered SOTA for modelling functions on irregular time series, outperforming other ODE benchmarks (ODE-RNN, GRU-ODE-Bayes) in offline prediction tasks. However, current implementations are not suitable to be used in online prediction tasks, severely restricting the domains of applicability of this powerful modeling framework. We identify such limitations with previous implementations and show how said limitations may be addressed, most notably to allow for online predictions. We benchmark our online Neural CDE model on three continuous monitoring tasks from the MIMIC-IV ICU database, demonstrating improved performance on two of the three tasks against state-of-the-art (SOTA) non-ODE benchmarks, and improved performance on all tasks against our ODE benchmark.

Joint work with Patrick Kidger, Lingyi Yang, and Terry Lyons.

The index of a saddle point of a smooth function is the number of descending directions of the saddle. While the index can usually be retrieved by counting the number of negative eigenvalues of the Hessian at the critical point, we may not have the luxury of having second derivatives in data deriving from practical applications. To address this problem, we develop a computational pipeline for estimating the index of a non-degenerate saddle point without explicitly computing the Hessian. In our framework, we only require a sufficiently dense sample of level sets of the function near the saddle point. Using techniques in Morse theory and Topological Data Analysis, we show how the shape of saddle points can help us infer the index of the saddle. Furthermore, we derive an explicit upper bound on the density of point samples necessary for inferring the index depending on the curvature of level sets. 

Coronary heart disease is characterised by the formation of plaque on artery walls, restricting blood flow. If a plaque deposit ruptures, blood clot formation (thrombosis) rapidly occurs with the potential to fatally occlude the vessel within minutes. Von Willebrand Factor (vWF) is a shear-sensitive protein which has a critical role in blood clot formation in arteries. At the high shear rates typical in arterial constrictions (stenoses), vWF undergoes a conformation change, unfolding and exposing binding sites and facilitating rapid platelet deposition. 

To understand the effect of  stenosis geometry and blood flow conditions on the unfolding of vWF and subsequent platelet binding, we developed a continuum model for the initiation of thrombus formation by vWF in an idealised arterial stenosis. In this talk I will discuss modelling proteins in flow using viscoelastic fluid models, the insight asymptotic reductions can offer into this complex system and some of the challenges of studying fast arterial blood flows. 

Topological data analysis is a growing area of research where topology and geometry meets data analysis. Many data science problems have a geometric flavor, and thus computational tools like persistent homology and Mapper were often found to be useful. Domains of applications include cosmology, material science, diabetes and cancer research. We will discuss some main tools of the field and some prominent applications.

Diffusion processes are widely used to model the evolution of random values over time. In many applications, the diffusion process is constrained to a finite domain. We consider the estimation problem of a diffusion process constrained by a polytope, i.e. intersection of finitely many (hyper-)planes, given a discretely observed time series data. Since the boundary behaviours of a diffusion process are characterised by its drift and diffusion functions, we derive sufficient conditions on the drift and diffusion functions for the nonattainablity of a polytope. We use deep learning to estimate the drift and diffusion, and ensure that their constraints are satisfied throughout the training.

Multivariate point processes are used to model event-type data in a wide range of domains. One interesting application is to model the emission of electric impulses of biological neurons. In this context, the point process model needs to capture the time-dependencies and interactions between neurons, which can be of two kinds: exciting or inhibiting. Estimating these interactions, and in particular the functional connectivity of the neurons are problems that have gained a lot of attention recently. The general nonlinear Hawkes process is a powerful model for events occurring at multiple locations in interaction. Although there is an extensive literature on the analysis of the linear model, the probabilistic and statistical properties of the nonlinear model are still mainly unknown. In this paper, we consider nonlinear Hawkes models and, in a Bayesian nonparametric inference framework, derive concentration rates for the posterior distribution.  We also infer the graph of interactions between the dimensions of the process and prove that the posterior distribution is consistent on the graph adjacency matrix.

What do random spanning trees, graph embeddings, random walks, simplices and graph curvature have in common? As you may have guessed from the title, they are indeed all intimately connected to the effective resistance on graphs! While originally invented as a tool to study electrical circuits, the effective resistance has proven time and again to be a graph characteristic with a variety of interesting and often surprising properties. Starting from a number of equivalent but complementary definitions of the effective resistance, we will take a stroll through some classical theorems (Rayleigh monotonicity, Foster's theorem), a few modern results (Klein's metricity, Fiedler's graph-simplex correspondence) and finally discuss number of recent developments (variance on graphs, discrete curvature and graph embeddings).

Michael Coughlan (with Sam Howison, Ian Hewitt, Andrew Wells)

Arctic sea ice forms a thin but significant layer at the ocean surface, mediating key climate feedbacks. During summer, surface melting produces considerable volumes of water, which collect on the ice surface in ponds. These ponds have long been suggested as a contributing factor to the discrepancy between observed and predicted sea ice extent. When viewed at large scales ponds have a complicated, approximately fractal geometry and vary in area from tens to thousands of square meters. Increases in pond depth and area lead to further increases in heat absorption and overall melting, contributing to the ice-albedo feedback. 

Previous modelling work has focussed either on the physics of individual ponds or on the statistical behaviour of systems of ponds. In this talk I present a physically-based network model for systems of ponds which accounts for both the individual and collective behaviour of ponds. Each pond initially occupies a distinct catchment basin and evolves according to a mass-conserving differential equation representing the melting dynamics for bare and water-covered ice. Ponds can later connect together to form a network with fluxes of water between catchment areas, constrained by the ice topography and pond water levels. 

I use the model to explore how the evolution of pond area and hence melting depends on the governing parameters, and to explore how the connections between ponds develop over the melt season. Comparisons with observations are made to demonstrate the ways in which the model qualitatively replicates properties of pond systems, including fractal dimension of pond areas and two distinct regimes of pond complexity that are observed during their development cycle. 

Different perimeter-area relationships exist for ponds in the two regimes. The model replicates these relationships and exhibits a percolation transition around the transition between these regimes, a facet of pond behaviour suggested by previous studies. The results reinforce the findings of these studies on percolation thresholds in pond systems and further allow us to constrain pond coverage at this threshold - an important quantity in measuring the scale and effects of the ice-albedo feedback.

We propose a subspace Gauss-Newton method for nonlinear least squares problems that builds a sketch of the Jacobian on each iteration. We provide global rates of convergence for regularization and trust-region variants, both in expectation and as a tail bound, for diverse choices of the sketching matrix that are suitable for dense and sparse problems. We also have encouraging computational results on machine learning problems.

In this talk we will discuss a problem that was worked on during MISGSA 2020, a Study Group held in January at The University of Zululand, South Africa.

We look at a communication network with two types of users - Primary users (PU) and Secondary users (SU) - such that we reduce the network to a set of overlapping sub-graphs consisting of SUs indexed by a specific PU. Within any given sub-graph, the PU may be communicating at a certain fixed frequency F. The respective SUs also wish to communicate at the same frequency F, but not at the expense of interfering with the PU signal. Therefore if the PU is active then the SUs will not communicate.

In an attempt to increase information throughput in the network, we instead allow the SUs to communicate at a different frequency G, which may or may not interfere with a different sub-graph PU in the network, leading to a multi-objective optimisation problem.

We will discuss not only the problem formulation and possible approaches for solving it, but also the pitfalls that can be easily fallen into during study groups.

Stress perfusion cardiac magnetic resonance (CMR) imaging has been shown to be highly accurate for the detection of coronary artery disease. However, a major limitation is that the accuracy of the visual assessment of the images is challenging and thus the accuracy of the diagnosis is highly dependent on the training and experience of the reader. Quantitative perfusion CMR, where myocardial blood flow values are inferred directly from the MR images, is an automated and user-independent alternative to the visual assessment.

This talk will focus on addressing the main technical challenges which have hampered the adoption of quantitative myocardial perfusion MRI in clinical practice. The talk will cover the problem of respiratory motion in the images and the use of dimension reduction techniques, such as robust principal component analysis, to mitigate this problem. I will then discuss our deep learning-based image processing pipeline that solves the necessary series of computer vision tasks required for the blood flow modelling and introduce the Bayesian inference framework in which the kinetic parameter values are inferred from the imaging data.

Many manufacturing processes require the use of robots to transport parts around a factory line. Some parts, which are very thin (e.g. car doors)
are prone to elastic deformations as they are moved around by a robot. These should be avoided at all cost. A problem that was recently raised by
F.E.E. (Fleischmann Elektrotech Engineering) at the ESGI 158 study group in Barcelona was to ascertain how to determine the stresses of a piece when
undergoing a prescribed motion by a robot. We present a simple model using Kirschoff-Love theory of flat plates and how this can be adapted. We
outline how the solutions of the model can then be used to determine the stresses. 

Numerous mathematical models have been proposed for modelling cancerous tumour invasion (Gatenby and Gawlinski 1996), angiogenesis (Owen et al 2008), growth kinetics (Wang et al 2009), response to irradiation (Gao et al 2013) and metastasis (Qiam and Akcay 2018). In this study, we attempt to model the qualitative behavior of growth, invasion, angiogenesis and fragmentation of tumours at the tissue level in an explicitly spatial and continuous manner in two dimensions. We simulate the effectiveness of radiation therapy on a growing tumour in comparison with immunotherapy and propose a novel framework based on vector fields for modelling the impact of interstitial flow on tumour morphology. The results of this model demonstrate the effectiveness of employing a system of partial differential equations along with vector fields for simulating tumour fragmentation and that immunotherapy, when applicable, is substantially more effective than radiation therapy.

Topology optimisation finds the optimal material distribution of a fluid or solid in a domain, subject to PDE and volume constraints. There are many formulations and we opt for the density approach which results in a PDE, volume and inequality constrained, non-convex, infinite-dimensional optimisation problem without a priori knowledge of a good initial guess. Such problems can exhibit many local minima or even no minima. In practice, heuristics are used to obtain the global minimum, but these can fail even in the simplest of cases. In this talk, we will present an algorithm that solves such problems and systematically discovers as many of these local minima as possible along the way.

Recent advances in experimental imaging techniques have allowed us to observe the fine details of how droplets behave upon impact onto a substrate. However, these are highly non-linear, multiscale phenomena and are thus a formidable challenge to model. In addition, when the substrate is deformable, such as an elastic sheet, the fluid-structure interaction introduces an extra layer of complexity.

We present two modeling approaches for droplet impact onto deformable substrates: matched asymptotics and direct numerical simulations. In the former, we use Wagner's theory of impact to derive analytical expressions which approximate the behavior during the early time of impact. In the latter, we use the open source volume-of-fluid code Basilisk to conduct simulations designed to give insight into the later times of impact.

We conclude by showing how these methods are complementary, and a combination of both can give a thorough understanding of the droplet impact across timescales. 

We consider the problem of global minimization with bound constraints. The problem is known to be intractable for large dimensions due to the exponential increase in the computational time for a linear increase in the dimension (also known as the “curse of dimensionality”). In this talk, we demonstrate that such challenges can be overcome for functions with low effective dimensionality — functions which are constant along certain linear subspaces. Such functions can often be found in applications, for example, in hyper-parameter optimization for neural networks, heuristic algorithms for combinatorial optimization problems and complex engineering simulations.

Extending the idea of random subspace embeddings in Wang et al. (2013), we introduce a new framework (called REGO) compatible with any global min- imization algorithm. Within REGO, a new low-dimensional problem is for- mulated with bound constraints in the reduced space. We provide probabilistic bounds for the success of REGO; these results indicate that the success is depen- dent upon the dimension of the embedded subspace and the intrinsic dimension of the function, but independent of the ambient dimension. Numerical results show that high success rates can be achieved with only one embedding and that rates are independent of the ambient dimension of the problem.

Introducing cheap function proxies for quickly producing approximate random numbers, we show convergence of modified numerical schemes, and coupling between approximation and discretisation errors. We bound the cumulative roundoff error introduced by floating-point calculations, valid for 16-bit half-precision (FP16). We combine approximate distributions and reduced-precisions into a nested simulation framework (via multilevel Monte Carlo), demonstrating performance improvements achieved without losing accuracy. These simulations predominantly perform most of their calculations in very low precisions. We will highlight the motivations and design choices appropriate for SVE and FP16 capable hardware, and present numerical results on Arm, Intel, and NVIDIA based hardware.

In a robust decision, we are pessimistic toward our decision making when the probability measure is unknown. In particular, we optimise our decision under the worst case scenario (e.g. via value at risk or expected shortfall).  On the other hand, most theories in reinforcement learning (e.g. UCB or epsilon-greedy algorithm) tell us to be more optimistic in order to encourage learning. These two approaches produce an apparent contradict in decision making. This raises a natural question. How should we make decisions, given they will affect our short-term outcomes, and information available in the future?

In this talk, I will discuss this phenomenon through the classical multi-armed bandit problem which is known to be solved via Gittins' index theory under the setting of risk (i.e. when the probability measure is fixed). By extending this result to an uncertainty setting, we can show that it is possible to take into account both uncertainty and learning for a future benefit at the same time. This can be done by extending a consistent nonlinear expectation  (i.e. nonlinear expectation with tower property) through multiple filtrations.

At the end of the talk, I will present numerical results which illustrate how we can control our level of exploration and exploitation in our decision based on some parameters.