Stochastic Delay Differential Equations: A Comprehensive Approach for Understanding Biosystems with Application to Disease Modelling

: Mathematical models have been of great importance in various ﬁelds, especially for understanding the dynamical behaviour of biosystems. Several models, based on classical ordinary differential equations, delay differential equations, and stochastic processes are commonly employed to gain insights into these systems. However, there is potential to extend such models further by combining the features from the classical approaches. This work investigates stochastic delay differential equations (SDDEs)-based models to understand the behaviour of biosystems. Numerical techniques for solving these models that demonstrate a more robust representation of real-life scenarios are presented. Additionally, quantitative roles of delay and noise to gain a deeper understanding of their inﬂuence on the system’s overall behaviour are analysed. Subsequently, numerical simulations that illustrate the model’s robustness are provided and the results suggest that SDDEs provide a more comprehensive representation of many biological systems, effectively accounting for the uncertainties that arise in real-life situations.


Introduction
Over the years, many models have been formulated to closely represent the behaviour of biological systems using different forms of differential equations [1,2].Until recently, many such models have ignored the potential impact of randomness within the system, which may be due to the complexity inherent in solving these models.Similarly, several models have not adequately included the effect of time delays in the system representations, and it is only recently that there has been a growing interest in understanding the implications of delays in dynamic systems.However, there is still relatively less work on models that combine both the influences of randomness and delay within the system, a gap that warrants attention [3].Notably, the limited application of stochastic delay differential equations (SDDEs) is partly attributed to their inherent complexity in finding solutions, which can be particularly challenging when real-time simulation or rapid iterative modelling is required [4][5][6].Moreover, SDDEs may demand significant computational resources and specialised numerical methods to ensure accuracy, rendering them less accessible to researchers with limited computational capabilities.These challenges underscore the imperative need for ongoing research aimed at developing more efficient numerical techniques and strategies to enhance the usability of SDDEs, particularly in large-scale and computationally demanding scenarios.Nevertheless, it is known that many real-life physical and biological systems are often exposed to the influences of several conditions that are not fully understood, making the explicit modelling of such a scenario less feasible.As such, there is an increasing need to extend the existing deterministic models into a more rigorous model that could incorporate the variations in the system dynamics.One of the ways of modelling such dynamics is an inclusion of the stochastic influence often referred to as the noise [7].The natural extension of the deterministic model based on the ordinary differential equations (ODEs) could be an inclusion of delay or noise which leads to the system of stochastic or delay differential equation where the basic difference is that the SDEs account for the influence of randomness while DDEs account for the impact of pauses within the system [8].Similarly, the standard extension of the DDEs-based model would therefore be a stochastic delay differential equation where these equations incorporated the effect of both noise and delay.However, the SDDEs are not just considered as the extension of the DDEs or SDEs, it is widely regarded as the generalisation of both the deterministic delay differential equations and the stochastic ordinary differential equations (SODEs) [9,10].To date, it has been widely accepted that stochastic differential equations play a significant role in the modelling of several physical phenomena such as population dynamics, epidemiology, ecology, and many more.In the same way, the delay differential equation has been greatly used for the analysis and predictions in various areas of life sciences.The time delays included in the models often account for the dependency of the present state on the history [10].Without loss of generality, both noise and delays are increasingly inevitable in biological systems and, due to their importance, an important part of the literature is now devoting attention towards the development and solutions.Some of the early references that discussed the SDDEs are the books by [11] which presented the small delay approximation of the equation and considered the effects of noise on the model equation.The basic introduction to the numerical analysis of these types of equations was presented by [3], where they considered Itô's form of the SDDE and provided a numerical simulation of the solution using the Euler-Maruyama scheme.The survey of results involving the theory and stability of the differential equations involving perturbation and delays are documented in [12].In their work, they discussed the survey of the existence and uniqueness of the solution of SDDEs together with the Markov properties, stochastic stability, and parameter estimations.Many more authors further extend the popular stochastic differential equation to incorporate delays and are used to describe technical devices such as control circuits, where the delay represents the time taken by a signal to travel to the control object together with the reaction time, while the noise accounts for the significant impact of randomness in the circuit performance that could influence the application with high precision.In finance, this equation could be used in developing an option pricing formula which will lead to the delayed Geometric Brownian Motion (GBM) and subsequently help obtain a modified Black-Scholes formula popularly used in economics for pricing European options.In summary, a delay is primarily introduced due to hidden processes that are not fully understood but are known to cause time lag [13,14], while the stochastic part is often included to account for the effect of randomness or behavioural representation which cannot be easily determined, and this could affect the overall dynamics of the system.In the subsequent section, we shall present the mathematical representation of this form of equation with their illustration in various physical modellings.

