Dynamic analysis of delayed vaccination process along with impact of retrial queues

Sudipa Chauhan; Shweta Upadhyaya; Payal Rana; Geetika Malik

doi:10.1515/cmb-2022-0147

Article Open Access

Dynamic analysis of delayed vaccination process along with impact of retrial queues

Sudipa Chauhan , Shweta Upadhyaya , Payal Rana and Geetika Malik

Published/Copyright: May 22, 2023

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information

From the journal Computational and Mathematical Biophysics Volume 11 Issue 1

Abstract

An unprecedented and precise time-scheduled rollout for the vaccine is needed for an effective vaccination process. This study is based on the development of a novel mathematical model considering a delay in vaccination due to the inability to book a slot in one go for a system. Two models are proposed which involve a delay differential equation mathematical model whose dynamical analysis is done to show how the delay in vaccination can destabilize the system. Further, this delay led to the formulation of a queuing model that accounts for the need to retry for the vaccination at a certain rate as delay in vaccination can have negative repercussions. The transition rates from one stage to another follow an exponential distribution. The transient state probabilities of the model are acquired by applying the Runge-Kutta method and hence performance indices are also obtained. These performance measures include the expected number of people in various states. Finally, numerical analysis is also provided to validate both models. Our results would specifically focus on what happens if the delay time increases or if the retrial rate increases (delay time decreases). The results reveal that a delay in being vaccinated by the first dose (i.e., 80 days) leads to an unstable system whereas there exists a delay simultaneously in getting vaccinated by both doses that destabilize the system early (i.e., 80 and 120 days for dose one and two, respectively). The system destabilizes faster in the presence of a delay for slot booking for both doses as compared to one dose delay. Further, the numerical results of queuing models show that if the retrial rate increases in this delay time to book the slots, it not only increases in the vaccinated class but also increases the recovered population.

Keywords: delay; dynamic analysis; retrial queues; RK method of fourth order; vaccination process

MSC 2010: 35A01; 65L10; 65L12; 65L20; 65L70

1 Introduction

Through analysis of the model, we make an attempt to understand today’s pandemic scenario that the entire human race is suffering from. To tackle the destructive COVID-19, many remarkable developments in the vaccine domain are being made to obtain the maximum number of people vaccinated and reduce the infectiousness of the whole population. Because of the global widespread of infection, there has been a great demand for vaccination. Countries all over the world struggle to implement vaccination schemes against SARS-CoV-2 with the hurdles like limited production of doses or mismanagement in planning vaccine acquirement and implementations, which further hamper the process to control the infection. Our considered model is very well relatable in this scenario where almost every individual is desiring for vaccination. It is also the answer to how retrial queueing models can be applied in solving congestion issues of various vaccination centers in hospitals.

Mathematical modeling is playing an important part in understanding the dynamics of Covid-19, which further gives an insight in determining the best strategies to undertake to control the spread of infection [3,18,19]. With the onset of pandemics, the focus has been shifted to SIR (susceptible-infected-recovered) models consisting of vaccination classes that may either contain different types of interventions or even consider time delays related to latency, immunity, or individuals responses to public health-related information [2,14, 20,21]. For rota-virus incorporating time delay, a model of delay differential equation and the effects of vaccination are discussed in [1]. A model consisting of delayed vaccination regimens in [10] is analyzed to discuss the better strategies required for vaccination campaigns. The SIR model with pulse vaccination is discussed in [16], and for impulsive differential equations with delay, they have used new computational techniques in it. A model in [5] is discussed, which implied that the current delay in the vaccination schedules can have serious consequences in terms of mortality. The effect of vaccination in the disease progression in Covid-19 has also been modeled and discussed with a delay term in [7]. But there is a need to understand the effect of vaccines and delays related to them in the vaccination system as well. Hence, emphasizing on the need to shift focuses on the management to avoid delays and fast inoculation processes.

The queueing theory is a powerful tool that is used while resolving the complications occurring in the telecommunication and networking sector. However, recently, it is also applied in the medical field so that the health care resources are well utilized along with minimal waiting time for the patients. An appropriate outpatient queueing model for proper appointment scheduling was studied in [17]. Their study is helpful in reducing the waiting time of patients and doctors’ idle time as well. Ref. [4] designed the flow of various classes of patients into a hospital and analyzed the same. They classified the patients into service groups and chose different service policies for sequencing the patients through hospital units. Recently, [27] employed queueing theory in managing the ventilator capacity of a hospital under different COVID-19 scenarios. Another remarkable work done is by [8]. They applied transient queueing analysis for managing emergency circumstances in a hospital.

Retrial queues, a significant segment of queueing theory involve those queueing systems wherein if the service provider is not accessible on the arrival of a customer, then that customer joins a virtual pool (called the orbit) and retries from there to avail the service. A lot of research is carried out in this sphere to sort out the congestion and economic issues of analyzers in information transmitting systems (cf. [15,23,24]). Healthcare models formed by amalgamating such types of queues can be quite helpful for medical experts. Ref. [11] carried out performance analysis of cloud computing in health care systems via multi-server queues and single-server retrial queues. The significance of recurrent systems for healthcare patient flow operations management is very well depicted by [26]. Further, [6] has explained how an unreliable priority retrial model is suitable for a system with a single doctor wherein ordinary as well as emergency patients arrive for the treatment.

In our investigation, we have considered an infestation retrial queueing model involving healthy individuals. The entire population is susceptible to infection, which can only be prevented by taking both doses of vaccination. As the population size is huge, all individuals cannot get the vaccination slot at the first attempt. Thus, they need to retry to avail of the service of vaccination. The Runge-Kutta method of fourth order (RK4) is employed to obtain the numerical results. Although some other methods like the Laplace transform can also be employed for solving the differential–difference equations, we have used RK4 method because of its high convergence. Many queueing researchers are admiring this method now as it is unlike multi-step complicated methods and for the fact that it is economically friendly. A model was analyzed by [22] with alert, infection, and vaccination states via RK4 method. This method was also applied by [13] wherein they carried out transient performance analysis of a single server Markovian queueing model with balking and discussed the impact of the model on the health care queueing management system.

The individuality of our work lies in the concept of developing our model with delay and involving retrial phenomenon for booking the slots of vaccination through queuing model. In the dynamic model whenever delay is taking place, then there is a time lag. During that time lag, the concept of retrial is applicable. So, we have incorporated queueing model too in our work, and thus, both the considered models are inter-related. This makes our work more practical in today’s scenario, and as observed, there is no such work done till now.

2 Dynamical model

We have studied our model in five states, namely, the healthy state, dose 1 of vaccination, after effects stage, dose 2 of vaccination, and the recovery stage. All the individuals try to book slots for vaccination, but there is a delay in doing so which leads to the concept of a retrial queueing system. We shall now formulate a mathematical model to study the consequences of delay in the vaccination schedules due to the inability to book a slot for vaccination (dose one and dose two). The model is formulated as per the following assumptions:

H is the healthy population going for vaccination during the pandemic to attain immunity.
V 1 are the healthy people vaccinated with the first dose, V A are people with mild symptoms, and V 2 are people vaccinated with the second dose.
The system consists of healthy and vaccinated individuals only and not infected individuals to keep the focus on the vaccination system during the pandemic.
R is the class consisting of the population who have successfully completed their vaccination and attained immunity against the disease.
Then two doses of vaccinations have been considered in context to countries like India where initially mainly two-dose vaccines (Covishield/AstraZeneca) were promoted by responsible stakeholders with a huge gap of 9 months between both doses.
We assume non-fatal implications of people going for vaccination (one dose, two doses, mild symptoms) or negligible death rates and considered the death rate of only healthy and recovered individuals.
Individuals after getting the first dose of the vaccine develop mild symptoms.
There is a delay in getting vaccine slots for both doses ( τ 1 , τ 2 ), which also gave rise to the retrial queues by people in between this delay time.
The system we are assuming is a vaccination system consisting of healthy and vaccinated individuals (compartments) only to see how the dynamics of the system change with change in this delay time.

The dynamical model for the system is given by:

(1) d H d t = ξ − ψ 1 H ( t − τ 1 ) − μ H d V 1 d t = ψ 1 H ( t − τ 1 ) − α V 1 d V A d t = α V 1 − ψ 2 V A ( t − τ 2 ) d V 2 d t = ψ 2 V A ( t − τ 2 ) − β V 2 d R d t = β V 2 − μ R

whose schematic diagram is shown in Figure 1, and the variables and parameters are defined in Table 1.

Figure 1

Schematic diagram of model.

Table 1

Parameters/variables meaning

Parameters/variables	Meaning
ξ	Intrinsic growth rate due to birth or immigration
ψ 1	Rate of healthy population getting vaccinated by first dose
μ	Natural death rate of healthy people
α	Rate at which people are joining the class V A due to after effects
ψ 2	Rate of V A class getting vaccinated by second dose after recovering
β	Recovery rate after attaining immunity after two doses
τ 1	Delay in vaccination by first dose
τ 2	Delay in vaccination by second dose

The system (1) has the non-trivial equilibrium point E ∗ ( H ∗ , V 1 ∗ , V A ∗ , V 2 ∗ , R ∗ ) , where:

H ∗ = ξ ψ 1 + μ , V 1 ∗ = ψ 1 ξ α ( ψ 1 + μ ) , V A ∗ = ψ 1 ξ ψ 2 ( ψ 1 + μ ) , V 2 ∗ = ψ 1 ξ β ( ψ 1 + μ ) , R ∗ = ψ 1 ξ μ ( ψ 1 + μ ) .

It is easy to see from the equilibrium points that people will remain in the healthy state for the time 1 ψ 1 + μ . ψ 1 ξ α ( ψ 1 + μ ) healthy people get vaccinated per unit time by first dose. In the same manner, ψ 1 ξ ψ 2 ( ψ 1 + μ ) are first dose vaccinated people showing mild symptoms per unit time, ψ 1 ξ β ( ψ 1 + μ ) people are second dose vaccinated people, and ψ 1 ξ μ ( ψ 1 + μ ) are the recovered people (immune people) after both the doses per unit time. It may be also noted that if the vaccination rate (first dose) ψ 1 or death rate of healthy people μ is sufficiently large or if intrinsic growth rate ξ is sufficiently small, then the non-trivial equilibrium point E ∗ will not exist.

2.1 Stability and local Hopf bifurcation analysis

We shall discuss the stability and Hopf bifurcation [12,25] of the positive equilibrium point E ∗ of system (1) as follows:

The community matrix of system at point E ∗ is given by J = A + B e − τ 1 λ + C e − τ 2 λ , where

A = − μ 0 0 0 0 0 − α 0 0 0 0 α 0 0 0 0 0 0 − β 0 0 0 0 β − μ , B = − ψ 1 0 0 0 0 ψ 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 , C = 0 0 0 0 0 0 0 0 0 0 0 0 − ψ 2 0 0 0 0 − ψ 2 0 0 0 0 0 0 0 .