General Formulation
Since SDDEs emerge from both delay and stochastic differential equations, it is expected that some of the properties of stochastic and delay differential equations would apply to stochastic delay differential equations.Thus, we define an SDDE as a differential equation that includes at least one stochastic process term denoted by W(t) and at least one time delay represented by τ.Mathematically, let 0 ≤ t 0 < t 1 , • • • , t n < T < ∞, and let W(t) be a one-dimensional Wiener process defined on a filtered probability space (Ω, F , P), where Ω represents the sample space, F denotes the filtration, and P is the probability.Then, a stochastic delay differential equation can be represented as follows: W(t) represent the Wiener process with zero mean and variance σ 2 , while the drift coeffi- are d-dimensional and assumed to be Lipschitz continuous.Furthermore, the equation given in (1) can be expressed in an integral form as ( The second integral in ( 2) is often referred to as the stochastic integral in the Itô sense [15].

Numerical Scheme for SDDEs
In this section, we outline the numerical approach for solving Equation (2).To begin, we establish a set of equidistant mesh points within the interval [0, T], denoted as h = T N and t n = nh, where n ranges from 0 to N. Additionally, we assume that, for a given step size h, there exists an integer value, denoted as m, such that the time delay can be expressed as τ = mh.For all instances where the indices satisfy n − m ≤ 0, we set ỹn−m equal to Φ(t n − τ).In all other cases, the numerical approximation of (1) takes the following form: ( We note that the increment function φ(h, ỹn , ỹn−m , I φ ) represents the computation which updates the numerical approximation of the solution at the next time step based on the information available at the current time step n and the delayed information from time step n − m.Similarly, the increment function φ(h, ỹn , ỹn−m , I φ ): R × R −→ R includes a finite number of multiple Itô-integrals [9,16,17] of the form: where j i ∈ {0, 1} and dW 0 (t) = dt , and with t = t n for (3), we denote I φ the collection of Itô-integrals required to compute the increment function φ.

Euler-Maruyama Scheme
The Euler-Maruyama scheme is a popular method for numerically solving Stochastic Differential Equations (SDEs).It is also widely used for solving Stochastic Delay Differential Equations (SDDEs).The basic idea behind the method is to discretise the time domain into small intervals and approximate the solution at the end of each interval by adding the deterministic and stochastic terms.The method is derived from the Taylor expansion of the solution, and it is a first-order weak scheme that converges in the mean-square sense.Suppose we consider the SDDE given in (1), where Y(t) is the state variable at time t, f (t, Y(t), Y(t − τ)) is the drift term, g(t, Y(t), Y(t − τ)) is the diffusion term, W(t) is the Wiener process, and τ is the delay time.The delay term makes the numerical solution more challenging since it involves the state variable at a previous time.To apply the Euler-Maruyama scheme to SDDEs, we first discretise the time domain into N equal intervals of length ∆t = (T − t 0 )/N, where T is the final time and t 0 is the initial time.Then, we can approximate the solution at the end of each interval using the following formula: Y(t n ) is the numerical solution at time t n , and ∆W(t n ) is the Wiener increment over the time interval [t n , t n + 1].The Wiener increments are normally distributed with mean 0 and variance ∆t, which implies that ∆W(t n ) ∼ N(0, ∆t).To obtain the numerical solution, we start from the initial condition Y(t 0 ) = Y 0 and compute Y(t n + 1) recursively using the above formula.We generate the Wiener increments using a random number generator, and we use interpolation or extrapolation techniques to evaluate the delay term ).The Euler-Maruyama scheme is a simple and efficient method for solving SDDEs, but it is only a first-order weak scheme that may not be accurate for complex SDDEs.However, it can be extended to higher-order schemes such as the Milstein method to improve the accuracy.

Milstein Scheme
The Milstein scheme is a popular numerical method for approximating the solutions of stochastic delay differential equations.It is an extension of the Euler-Maruyama method, which is a simple but often inaccurate method for approximating solutions of stochastic differential equations [18].The Milstein scheme improves on the Euler-Maruyama method by including a correction term that takes into account the second-order moments of the Wiener process.
Consider the one-dimensional stochastic delay differential equation given in (1).The Milstein approximation for the solution Y(t) at a time step t n with time increment ∆t is given by where ∆W n = W t n+1 − W t n is the increment of the Wiener process over the time step ∆t.
The first term on the right-hand side of (6) represents the deterministic drift term, which accounts for the delayed influence of the process.The second term is the random term, which accounts for the stochastic fluctuations in the process.The third term is the correction term, which accounts for the second-order moments of the Wiener process.The fourth term is a further correction term that accounts for the covariance between the Wiener process at different time steps.We present the algorithm used for Milstein scheme in Listing 1.The above pseudo-code implements the numerical scheme for a one-dimensional stochastic delay differential equation with the delayed term and the noise term.The solution array Y is initialised with the initial value Y 0 , and then the scheme is applied iteratively to compute the solution values at each time step.For the Milstein scheme, it includes four terms: the deterministic drift term, the random term, the correction term, and the covariance correction term.Each term is computed using the current solution value and the Wiener process increment at the current time step.The scheme has a strong convergence rate of order 1.5 for linear equations and order 1 for non-linear equations [19].

The Evolution of Modelling to SDDEs
Here, we outline the evolution of mathematical modelling, progressing from ordinary differential equations to the system of stochastic delay differential equations.We emphasise the expansion of each model to account for the complexity of overall system behaviour.To begin, we investigate population dynamics.Subsequently, we delve into the applications of SDDEs in various modelling domains, including epidemiology and ecology.Finally, we present the development of the model equations and discuss the rationale behind modifying existing models.

Population Dynamics
It is well-established that population dynamics discusses the changes in the size and composition of populations over time.The study of population dynamics has a long and fascinating history, from early mathematical models of exponential growth to more recent models that incorporate delays and stochasticity.Here, we will explore the evolution of population dynamics models from exponential growth to stochastic delay differential equations.The simplest model of population growth is the exponential growth model.This model assumes that the rate of population growth is proportional to the size of the population.Mathematically, this can be represented as: where P is the population size, t is time, and r is the intrinsic growth rate.This model predicts that populations will grow without bounds unless there are limiting factors, such as resource availability [20].As a result of the limitation of this equation, Verhulst proposed the logistic equation, which is a modification of the exponential growth model that takes into account limiting factors.The logistic equation assumes that as the population size approaches a carrying capacity, the rate of growth slows down.Mathematically, this can be represented as The parameter K represents the carrying capacity and the logistic equation predicts that populations will stabilise at a carrying capacity.The equation has been widely used to model growth in different phenomena and was accepted to better explain population dynamics.Even though this model is widely used, it is worth noting that the logistic equation assumes that the effects of limiting factors on population growth are instantaneous.However, in reality, there are often delays between changes in limiting factors and changes in population growth rates.This led to the further modification by Hutchinson [21] to have a delay logistic equation that takes into account these delays.Mathematically, this is represented as The difference is the inclusion of the delay parameter τ, which thereby allows for the equation to predict that populations will exhibit oscillatory behaviour before stabilising at a carrying capacity.We note that the delay logistic equation can be transformed to become dy dt = ry(1 − y(t − τ)).By considering the equation, we could observe that the equation assumes that the effects of limiting factors on population growth are deterministic.In reality, there are often random fluctuations in the environment that affect population growth rates.Hence, the impact of noise could also be introduced into the system such that we will have r(t) = r(t) + η(t), where η(t) is the Gaussian white noise with a time-varying intensity σ 2 (t) [22].Thus, we could write η(t)dt = σ(t)dW t , with W t representing the Wiener process.Therefore, the stochastic version of the delay logistic equation could be written as which can therefore predict that populations will exhibit both deterministic and stochastic behaviour.In Figure 1, we present the effect of delay and noise in population dynamics using logistic equation.

Epidemiology
In infectious disease modelling, the population is often classified into different compartments based on the population status.The epidemic SIR (Susceptible-Infectious-Recovered) model has been widely used to study the spread of infectious diseases in populations.The classical model assumes that the population is well-mixed and that the transmission rate is constant over time.However, in real-world scenarios, there are often time lags between the exposure of individuals to the disease and their infectiousness.The classical SIR model is based on ordinary differential equations (ODEs) and assumes that the transition between compartments occurs instantaneously.However, this assumption is not always realistic, especially for diseases with a long incubation period or with a significant delay in the onset of infectiousness.As a result, the SIR model has been extended to delay differential equations (DDEs) which are used to model such delays in the transmission process [9].Similarly, SDEs are used to capture the randomness in the transmission process but much work does not consider the scenario wherein there are both delay and randomness in the transmission process.Hence, we shall present the stochastic delayed model for the widely used epidemic SIR model.Mathematically, the classical SIR model is a set of three coupled ODEs that describe the evolution of the three compartments of a population: susceptible (S), infectious (I), and recovered (R).The equations are given by where β is the transmission rate and γ is the recovery rate.The total population size N = S + I + R is assumed to be constant.The delayed version of this model indicates that there may be a time delay between infection and infectiousness.The simple form of the delay inclusion could be given as where τ is the time delay between infection and infectiousness.The delay SIR model has been used to model the spread of diseases such as tuberculosis [23], COVID-19 [10,[24][25][26] and HIV [27,28], which have a long incubation period.We note that the delay SIR model assumes that the time delay is constant and deterministic.However, in many cases, the time delay may be random and unpredictable.Therefore, to incorporate the randomness in the delay, there is a need for the representation of the system using the stochastic delay SIR model.Hence, the model could take the following form: where W t represent independent Wiener processes and σ is the strength of the noise, which are motivated by the randomness and uncertainties inherent in the spread of the disease.We present the illustration of the effect of delay and noise in population dynamics using SIR model in Figure 2.

Ross-Macdonald Malaria Model
The Ross-Macdonald model is a well-known mathematical model used to study the transmission dynamics of malaria.The model started with Ross in 1911 and was designed to describe the transmission of malaria.The concepts were a source of inspiration for a group of medical entomologists, and beginning in 1950, there was a significant advancement in the field.This was largely due to the combination of the theoretical work of George Macdonald and the empirical work of his colleagues [29].Their work provided a complete understanding of transmission dynamics and control, and successfully linked the models and the metrics.To date, this model serves as the basis of countless models for vectorborne diseases.However, we consider the Aron and May version of the Ross-Macdonald model [30].This famous model used the ordinary differential equation, which is given as where x(t) and y(t) are the proportion of infected humans and mosquitoes at time t is the prevalence of malaria, m is the number of mosquitoes per human host, the pathogen's host is denoted by a, the proportion of bites by infectious mosquitoes that infect a human is represented by b.Altogether, the product mabz is the force of infection or "happenings" rate for human infections and the instantaneous death rate is denoted by ν.Further explanation of each parameter can be found in [29,31].By considering this model, it could be observed that the influence of the time-varying nature of the vector and host populations was not considered.To account for this, there is therefore a need to incorporate a time delay.
According to [31], the feedback dynamics from mosquito to human and back to mosquito involve considerable time delays due to the incubation periods of the parasites.Therefore, a delayed model is required to explicitly account for the incubation periods of parasites within the human and the mosquito.Thus, this leads to an extension of the model to a delayed Ross-Macdonald model which is given as [31] where τ 1 , τ 2 represents the delays.These delays represent the time it takes for an infected mosquito to transmit the disease to a human, and the time it takes for an infected human to become infectious.Ruan et al. [31] used this model to calculate the basic reproduction number R 0 and carry out the sensitivity analysis of the R 0 on the incubation periods.As the model adequately accounts for the effect of time delays on the basic reproduction number, it does not include the influence of additional disturbance and randomness.Hence, there is need for a further modification that would include the impact of delays and fluctuations.The delay would capture the fact that infected individuals do not become infectious immediately after being bitten by a mosquito, but instead take a certain amount of time to develop symptoms and become infectious themselves.Similarly, recovered individuals do not immediately lose their immunity to the disease, but instead, retain it for a certain amount of time, while the addition of a stochastic part will account for the random fluctuations in the transmission and recovery rates.This will therefore lead to a stochastic delay differential equation (SDDE) model of the form where W t are independent Wiener processes representing the random fluctuations within the population and σ i denoting the corresponding deviations.The illustration of the effect of delay and noise using the malaria model ( 15) is presented in Figure 3.

Numerical Illustration
In this section, we present the numerical simulation of the SIR model with vaccination and quarantine which is presented in [32].The model considered a population of five compartments which are grouped into Susceptible S(t), Vaccinated class V(t), Infected class I(t), Quarantine class Q(t), and recovered class R(t).Total population is expressed as N(t) = S(t) + V(t) + I(t) + Q(t) + R(t).Assuming the population is constant, the birth rate is denoted by µ, while β represents the transmission rate (the rate at which susceptible individuals become infected).γ is the transmission rate and the rate at which susceptible individuals become vaccinated is represented by ω 1 .Similarly, the quarantine rate is denoted by θ 1 , while the rate at which quarantined individuals recovered is represented by α 2 , while α 1 is the rate at which infected individuals recover and become immune.Hence, the model is presented as The description of parameters in the model ( 16) is presented in Table 1.Firstly, we shall present the quantitative analysis of this model after which we will incorporate the delay parameter and randomness in the model.We will discuss the positivity condition and then examine the equilibria as well as the reproduction number.

Positivity of the Solution
Here, we show that the model ( 16) satisfies the positivity condition provided that N(t) = S(t) + V(t) + I(t) + Q(t) + R(t).Then, by considering each compartment of ( 16 Hence, N(t) at any time t is given by, where lim

Equilibrium
We shall determine both the disease-free equilibrium and the endemic equilibrium.To find the disease-free equilibrium, we take I = Q = R = 0, which means that there are no infected, quarantined, or recovered individuals in the population.By substituting into the equations, we have From the second equation, we obtain V = ω 1 γ S. Therefore, Now, substituting this expression for S into the equation for V, we obtain Thus, the disease-free equilibrium of the system is , 0, 0, 0 .
Similarly, for the endemic equilibrium, we need to solve the system of equations obtained by setting the derivatives to zero: By Substituting each variable into the equation or alternatively, the equilibrium point could be obtained using the MATLAB command "vpasolve".Then, we obtain Substituting this expression for S into the equation for V, we obtain Using the equation for I and the expression for Q obtained by setting the derivative to zero, we have ; where In the same way, we have . Finally, where Therefore, the equilibrium points of the system are given by

Reproduction Number (R 0 )
To calculate the basic reproduction number for this model, we need to determine the expected number of secondary infections caused by one infected individual in a completely susceptible population.The derivation of R 0 involves finding the spectral radius of the next-generation matrix (NGM), which is a matrix that describes the expected number of secondary infections caused by each compartment.The NGM for this model is Next, we find the spectral radius of this matrix where we need to find the eigenvalues and choose the one with the largest absolute value.Thus, we have the eigenvalues as The spectral radius of the NGM, which represents the basic reproduction number R 0 for the SVIQR model, is given by

Numerical Illustration
We display the model behaviour for each compartment, and to further understand the contributions of the parameters on R 0 , we shall perform the sensitivity analysis.
Figure 4 shows the behaviour of individuals within the population.It shows that most of the population is in the susceptible compartment, meaning they have not been exposed to the disease and are susceptible to infection.By considering the susceptible curve, it indicates that, as time progresses, the number of infected individuals increases as the disease spreads from person to person.This increase in the number of infected individuals leads to a decrease in the number of susceptible individuals as some people become infected and move from the susceptible compartment to the infected compartment.Eventually, the number of infected individuals reaches a peak and begins to decline as some individuals recover from the disease or die.As the number of infected individuals decreases, the number of removed individuals (those who have recovered from the disease or died) increases.In addition, the result also shows that, as the number of vaccinated individuals increases over time, there is a slow spread of the disease, which therefore reduces the number of susceptible individuals.From a biological perspective, this result shows how vaccination can be an effective strategy for reducing the spread of infectious diseases.By increasing the number of vaccinated individuals, we can reduce the number of susceptible individuals and slow the spread of the disease.Furthermore, it also suggests how the number of infected individuals can peak and decline over time as individuals recover from the disease or die.
Next, we plot the reproduction number against the dependent parameters to understand how the R 0 changes with a change in parameter values.The output displayed in Figure 5 shows how R 0 varies with different values of [γ, θ 1 , α 1 ].As the value of these parameters increases, the transmission potential of the disease decreases, which leads to a lower value of R 0 .The result exhibits a curve that starts high at low values of the parameters, and then gradually decreases as the parameter increases.This behaviour of the curve can be interpreted as a trade-off between the transmission potential of the disease and the ability of the population to mount an immune response.At low values of [γ, θ 1 , α 1 ], the transmission potential is high, but the population has a weak immune response, so R 0 is high.Similarly, at high values of [γ, θ 1 , α 1 ], the transmission potential is low, which implies that the population has a strong immune response, and therefore R 0 will be low.Hence, the minimum value of R 0 would correspond to the optimal balance between transmission potential and immune response.Figure 5 shows the relationship between the reproduction number (R 0 ) and the transmission rate (β) for different values of the recovery rate (γ), quarantine rate (θ 1 ), and quarantine recovery rate (α 1 ).The result shows that R 0 decreases with γ, θ 1 , and α 1 .In addition, it shows that the effect of α 1 on R 0 is more pronounced at lower values of γ and θ 1 .This implies that when γ and θ 1 are low, infected individuals are more likely to remain infectious and transmit the disease.

Model with Incorporated Delay
Here, we present the modification of this model to incorporate the time delay with the assumption that it takes a time period for the population to move from one compartment to the other.Hence, a modified version of (16), is given as In the model ( 20), we have added a time delay τ to the terms involving the past values of S, I, and Q.The time delay is introduced to account for the fact that there is a lag between the time a person is infected with the disease and the time they become infectious.The delay also allows for the time it takes for an infected person to seek treatment and be moved to a quarantine centre.The modified model for the susceptible population includes a delay term βS(t − τ)I(t − τ), which accounts for the time delay between the time an individual becomes infected and the time they become infectious.The recovered populations also include delay term α 1 I(t − τ), which represents the delay between the time an individual becomes infected and the time they recover.θ 1 I(t − τ) accounts for the delay between the time an individual is diagnosed with the disease and the time they are quarantined.The vaccinated population does not include a delay term, since vaccination is assumed to take effect immediately.Additionally, it is important to note that the choice of time delays depends on the specific disease and its transmission dynamics, and may require careful analysis or calibration to obtain accurate results.Additionally, incorporating time delays can make the model more complex and may require additional assumptions or parameters.
By comparing the output of Figure 4 with the result obtained in Figure 6, it could be observed that the delay has a significant impact on the behaviour of the system.In the model which does not incorporate any delay, the epidemic peaks rapidly and then gradually declines as the population becomes immune.The overall duration of the epidemic is relatively short, and the peak of the epidemic occurs at a relatively early time point.In contrast, the model which incorporates a delay in the transmission of the disease shows that the epidemic peak is delayed, and the epidemic lasts for a longer time.This delay is due to the time required for individuals to become infectious after being exposed to the disease.As a result, the epidemic spreads more slowly and gradually, and the peak of the epidemic occurs at a later time point.Therefore, we can conclude that incorporating delay in the transmission of the disease has a significant impact on the behaviour of the epidemic and can lead to different outcomes compared to a system without delay.

Model with Incorporated Delay and Additive Noise
The idea of incorporating noise into the system is born from the understanding that many real-world systems exhibit complex behaviour and are inherently noisy, with many exhibiting stochastic behaviour due to the presence of random fluctuations.In infectious disease models, the spread of disease is influenced by a multitude of factors that are difficult to model precisely.To account for some of the uncertainty in the system, noise can be added to the model.Moreover, many of the parameters in infectious disease models are uncertain or difficult to estimate precisely, adding further uncertainty to the model.By adding noise to the model, this uncertainty can be accounted for, resulting in a more realistic representation of the system.Additionally, non-linear dynamics can also contribute to the complexity of the system, even in the absence of external noise.The interactions between the various components of the system can give rise to hidden patterns and structures that may be difficult to discern otherwise.Hence, adding noise to the model can reveal these patterns and structures, allowing for a deeper understanding of the system.However, there are some challenges in adding noise to a model as there could be difficulty in determining the appropriate level of noise to add without overfitting or underfitting the data.Additionally, the choice of noise distribution can also have an impact on the model results.Therefore, when adding noise to infectious disease models, it is important to consider the characteristics of the disease being modelled and the specific research question being addressed.The type of disease, the available data, and the research objectives all play a role in determining the appropriate level and type of noise to incorporate into the model.Here, we are not looking into a specific dataset and we will only add a sample representation of noise into our system to show how the system behaviour changes when noise is included.Hence, we will represent our model as σ S,V,I,Q,R represent the noise associated to each model category.The assumption adopted here is that there is a level of uncertainty in the spread of the diseases and each compartment does exhibit random behaviour.A biological justification for noise inclusion in vaccination could be attributed to the existing knowledge that the efficacy of vaccines could vary based on factors such as the age and health status of the individual, the timing and number of vaccine doses, and the prevalence of different strains of the disease, etc.Thus, including noise in the vaccinated compartment V(t) can help account for this uncertainty and provide a more realistic representation of the effects of vaccination on the spread of the disease.A similar argument can be generated for other compartments depending on the type of disease under consideration.Figure 7 displays the behaviour of the model with the noise.We could note that without the noise term, the model would be deterministic and assume that all individuals follow the same path of infection and recovery.However, in reality, there are always variations and uncertainties in how individuals are infected, vaccinated, or recover.The result here, therefore, captures the level of randomness exhibited in the system.Also, the result of the inclusion of noise in the model allows for the assessment of the variability in the predicted outcomes of the model, which can provide valuable insights into the robustness of the system and help identify areas of improvement for disease control measures.

Conclusions
Mathematical models have played a crucial role in the understanding of various biosystems.However, classical approaches based on ordinary differential equations (ODEs), delay differential equations (DDEs), or stochastic differential equations (SDEs) may not fully capture the real-life scenarios of these systems.Thus, the need arises for a more comprehensive approach that combines the features of these classical approaches, and this leads to the development of Stochastic Delay Differential Equations (SDDEs).In this work, we have investigated the application of SDDE-based models in understanding the behaviour of biosystems, particularly in epidemiology.We have presented a numerical technique for solving the model and shown that the SDDEs provide a more robust representation of real-life scenarios.The analysis of the quantitative and qualitative role of delay and noise reveals their significant influence on the overall behaviour of the system.The inclusion of delay terms in the model helps to capture the dynamics of the system, particularly in cases where there is a time lag between the infection and the manifestation of symptoms or the recovery process.
Similarly, the incorporation of stochasticity accounts for the uncertainty that arises in real-life scenarios and provides a more realistic representation of the system.We have presented numerical simulations to illustrate the robustness of the model, and the results suggest that the SDDEs provide a richer representation of many biological systems.Furthermore, the effects of delay and noise on the system have been investigated, and our results show that both delay and noise significantly impact the system's behaviour.Delay leads to oscillatory behaviour in the system, and the amplitude of the oscillations depends on the size of the delay.On the other hand, the stochasticity introduced by noise leads to fluctuations in the system, and the size of the fluctuations depends on the strength of the noise.The combined effect of delay and noise leads to complex behaviour in the system, with a greater degree of variability and unpredictability than in deterministic systems.
Hence, the model devised in this study stands poised to significantly enhance our understanding of biosystem dynamics, with a particular focus on epidemiology.The utilisation of SDDE-based models offers a substantial leap towards crafting a more resilient representation of real-world scenarios by adeptly encompassing the dual factors of delay and stochasticity.Our numerical simulations robustly affirm that the model not only furnishes precise and efficient solutions but also holds substantial potential in predicting the intricate behaviours exhibited by biological systems.Furthermore, the insights garnered from this model can be instrumental in advancing the development of intervention strategies for combating diseases, thereby fortifying our ability to safeguard public health.
For future research directions, a pivotal area of exploration is the development of robust parameter estimation techniques tailored specifically for SDDEs.The accurate determination of model parameters from empirical data is critical for the practical application of SDDEs in epidemiology and other fields.Therefore, innovative approaches and methodologies for estimating these parameters will undoubtedly contribute to the model's refinement and its ability to address real-world challenges effectively.In addition, we acknowledge that, in more practical scenarios, individual variables may indeed have different time delays.Future extensions of our work would be to explore the variations in delay values to provide a more nuanced understanding of the system's behaviours.

Listing 1 . 1 .
Implementation algorithm.Define the initial condition Y[0] and the time interval [0,T].2. Define the time step and the number of time steps N = (T−t_0)/dt 3. Initialize the solution array Y with Y[0] and the Wiener process array W with increments W(t_n).4. For n = 1 to N do the following: a Compute the drift term f(t_n,Y(t_n),Y(t_n−tau)) and the diffusion term g(t_n,Y(t_n),Y(t_n−tau)).b Generate a Wiener increment W(t_n) using a random no generator.c Compute the approximation Y(t_{n+1}) using the Euler-Maruyama or Milstein formula.d Append Y(t_{n+1}) to the solution array Y and W(t_n) to the Wiener process array W. 5.Return the solution array Y and the Wiener process array W.

Figure 2 .
Illustration of the effect of delay and noise in population dynamics.(a) SIR model without noise and delay with β = 0.5, γ = 0.1.(b) SIR model incorporated with noise and delay with

Table 1 .
Description of the model parameters.
I(t)Infectious population at time t.The group of individuals currently infected and capable of spreading the disease.Q(t)Quarantined population at time t.Individuals who are quarantined due to infection.R(t) Recovered population at time t.Includes individuals who have recovered from the disease and are no longer infectious.µ Recruitment rate of the susceptible population, i.e., the rate at which individuals enter the susceptible group.β Transmission rate of the disease which is the rate at which susceptible individuals become infected when they come into contact with infectious individuals.τ Time delay indicating that newly infected individuals do not become infectious immediately but after a delay of τ time units.γ Recovery rate, i.e., the rate at which infected individuals recover and move to the recovered compartment.σ S , σ V , σ I , σ Q , σ R Represent the level of fluctuation within each model compartment.