The characteristic equation for the linearized system at the endemic equilibrium E ∗ for our model [9] is given as below:

(2) P ( λ ) + Q ( λ ) e − τ 1 λ + R ( λ ) e − τ 2 λ = 0 ,

where P ( λ ) = λ ( λ + μ ) ( λ + α ) ( λ + β ) ( λ + μ ) , Q ( λ ) = ( λ + ψ 1 ) λ 4 , R ( λ ) = ( λ + ψ 2 ) λ 4 .

Here, for the two different delays τ 1 and τ 2 , we have the following cases:

Case (i): τ 1 = 0 , τ 2 = 0

Theorem 1

E ∗ is locally asymptotically stable for τ 1 = 0 , τ 2 = 0 provided K 14 K 13 − K 12 > 0 , K 14 ( K 13 K 12 − K 11 K 14 ) > K 12 K 12 , K 14 ( K 13 K 12 − K 11 K 14 ) > K 12 K 12 , K 14 ( K 13 K 12 − K 11 K 14 ) > K 12 2 , and K 14 ( K 13 K 12 − K 11 2 ) > K 12 2 .

Proof

The reduced characteristic equation of (2) is given as follows:

λ 5 + K 14 λ 4 + K 13 λ 3 + K 12 λ 2 + K 11 λ = 0 ,

where K 14 = ( β + α + 2 μ + ψ 1 + ψ 2 ) 5 , K 13 = ( β + μ + β μ + μ α ) 5 , K 12 = ( β + μ + μ α ( β + μ ) ) 5 , and K 11 = μ 2 β α 5 .

By Routh–Hurwitz criterion Re ( λ ) < 0 if and only if

det 1 = K 14 1 K 12 K 13 > 0

provided K 14 K 13 − K 12 > 0 ,

det 2 = K 14 1 0 K 12 K 13 K 14 0 K 11 K 12 > 0

provided K 14 ( K 13 K 12 − K 11 K 14 ) > K 12 K 12 , and

det 3 = K 14 1 0 0 K 12 K 13 K 14 1 0 K 11 K 12 K 13 0 0 0 K 11 > 0

provided K 14 ( K 13 K 12 − K 11 K 14 ) > K 12 K 12 , K 14 ( K 13 K 12 − K 11 K 14 ) > K 12 2 , and

det 4 = K 14 1 0 0 0 K 12 K 13 K 14 1 0 0 K 11 K 12 K 13 K 14 0 0 0 K 11 K 12 0 0 0 0 K 11 > 0

provided K 14 ( K 13 K 12 − K 11 2 ) > K 12 2 .□

Case (ii): τ 1 ≠ 0 , τ 2 = 0

Theorem 2

For τ 1 ≠ 0 , τ 2 = 0 , the endemic equilibrium E ∗ is locally asymptotically stable for the system when τ 1 ∈ [ 0 , τ 10 ) (where τ 10 is a threshold) and the system undergoes a Hopf bifurcation at E ∗ when τ 1 = τ 10 a family of periodic solutions bifurcate from E ∗ .

Proof

The reduced characteristic equation of (2) is as follows:

(3) 2 λ 5 + ( β + α + 2 μ + ψ 1 ) λ 4 + ( β + μ + β μ + μ α ) λ 3 + ( β + μ + μ α ( β + μ ) ) λ 2 + λ μ 2 β α + e − τ 1 λ ( λ 5 + ψ 1 λ ) = 0 .

By considering λ = i ω , we obtain

(4) ( ω 5 + ψ 1 ω ) sin ( τ 1 ω ) = ω 2 ( β + μ + μ α ( β + μ ) ) − ω 4 ( β + α + 2 μ + ψ 1 ) ( ω 5 + ψ 1 ω ) cos ( τ 1 ω ) = ω 3 ( β + μ + β μ + μ α ) − ω μ 2 β α − 2 ω 2 ,

through which we obtain

(5) 3 ω 10 + M 1 ω 8 + M 2 ω 6 + M 3 ω 4 + M 4 ω 2 = 0 ,

where M 1 = − ( β + α + 2 μ + ψ 1 ) 2 − 4 , M 2 = 2 ψ 1 − 2 ( β + α + 2 μ + ψ 1 ) ( β + μ + μ α ( β + μ ) ) + ( β + μ + β μ + μ α ) 2 + 4 μ 2 β α , M 3 = ( β + μ + μ α ( β + μ ) ) 2 − 2 ( β + μ + β μ + μ α ) μ 2 β α , M 4 = ψ 1 2 + ( μ 2 β α ) 2 . Let v = ω 2 , so we obtain

(6) 3 v 5 + M 1 v 4 + M 2 v 3 + M 3 v 2 + M 4 v = 0

J 21 : equation (6) has a positive root v 0 . Then, equation (5) has a positive root ω 0 = v 0 . We thus obtain:

τ 10 = 1 ω 0 × arccos M 21 ( ω 0 ) M 22 ( ω 0 ) ,

where M 21 ( ω 0 ) = ( β + μ + β μ + μ α ) ω 3 − μ 2 β α ω − 2 ω 5 and M 22 ( ω 0 ) = ω 5 + ω ψ 1 . Now we shall differentiate equation (3) with respect to τ 1 and putting λ = i ω 10 , and we obtain:

Re d λ d τ 1 λ = i ω 10 − 1 = g 1 ′ ( v 10 ) M 22 ( ω 10 2 ) .

Hence, Re d λ d τ 1 λ = i ω 10 − 1 ≠ 0 if the condition J 22 : g 1 ′ ( v 10 ) ≠ 0 holds, where g 1 ( v 0 ) = 3 v 5 + M 1 v 4 + M 2 v 3 + M 3 v 2 + M 4 v . Based on the Hopf bifurcation theorem, system undergoes a Hopf bifurcation when conditions J 21 and J 22 hold.□

Case (iii): τ 1 = 0 , τ 2 ≠ 0

Theorem 3

For τ 2 ≠ 0 , τ 1 = 0 , the endemic equilibrium E ∗ is locally asymptotically stable for the system when τ 2 ∈ [ 0 , τ 20 ) (where τ 20 is a threshold), and the system undergoes a Hopf bifurcation at E ∗ when τ 2 = τ 20 a family of periodic solutions bifurcate from E ∗ .

Proof same as above.

Case (iv): τ 1 ≠ 0 , τ 2 ≠ 0

In this case for τ 1 > 0 and τ 2 ∈ [ 0 , τ 20 ) (where τ 20 is a threshold), we shall consider τ 1 as bifurcation parameter and τ 2 in stable interval and consider the characteristic equation as follows:

(7) λ 5 + I λ 4 + J λ 3 + K λ 2 + L λ + e − τ 1 λ ( λ 5 + ψ 1 λ ) + e − τ 2 λ ( λ 5 + ψ 2 λ ) = 0 ,

where I = ( β + α + 2 μ ) , J = ( β + μ + β μ + μ α ) , K = ( β + μ + μ α ( β + μ ) ) , and L = μ 2 β α . Now let λ = i ω ∗ as a root of equation (7), we obtain:

(8) ω ∗ 5 sin ( τ 1 ω ∗ ) + ψ 1 ω ∗ 4 cos ( τ 1 ω ∗ ) = M 31 , ω ∗ 5 cos ( τ 1 ω ∗ ) + ψ 1 ω ∗ 4 sin ( τ 1 ω ∗ ) = M 32 ,

where M 31 = I ω ∗ 4 − K ω ∗ 2 + ω ∗ 5 sin ( τ 2 ω ∗ ) + ψ 2 ω ∗ 4 cos ( τ 2 ω ∗ ) and M 32 = ω ∗ 5 − J ω ∗ 3 + L ω ∗ + ω ∗ 5 cos ( τ 2 ω ∗ ) − ψ 2 ω ∗ 4 sin ( τ 2 ω ∗ )

H 31 : equation (7) has a positive root v ∗ . Then, equation (5) has a positive root ω ∗ = v ∗ . We thus obtain:

τ 1 = 1 ω ∗ × arccos ψ 1 ω ∗ 4 M 31 + ω ∗ 5 M 32 ( ψ 1 ω ∗ 4 ) 2 + ω ∗ 10 .

Now we shall differentiate equation (7) with respect to τ 1 and putting λ = i ω ∗ , we obtain:

Re d λ d τ 1 λ = i ω ∗ − 1 = g ( ω ∗ ) ( ψ ω ∗ 4 ) 2 + ( ω ∗ 5 ) 10 ,

where g ( ω ∗ ) = − 4 ω ∗ 3 I + 2 K ω ∗ − 4 ψ 2 ω ∗ 2 cos ( ω ∗ τ 2 ) + 5 ω ∗ 4 sin ( ω ∗ τ 2 ) − τ 2 ω ∗ 4 cos ( ω ∗ τ 2 ) − τ 2 ψ 2 ω ∗ 4 sin ( ω ∗ τ 2 ) + 5 ω ∗ 3 sin ( ω ∗ τ 1 ) − ψ 1 4 ω ∗ 2 cos ( ω ∗ τ 1 ) − τ 1 ω ∗ 4 cos ( ω ∗ τ 1 ) + τ 1 ψ 1 ω ∗ 3 sin ( ω ∗ τ 1 ) . Thus, Re d λ d τ 1 λ = i ω ∗ − 1 ≠ 0 if the condition H 22 : g ( v ∗ ) ≠ 0 holds. Based on the Hopf bifurcation theorem, system undergoes a Hopf bifurcation when conditions H 21 and H 22 hold.

In the next section, we would propose a retrial queuing model to show that in the delay time, people will keep trying, again and again, to book their slots for both doses.

3 Queueing model

For the particular queueing model, we have considered a population size of K individuals. All healthy individuals try to take up vaccination, which comprises two phases, namely, dose 1 and dose 2. As the population size is huge, so it is difficult for them to book the vaccination slot on the first attempt. Thus, they retry with a certain retrial rate for taking up the vaccine. Further, once the first dose is completed, there are certain after-effects such as slight fever, vomiting, etc. After overcoming these effects, the individuals try for the second dose of vaccination with the same retrial rate. This shows that there is a time lag between two doses of vaccination; due to this retrial, the phenomenon can be seen in the same model whose dynamical analysis is discussed in Section 2. Each person is supposed to be completely recovered if they avail both doses of vaccination. So we will use the same definitions for a few symbols defined in Section 2 for a dynamical model, which are H , V 1 , V 2 , and R , but now here, we will use these for queueing modelling in transient state only. Further, we will define the following assumptions for queueing model:

There is a certain birth rate ( ξ ) and death rate ( μ ) in healthy/recovery state.
Healthy persons retry to get vaccinated with retrial rate γ .
The rate of arrival of individuals for vaccination of dose 1 and dose 2, are respectively, ϕ 1 and ϕ 2 ( v and v ′ in developed retrial queueing model).
After the first dose of vaccination, there are certain aftereffects. The transition rate of the vaccination (dose 1) class to after effects class is denoted by α .
After overcoming the consequences of dose 1, individuals attempt for the second dose with the same retrial rate.
The entire population recovers with rate β after getting a second dose of vaccination.
The vaccination completion rates for dose 1 and dose 2 are w and w ′ , respectively.

The states of the system (for j = 1 , 2 , 3 , … , K ) are defined as follows:

( H , j ): The state when there are j number of healthy individuals in the system.

( V 1 , j ): The state when j number of healthy individuals are under the first dose of vaccination.

( V A , j ): The state when j number of individuals undergo after-effects of dose 1 of vaccination.

( V 2 , j ): The state when j number of individuals are under the second dose of vaccination after overcoming the after-effects of dose 1.

( R , j ): The state when j number of individuals are completely recovered.

We further assume that P k ( t ) be the transient state probability that system is in the kth state, at time t where, [ k = ( H , j ) , ( V 1 , j ) , ( V A , j ) , ( V 2 , j ) , ( R , j ) ; j = 1 , 2 , 3 , … , K ]. All the aforementioned assumptions are summarized in Figure 2. We can see that all five states are depicted in this figure. The dots here represent the population size which goes from 1 to k . After getting both doses of vaccination, the individual has a certain immunization and is thus considered as recovered. This figure also illustrates all the transition rates and the retrial rate among different states.

Figure 2

State transition diagram.

The transient state equations for the model described earlier are as follows:

(9) d P H , 1 ( t ) d t = − ξ P H , 1 ( t ) + μ P H , 1 ( t ) ,

(10) d P H , j ( t ) d t = − ( μ + γ + ξ ) P H , 1 ( t ) + ξ P H , j − 1 ( t ) + μ P H , j + 1 ( t ) ; 2 ≤ j ≤ K ,

(11) d P H , K ( t ) d t = ξ P H , K − 1 ( t ) − ( μ + γ ) P H , 2 ( t ) ,

(12) d P V 1 , 1 ( t ) d t = − ( v + α ) P V 1 , 1 ( t ) + w P V 1 , 2 ( t ) + γ P H , 2 ( t ) ,

(13) d P V 1 , j ( t ) d t = − ( v + α + w ) P V 1 , j ( t ) + w P V 1 , j + 1 ( t ) + v P V 1 , j − 1 ( t ) + γ P H , j + 1 ( t ) ; 2 ≤ j ≤ K ,

(14) d P V 1 , K ( t ) d t = − ( α + w ) P V 1 , K ( t ) + v P V 1 , K − 1 ( t ) ,

(15) d P V A , 1 ( t ) d t = α P V 1 , 1 ( t ) ,

(16) d P V A , j ( t ) d t = α P V 1 , j ( t ) − γ P V A , j ( t ) ,

(17) d P V 2 , 1 ( t ) d t = − ( V 2 + β ) P V 2 , 1 ( t ) + w ′ P V 2 , 2 ( t ) + γ P V A , 2 ( t ) ,

(18) d P V 2 , j ( t ) d t = − ( v ′ + β + w ′ ) P V 2 , j ( t ) + w ′ P V 2 , j + 1 ( t ) + v ′ P V 2 , j − 1 ( t ) + γ P V A , j + 1 ( t ) ; 2 ≤ j ≤ K ,

(19) d P V 2 , K ( t ) d t = − ( w ′ + β ) P V 2 , K ( t ) + v ′ P V 2 , K − 1 ( t ) ,

(20) d P R , j ( t ) d t = β P V 2 , j ( t ) ; 1 ≤ j ≤ K .

4 Performance measures

As we have obtained all transient state probabilities of system states using the Runge-Kutta method of the fourth order, we can now obtain various fruitful performance metrics to apply it in any healthcare system with a single server. A queueing model is not considered to be complete if the system manager is not aware of its performance as what queue length or waiting time does it anticipate? For that purpose, performance indices play a major role and help to build up a model which can satisfy the needs of service providers and customers as well. Thus, the eminent performance indices of this model are obtained via the Runge-Kutta method of fourth order, which is implemented via the MATLAB’s “ode45” function. Once the performance indices are acquired, then only we can numerically study the model. Some of the important metrics of the developed queueing model are discussed here:

The expected number of individuals in a healthy state is expressed as follows:
(21) E [ H ] = ∑ j = 1 K j P H , j ( t ) .
The expected number of individuals under vaccination dose 1 state is expressed as follows:
(22) E [ V 1 ] = ∑ j = 1 K j P V 1 , j ( t ) .
The expected number of individuals under the after effects state is expressed as follows:
(23) E [ V A ] = ∑ j = 1 K j P V A , j ( t ) .
The expected number of individuals under vaccination dose 2 state is expressed as follows:
(24) E [ V 2 ] = ∑ j = 1 K j P V 2 , j ( t ) .
The expected number of individuals under recovery is expressed as follows:
(25) E [ R ] = ∑ j = 1 K j P R , j ( t ) .

5 Numerical discussion

We shall introduce some numerical simulations by means of Matlab in this section, which would illustrate our obtained theoretical results for both the developed dynamical model and the queueing model and have used the parametric values (hypothetical) as mentioned in Table 2 in both cases.

Table 2

Parameter values

Parameter	Value
ξ	10
ψ 1	0.9
μ	0.09
α	0.2
ψ 2	0.09
β	0.09
γ	0.55
w	0.28
w ′	0.28

5.1 Numerical simulations for delayed ODE model

In this subsection, we would discuss the numerical simulation of our delay-based ODE model to show the effect of delay in obtaining the first and second doses of vaccination.

In Figure 3, we took the delay τ 1 = 43 ( τ 1 ≠ 0 , τ 2 = 0 ), which is an example of a stable behavior of the system for case (ii). The system remains stable for τ 1 ∈ [ 0 , τ 10 = 80 ) . At τ 10 = 80 , the system starts showing minor periodic solutions as shown in Figure 4. All the graphs are plotted at initial conditions ( 100 , 40 , 40 , 30 , 40 ) . To validate the theoretical results, we have also obtained the equation 3 v 4 − 5.8769 v 3 − 101.2404 v 2 + 0.0335 v + 0.81 = 0 , which is in line with Theorem 2 and obtain a positive root 0.0894 for which we obtain our threshold τ 10 . The availability of the first dose of vaccine is the first step towards ensuring a disease-free system. Hence, the first dose of vaccine should be made available within a certain time period so that it can ensure a system with maximum vaccinated individuals without delay (Figure 5).
In the presence of both the delays ( τ 1 ≠ 0 , τ 2 ≠ 0 ), now we keep τ 1 = 70 > 0 , which is in the stable region and start variating τ 2 . At τ 2 = 1 , we obtain stable dynamics as shown in Figure 6. Now, as we keep increasing τ 2 , the system remains stable till we obtain hopf bifurcation for τ 1 = 79 and τ 2 = 120 . With regard to the aforementioned observations, we can conclude that the vaccination delay (both doses) is a key parameter in the dynamic of our model and can even destabilize the number of vaccinated and non-vaccinated people. Thus, a delay in getting both doses can proceed to make the system unstable, which can be indicative of an uncertain vaccination drive.

Vaccines (doses one and two) are launched with the aim to make the system stable, but a delay in getting the doses creates an instability in the number of vaccinated people with both doses and consequently can hamper the vaccination strategy and the maximum benefit of the drive. The periodic behavior can lead to extremely fluctuating values for the population during the vaccination process, which would cause problems in forming future strategies for vaccination.

$Figure 3 Stable dynamics of the system at τ 1 = 43 {\tau }_{1}=43 .$

Figure 3

Stable dynamics of the system at τ 1 = 43 .

$Figure 4 Dynamics of the system at τ 1 = 80 {\tau }_{1}=80 .$

Figure 4

Dynamics of the system at τ 1 = 80 .

$Figure 5 Dynamics of the system at τ 1 = 80 {\tau }_{1}=80 and τ 2 ∈ [ 0 , τ 20 ) {\tau }_{2}\in \left[0,{\tau }_{20}) , where τ 20 = 120 {\tau }_{20}=120 .$

Figure 5

Dynamics of the system at τ 1 = 80 and τ 2 ∈ [ 0 , τ 20 ) , where τ 20 = 120 .

$Figure 6 Stable dynamics of the system at τ 1 = 80 {\tau }_{1}=80 and τ 2 = 1 {\tau }_{2}=1 .$

Figure 6

Stable dynamics of the system at τ 1 = 80 and τ 2 = 1 .

5.2 Numerical illustration for developed queuing model

In this section, we have numerically verified the results provided for the queuing model. The expected value of individuals in different states are studied, and hence, the respective graphs are plotted with respect to equally slotted time interval. The Runge-Kutta (RK) method in Section 5.2 of the fourth order is employed for this purpose and is implemented via MATLAB’s “ode45” function. The default values of the system parameters are the same as provided in Table 2. The remaining parameters have default values γ = 0.55 , w = w ′ = 0.28 (c.f. we have referred to [22] for the values of queueing parameters).

We have studied the trend of the expected number of individuals in the vaccination dose 1 state, vaccination dose 2 state, and recovery state with respect to time. We observe in Figure 7 that as there is an increase in v and v ′ , then E [ V 1 ] and E [ V 2 ] also rise, respectively. This matches with the real-life scenario that a rise in the arrival rate for vaccination implies a rise in the number of individuals in that respective state. Further, in Figure 8, the trend of E [ V 2 ] and E [ R ] is observed by varying the retrial rate. We see that a rise in γ implies a rise in E [ V 2 ] . The behavior of these graphs are studied on a longer simulation over time. Further, we can conclude that when the value of γ increases, then it directly increases E [ V 2 ] and E [ R ] . It is quite obvious that since more people retrying for the vaccination implies that more individuals will avail the vaccine and hence will get recovered. We can say that the trend of E [ V 2 ] and E [ R ] remains same via varying retrial rate even if we increase the time span. So, the congestion issues in distinct vaccination centers can be reduced when the queueing parameters are chosen wisely. Other health care departments, for instance, blood testing or ultrasound department can also reduce the waiting time of their patients via a suitable retrial queueing model. Overall we conclude that the number of people vaccinated depend directly on the retrial rate, which depicts that as the willingness in getting vaccinated by the people increases (by increasing retrial rate), more people get vaccinated which in turn helps in vaccinating a larger population.

$Figure 7 (a) Effect of v v on E [ V 1 ] E\left[{V}_{1}] . (b) Effect of v ′ v^{\prime} on E [ V 2 ] E\left[{V}_{2}] .$

Figure 7

(a) Effect of v on E [ V 1 ] . (b) Effect of v ′ on E [ V 2 ] .

$Figure 8 (a) Effect of γ \gamma on E [ V 2 ] E\left[{V}_{2}] . (b) Effect of γ \gamma on E [ R ] E\left[R] .$

Figure 8

(a) Effect of γ on E [ V 2 ] . (b) Effect of γ on E [ R ] .

6 Conclusion

To help guide decision-making, we propose here a novel model and analyze the dynamics in regard to the first and second vaccine dose delays/retrials. To prevent the collapse of the complex inoculation system, the vaccination of the general public must be done efficiently. A queuing model with five states, namely, healthy, vaccination (dose 1 and dose 2), after effects, and recovery is analyzed. We found the equilibrium point E ∗ for our ODE system. For local stability, we have obtained sufficient conditions by analyzing the associated characteristic equation. Then, the governing equations are formed, and the fruitful performance measures are acquired by applying the RK method of the fourth order. Then, through numerical simulations, we saw that with the increase of delay ( τ 1 = 80 ) for getting a slot for the first dose of vaccine, the system becomes more unstable. And with the delay ( τ 1 = 80 and τ 2 = 120 ) in getting the second dose of vaccine as well, we saw the system become unstable and showed hopf bifurcation. Thus due to the delay of both doses (one and two), it can make the system go to instability rather than a more ideal stable system in order to suppress the spread of the disease. This should serve as an alert to health authorities to speed up the vaccination process so that a maximum number of people would be immunized in the shortest period of time, i.e., delay be as less as possible for a stable system. For the delay in getting a vaccination slot, a retrial by an individual for booking the slot can help stabilize the system. That is an increase in the rate of retrial for slot booking can help in reducing the delay in the system and ensures a timely vaccination drive. We also observe that the system parameters directly influence the expected queue size of the different stages considered in the queuing model. These results are also verified with the help of numerical illustrations. We can conclude that the retrial rate directly influences the expected number of individuals in the vaccination (dose 2) state and recovery state as well. Both these states have more individuals if the retrial rate is increased. Thus, choosing an apt rate of retrial can lead us to the maximum number of healthy individuals. This study can be helpful in solving congestion issues in the medical field as a greater number of people can be vaccinated in less time. We can say that the battle against a pandemic can be won if appropriate system parameters of the queuing model are chosen. Responsible stakeholders should make sure that successful and timely vaccination of people is achieved through rigorous vigilance of registration and vaccine booking process. No delay in getting a vaccination slot will stimulate a successful vaccination drive and further motivate especially the government of India to vaccinate more people mainly in rural areas and children below 18 years of age to get vaccinated in regard to the easy vaccination process. In the future, we can work on the effect of the delay in the booster dose vaccination process started in India (but not at the proper rate) in reducing the congestion faced by healthcare systems.

# All authors contributed equally to this work.

Funding information: This research received no specific grant from any funding agency, commercial or non-profit section.
Conflict of interest: All authors declare no conflicts of interest in this paper.
Ethical approval: The conducted research is not related to either human or animal use.
Data availability statement: Data sharing is not applicable to this article as no datasets were generated or analyzed during the current study.

References

[1] Adongo, F. A., Lawrence, O. O., Bonyo, J., Lawi, G., & Elisha, O. A. (2020). A delayed vaccination model for rotavirus infection. European Journal of Pure and Applied Mathematics, 13(4), 840–851. 10.29020/nybg.ejpam.v13i4.3822Search in Google Scholar

[2] Agaba, G., Kyrychko, Y., & Blyuss, K. (2017). Dynamics of vaccination in a time-delayed epidemic model with awareness. Mathematical Biosciences, 294, 92–99. 10.1016/j.mbs.2017.09.007Search in Google Scholar PubMed

[3] Akman, O., Chauhan, S., Ghosh, A., Liesman, S., Michael, E., Mubayi, A., … , & Tripathi, J. P. (2022). The hard lessons and shifting modeling trends of covid-19 dynamics: Multiresolution modeling approach. Bulletin of Mathematical Biology, 84(1), 1–30. 10.1007/s11538-021-00959-4Search in Google Scholar PubMed PubMed Central

[4] Alenany, E., & El-Bas, M. A. (2017). Modelling a hospital as a queueing network: analysis for improving queueing performance. International Journal of Industrial and Manufacturing Engineering, 11, 1181–1187. Search in Google Scholar

[5] Amaku, M., Covas, D. T., Coutinho, F. A. B., Azevedo, R. S., & Massad, E. (2021). Modelling the impact of delaying vaccination against sars-cov-2 assuming unlimited vaccine supply. Theoretical Biology and Medical Modelling, 18(1), 1–11. 10.1186/s12976-021-00143-0Search in Google Scholar PubMed PubMed Central

[6] Bhagat, A. (2020). Unreliable priority retrial queues with double orbits and discouraged customers. Advancements in Mathematics and its Emerging Areas, 2214. 10.1063/5.0003372Search in Google Scholar

[7] Caga-Anan, R. L., Raza, M. N., Labrador, G. S. G., Metillo, E. B., delCastillo, P., & Mammeri, Y. (2021). Effect of vaccination to covid-19 disease progression and herd immunity. Computational and Mathematical Biophysics, 9(1), 262–272. 10.1515/cmb-2020-0127Search in Google Scholar

[8] Curry, G. L., Moya, H., Erraguntla, M., & Banerjee, A. (2021). Transient queueing analysis for emergency hospital management. IISE Transactions on Healthcare Systems Engineering, 12(1), 36–51. 10.1080/24725579.2021.1933655Search in Google Scholar

[9] Erman, S., Sevindir, H. K., & Demir, A. (2018). The stability analysis of a system with two delays. Boundary Value Problems, 2018(1), 1–11. 10.1186/s13661-018-1024-9Search in Google Scholar

[10] Gonzalez-Parra, G. (2021). Analysis of delayed vaccination regimens: A mathematical modeling approach. Epidemiologia, 2(3), 271–293. 10.3390/epidemiologia2030021Search in Google Scholar PubMed PubMed Central

[11] Kannan, S., & Ramakrishnan, S. (2016). Performance analysis of cloud computing in healthcare system using tandem queues. International Journal of Intelligent Engineering and Systems, 10, 256–264. 10.22266/ijies2017.0831.27Search in Google Scholar

[12] Karim, F., Chauhan, S., & Dhar, J. (2020). On the comparative analysis of linear and nonlinear business cycle model: Effect on system dynamics, economy and policy making in general. Quantitative Finance and Economics, 4, 172–203. 10.3934/QFE.2020008Search in Google Scholar

[13] Kumar, R., Soodan, B. S., & Sharma, S. (2021). Modelling health care queue management system facing patientas impatience using queueing theory. Reliability Theory and Applications, 16, 161–175. Search in Google Scholar

[14] Laarabi, H., Abta, A., & Hattaf, K. (2015). Optimal control of a delayed sirs epidemic model with vaccination and treatment. Acta Biotheoretica, 63, 87–97. 10.1007/s10441-015-9244-1Search in Google Scholar PubMed

[15] Malik, G., Upadhyaya, S., & Sharma, R. (2021). Particle swarm optimization and maximum entropy results for mx∕g∕1 retrial g-queue with delayed repair. International Journal of Mathematical, Engineering and Management Sciences, 6(2), 541–563. 10.33889/IJMEMS.2021.6.2.033Search in Google Scholar

[16] Meng, X., Chen, L., & Wu, B. (2010). A delay sir epidemic model with pulse vaccination and incubation times. Nonlinear Analysis: Real World Applications, 11(1), 88–98. 10.1016/j.nonrwa.2008.10.041Search in Google Scholar

[17] Obulor, R., & Eke, B. (2016). Outpatient queueing model development for hospital appointment system. International Journal of Scientific Engineering and Applied Sciences, 2, 15–22. Search in Google Scholar

[18] Rana, P., Chaudhary, K., & Chauhan, S. (2021). Dynamical analysis and optimal control problem of impact of vaccine awareness programs on epidemic system. Communications in Mathematical Biology and Neuroscience, 2021, 1–17. Search in Google Scholar

[19] Rana, P., Chauhan, S., & Mubayi, A. (2022). Burden of cytokines storm on prognosis of sars-cov-2 infection through immune response: dynamic analysis and optimal control with immunomodulatory therapy. The European Physical Journal Special Topics, 231, 1–19. 10.1140/epjs/s11734-022-00435-7Search in Google Scholar PubMed PubMed Central

[20] Rana, P., Jha, D., & Chauhan, S. (2022). Dynamical analysis on two dose vaccines in the presence of media. Journal of Computational Analysis & Applications, 30(2), 260–280.Search in Google Scholar

[21] Shim, E. (2011). Prioritization of delayed vaccination for pandemic influenza. Mathematical Biosciences and Engineering: MBE, 8(1), 95–112. 10.3934/mbe.2011.8.95Search in Google Scholar PubMed PubMed Central

[22] Singh, R. P., & Raina, A. (2018). Markovian epidemic queueing model with exposed infection and vaccination based on treatment. World Scientific News, 106, 141–150. Search in Google Scholar

[23] Upadhyaya, S. (2020). Investigating a general service retrial queue with damaging and licensed units: an application in local area networks. OPSEARCH, 57, 716–745. 10.1007/s12597-020-00440-1Search in Google Scholar

[24] Upadhyaya, S., & Khushwaha, C. (2020). Performance prediction and anfis computing for unreliable retrial queue with delayed repair under modified vacation policy. International Journal of Mathematics in Operational Research, 17(4), 437–466. 10.1504/IJMOR.2020.110843Search in Google Scholar

[25] Zhao, T., Zhang, Z., & Upadhyay, R. (2018). Delay-induced Hopf bifurcation of an sveir computer virus model with nonlinear incidence rate. Advances in Difference Equations, 256. 10.1186/s13662-018-1698-4Search in Google Scholar

[26] Zhu, H., & Liu, X. (2020). The analysis of feedback patient flow as a retrial queueing system in chinese outpatient department. Chinese Control and Decision Conference (pp. 3549–3552). Heifei, China: IEEE. 10.1109/CCDC49329.2020.9164178Search in Google Scholar

[27] Zimmerman, S., Rutherford, A. R., van der Waall, A., Norena, M., & Dodek, P. (2021). A queueing model for ventilator capacity management during the covid-19 pandemic. medRxiv. 10.1101/2021.03.17.21253488Search in Google Scholar

Received: 2022-03-05

Revised: 2023-03-15

Accepted: 2023-03-15

Published Online: 2023-05-22

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/cmb-2022-0147

Keywords for this article

delay; dynamic analysis; retrial queues; RK method of fourth order; vaccination process

Creative Commons

BY 4.0