Complexity analysis of cold chain transportation in a vaccine supply chain considering activity inspection and time-delay

The development of COVID-19 vaccine is highly concerned by all countries in the world. So far, many kinds of COVID-19 vaccines have entered phase III clinical trial. However, it is difficult to deliver COVID-19 vaccines efficiently and safely to the areas affected by the epidemic. This paper focuses on vaccine transportation in a supply chain model composed of one distributor and one retailer (clinic or hospital), in which the distributor procures COVID-19 vaccines from the manufacturer and then resells them to the retailer. Distributor detects the activity level of the vaccines, and retailer is responsible for transportation of the vaccines. Firstly, we establish a difference equations model with time-delay. Secondly, we investigate the impact of time-delay on the stability of vaccine supply chain. In addition, we explore the influence of decision adjustment speed of the distributor (or retailer) on the stability of vaccine supply chain. Finally, we verify the theoretical results by a two-dimensional bifurcation diagram, the largest Lyapunov exponent, entropy, and domain of attraction. The results show that when the decision delay-time or the adjustment speed of decision variables exceeds a certain threshold, it brings a negative impact on the stability of vaccine supply chain system. The stability domain of the system shrinks as customers’ sensitivity to cold chain transportation decreases and by contrast expends as customers’ sensitivity to vaccine prices decreases. When the vaccine supply chain is in a state of chaos, the effect of external control over the system is superior to that of internal control over the system.


Introduction
It is imperative to develop vaccines to prevent the rapidly spreading COVID-19 epidemic. Fortunately, up to now, many vaccine candidates of COVID-19 have been progressing at an unparalleled speed. 165 of them have been in the exploration or preclinical stage, 26 vaccine candidates have entered the clinical stage, and 6 of them with the fastest progress had even entered the clinical stage III. Three COVID-19 vaccines have been approved for clinical test III in China, among which Chen Wei's team focusses on the adenovirus vector vaccine.
It is very important how the vaccine is in time delivered to the patients infected with the virus. However, the transportation of vaccines is much more complicated than ordinary commodities. Transportation of the COVID-19 vaccines from the production plant to the final retailers within a few days is a complex system including storage factories, cargo stations, airplanes, warehouses, and so on. Therefore it is very urgent to build a safe global vaccine supply chain. Once COVID-19 vaccines are available, there must be a coordinated global strategy to ensure that they are delivered to millions of infected patients as quickly and safely as possible. Although the logistics industry is hard by the COVID-19 epidemic, the air cargo industries work hard to make the delivery process smoother. TIACA (The International Air Cargo Association) and the cross-industry cooperation platform Pharma.Aero have acted in advance to jointly compile guidelines for the global air cargo industries to realize the safe transportation of the COVID-19 vaccines. Based on the status quo that "the world is working hard to build a safe and efficient supply chain to ensure that the COVID-19 vaccines reach delivery place safely", the interesting conclusions obtained in this paper may provide some references for the transportation of the COVID-19 vaccines.
Vaccine supply chain has been studied from qualitative perspective. Evelot et al. [1] reviewed the relevant literature and found that the current research on vaccine supply chain mainly focuses on the following four aspects: vaccine quality, demand, allocation, and transport. Haidari et al. [2] found that an unmanned aerial system (UAS) to transport vaccines could improve the safety of vaccines and reduce costs. Veronica et al. [3] discovered that the inadvertent freezing of vaccine was an overlooked problem in the process of vaccine transportation. Sarley et al. [4] documented the work that the backward vaccine supply chain with Lagos State government was updated. Huang et al. [5] conducted a pilot experiment in which the health zone (HZ) was set to optimize the vaccine supply chain and assessed the incremental financial requirements for establishing a new system.
Many researchers have paid enormous attention to vaccine supply chain coordination from operation management perspective. Buyuktahtakm et al. [6] proposed a new epidemics-logistics mixed-integer programming (MIP) model to provide explicit intervention timing and intensity for these most affected countries. Abrahams et al. [7] presented a new binary integer programming model and a genetic algorithm to solve the complex scheduling problem in vaccine transportation. Arifoglu et al. [8] investigated the impact of self-interested consumers and yield uncertainty on the inefficiency of the influenza vaccine supply chain and revealed that government intervention could solve the inefficient problem of influenza vaccine supply chain. Cho [9] tried to determine the optimal composition of Influenza vaccines when vaccine production was uncertain. Dai et al. [10] investigated an influenza vaccine supply chain composed of a manufacturer and a retailer, and analyzed late deliveries and lost sales. Lin et al. [11] established a vaccine supply chain consisting of a distributor and a retailer, and discussed the impact of retailer's inspection at the end of transportation on the original decision-making of distributor. Niu et al. [12] constructed two representative transportation channel structures, established the game theory models, and derived the vaccine price equilibrium.
Some researchers have investigated in supply chain with methods of system dynamics. Ma and Xie [13] showed complex dynamic characteristics of the evolution process in a supply chain. Li et al. [14] analyzed the influence of different parameter values on the stability and utility of low-carbon supply chain system by means of a two-dimensional bifurcation diagram, parameter plot basin, the domain of attraction, and chaos attractor. Elsadany et al. [15] focused on the price and quantity competition in a mixed duopoly game and explored the dynamical behaviors of the models. Scholars have also found similar phenomena in the vaccine supply chain. Duijzer et al. [16] used the SIR model and nonlinear dynamics of epidemic to formulate disease progression and analyzed the best time for vaccination. Duijzer et al. [17] established a differential equation of epidemic time course and analyzed the relation between the "herd effect" (when people may escape infection without being vaccinated) and vaccination fraction. Dushoff et al. [18] established an epidemic mathematical model and found that a small change of the parameter value in the vaccine supply chain would make a huge change to the optimal vaccination strategy. The existing literature mainly focuses on the analysis of complex dynamic characteristics of vaccination rate in a supply chain. However, there is little literature on complex analysis of vaccine transportation in a supply chain.
In this paper, we investigate cold transportation in a vaccine supply chain consisting of a distributor and a retailer, who are of bounded rationality. The entropy theory is used to analyze the complex dynamic characteristics of vaccine supply chain according to different parameters, such as sensitivity coefficient of vaccine retail price, sensitivity coefficient of cold chain transportation, and the adjustment speed of make decision variables. This paper is organized as follows. In Sect. 2, we make description and assumptions of the problem. In Sect. 3, we formulate the problem discussed. In Sect. 4, we obtain the conditions for the bifurcation of the system caused by time delay. In Sect. 5, we discuss in detail the simulation results. In Sect. 6, we analyze control of the chaotic system, and Sect. 7 concludes the paper.

Description and assumptions of the problem 2.1 Model construction and assumptions
In this paper, we mainly consider cold transportation in a vaccine supply chain consisting of a distributor and a retailer in which the for-profit distributor procures vaccines from the pharmaceutical manufacturer and then resells them to the retailer. To ensure effective vaccination for every patient infected, distributor must carry out quality sampling of vaccines to be delivered, and the retailer is responsible for transportation of the vaccines and retails them to the infected patients. Because vaccines are sensitive to temperature, retailer must choose a cold chain to transport them.
Our main assumptions are as follows: (i) The distributor and the retailer are bounded rational.
(ii) We only consider one-time investment. If the quality sampling level for vaccines is y, then the distributor incurs the investment cost k 1 y 2 2 . If cold chain transportation level is b, then the retailer incurs the investment cost k 2 b 2 2 , where k 1 and k 2 are cost coefficients [19]. (iii) To ensure that the distributor and retailer can make normal profits, let p > w > c.

Symbolic description
The meanings of β, p, w, A, B, y, b, and c are described concisely in Table 1. Cold chain transport sensitivity. y The vaccine activity inspection level. b The cold chain transportation level.
The functional form of market demand can be written as follows:

Multiperiod decision-making game model with delay-time
The profit functions of the distributor and retailer are expressed as follows, respectively: From Eqs. (2) and (3), the decision variable of the distributor is the vaccine activity inspection level y, the retailer's decision variable is the cold chain transportation level b, the marginal profit function can be written as In the face of unfamiliar viruses, the research and development of vaccines need to invest much capital. Therefore the distributor and retailer must make an appropriate decision to reduce risk, that is, they adjust the next decision according to the current marginal profit. When the marginal profit of distributor or retailer is positive (or negative), this means that the increase of decision variables of distributor or retailer can increase (or reduce) the profit at the next period. Due to bounded rational distributor and retailer, the values of decision variables in period t + 1 are the values of decision variables at period t plus the changes of decision variables at period t [20]: where v 1 and v 2 are the adjustment speeds of decision variables of the distributor and retailer, respectively. Even though many scholars have investigated the vaccine supply chain, they assume that decision-makers make decisions instantaneously. However, this assumption in practice seems not realistic. The decision-makers cannot receive enough information in time when an epidemic breaks out. So we introduce the time-delay parameter into differential equations. The meaning of τ is the interval between the time when decision-maker should make a decision and the time when decision-maker takes a decision. We establish model I:

Existence and local stability of Neimark-Sacker bifurcation 4.1 Positive equilibrium points and characteristic equation of model I
By calculation we obtain two equilibrium points of model I: It is not hard to see that only E 2 is a positive equilibrium point. The point E 1 is nonpositive, and it may be an unstable equilibrium point. We can solve the Jacobian matrix of each equilibrium point to verify the observation.
We obtain the Jacobian matrix J of model I and judge its stability according to the magnitude of its eigenvalues. If all the eigenvalues of the Jacobian matrix are less than 1, then this equilibrium point is stable;otherwise, it is unstable.
The Jacobian matrix J of model I can be written as where The Jacobian matrix J 1 at the equilibrium point E 1 (0, 0) can be written as We can see that the eigenvalues of J 1 are not all less than 1, so the equilibrium point E 1 (0, 0) is an unstable equilibrium point. Similarly, it can be proved that only the equilibrium point E 2 is a stable equilibrium point of model I.
For brevity, let u 1 = y(t)y * , u 2 = b(t)b * , and we can transform the stability of the model I at the equilibrium point E 2 to the stability at the point (0,0), y(t) = u 1 + y * , b(t) = u 2 + b * . We use the Taylor expansion to expand Eq. (6) at the equilibrium point E 2 : where Next, the characteristic determinant of model I can be written as Then we obtain the characteristic equation of model I

Numerical simulation
In this section, we mainly verify the theoretical results. Set p = 4.8, w = 3.4, β = 0.26, c = 2, A = 3, B = 0.65, y = 2, b = 3, k 1 = k 2 = 0.5. We use the largest Lyapunov exponent and entropy to measure the features of system changes. The principle of the largest Lyapunov exponent is that the system is stable when the exponent is less than zero; otherwise, it is unstable. The rules of entropy to judge the system stability are as follows: the system is in a stable state when the entropy is zero; otherwise, it is unstable. Figures 1(a)-(c) show how the system stability changes with time-delay. As can be seen from Fig. 1(a), with the increase of τ , at first, Model I remains stable, then it gradually produces a Neimark-Sacker bifurcation, and, finally, it becomes chaotic. The critical point τ 0 of bifurcation is equal to 0.19. That is, when τ < 0.19, the equilibrium point E 2 (y * , b * ) is asymptotically stable; when τ increases and passes through τ 0 , the model I bifurcates and loses its stability.

Neimark-Sacker bifurcation diagram caused by time-delay
From Fig. 1(a) we can see that the vaccine transportation equilibrium point is locally asymptotically stable if the time-delay parameters is less than the critical value. When Figure 1 The impact of τ on the stability vaccine transportation is stable, the distributors and retailers can calmly respond to the epidemic and take reasonable decisions to ensure that the vaccine can reach its destination safely. However, once the time-delay parameter exceeds a certain threshold, the system undergoes a Neimark-Sacker bifurcation and goes into chaos. In this case, the distributor and retailer must bear the losses caused by the failure to respond to the epidemic in time. Figures 1(b) and (c) show the corresponding largest Lyapunov exponent and entropy of Fig. 1(a). When τ < 0.19, the largest Lyapunov exponent is less than zero, and the entropy is equal to zero, and the system keeps a stable state; otherwise, the largest Lyapunov exponent is greater than zero, the entropy of the system continues to increase, and the system falls into chaos state.

Bifurcation diagram caused by adjusting speed of decision variable
To further understand the dynamic characteristics of the system, in this section, we use the bifurcation diagram, largest Lyapunov exponent, and entropy to analyze the influence of the adjustment speed of decision variables on the system stability. Ceteris paribus, letting τ = 0.1, the influence of time-delay on the stability of the system can be eliminated. Assuming that v 2 = 0.015 and v 1 changes from 0 to 0.5, we analyze the effect of v 1 on the system stability. From Fig. 2(a) we see that when v 1 increases from 0 to 0.261, the vaccine activity inspection level y (indicated by yellow line) and the cold chain transportation level b (indicated by brown line) are asymptotically stable; when v 1 reaches 0.261, y and b start to bifurcate, and the system turns into stable cycles of period 2; as the v 1 continues to increase, the system eventually enters into a chaotic state. Figure 3(a) presents the effect of v 2 on the system stability, which is similar to Fig. 2(a). When v 2 ∈ [0, 0.146), the system is in a stable state; when v 2 increases to 0.146, the first bifurcation occurs, and the system turns into stable cycles of period 2; with the increase of v 2 , the system finally enters into chaotic state.
The largest Lyapunov exponent can further verify the dynamic characteristics of the system. Take Fig. 2(b) for an example: when the Lyapunov exponent first returns to the zero axis, the system appears the first bifurcation; when Lyapunov exponent is always greater than zero, then the system enters into a chaotic state, which is consistent with the dynamic characteristics shown in Fig. 2(a). Similarly, we also can illustrate the dynamic characteristics of the system through the change of entropy. From Fig. 2(c) we see that when v 1 ∈ (0, 0.261), the entropy of the model I is equal to zero, and at this time, y and b are asymptotically stable; when v 1 > 0.261,the entropy of model I is greater than zero, and the model I undergoes a period doubling bifurcation state; with the further increase of v 1 , the entropy also continues to increase. At this time, y and b are unstable and may take multiple possible values. Due to increase of the entropy, the distributor and retailer need more additional information to make reasonable decisions. Similar insights can be obtained in Figs. 3(b) and (c).
From Figs. 2-3 we see that the faster the adjustment speed of the vaccine activity inspection level (or the cold chain transportation level), the more chaotic the vaccine supply chain. From the perspective of entropy theory, when the vaccine supply chain falls into chaos, its entropy will be high. So the distributor and the retailer may obtain additional information to choose an appropriate adjustment speed. Therefore the distributor and retailer should collect abundant market information in advance and make rational decisions to prevent the vaccine supply chain from getting into chaos. Comparing Fig. 2(a) with Fig. 3(a), the critical value of system bifurcation caused by the adjustment speed v 2 of cold chain transportation level is less than that caused by the adjustment speed v 1 of activ- ity inspection level. This interesting phenomenon means that the reasonable adjustment range of activity inspection level is larger than that of cold chain transportation level, and consumers may be more sensitive to cold chain transportation level. Thus the improvement of the cold chain transportation level in the vaccine supply chain seems to be in the first place.

Global stability analysis
To further explore the dynamic characteristics of model I with the change of v 1 , v 2 , β, and B, ceteris paribus, letting y and b ∈ (0,1), we obtain that the domain of attraction of the model I is as in Fig. 4, in which the dark blue region denotes the stable attraction domain, light blue denotes the period-2 attraction region, and the deep red denotes the escape area. Figs. i (i = 5, 6,7,8) show the domain of attraction when v 1 , v 2 , β, B change in turn. Comparing Fig. 4 with Figs. 5 and 6, we can see that the stable attraction domains decrease with the increase of v 1 or v 2 , which is consistent with the conclusions reflected in Figs. 2 and 3. Comparing Fig. 4 with Fig. 7, the stable attraction domain decreases with the increase of β, which means that the lower the sensitivity of consumers to the vaccine price, the stabler the vaccine supply chain. Similarly, comparing Fig. 4 with Fig. 8, the stable attraction domain decreases with the decrease of B, which means that the higher the sensitivity of consumers to the cold chain transportation, the stabler the vaccine supply chain.
From Figs. 4-8 we see that the increase of consumers' sensitivity to vaccine prices has a negative impact on vaccine transportation. It is very necessary for the government to take measures to eliminate consumers' concerns about vaccine prices, such as providing subsi- dies to consumers. On the other hand, consumers' sensitivity to cold chain transportation is beneficial to the vaccine transportation system.

Chaos control
The vaccine transportation cannot be carried out smoothly by chaos in the vaccine supply chain, and the vaccines may be inactivated. Chaos in vaccine supply chain can be controlled by the adjustment parameter control method and variable feedback control method. The chaotic control effect of adjusting the parameters on model I is analyzed by numerical simulation. As mentioned before, ceteris paribus, when τ = 0.4, v 1 = 0.35, v 2 = 0.35, the vaccine supply chain is in chaos.

Adjustment parameter control method
The original vaccine supply chain system is The vaccine supply chain system after parameter adjustment control is as follows [23]: With the variation of the adjustment parameter u, the system changes as shown in Fig. 9. When u is less than the threshold (0.358), the system is in chaotic state, which indicates that the distributor and retailer did not take any effective joint measures to control chaos. With the increase of u, the system reaches a stable state, which indicates that the distributor and retailer can effectively reduce chaos by taking joint measures, such as signing contracts and so on.

Variable feedback control method
The main principle of variable feedback control method is using an equation variable as control signal to eliminate chaos. Compared with other control methods, this method has the advantages of simple controller design and strong timeliness. Therefore it is widely controllers is established as follows [24]: (31) Figure 10 shows that the chaotic system gradually returns to the stable state through variable feedback control method. When u is less than the threshold (0.118), the system is in chaotic state, which indicates that the government does not take effective control measures for the chaotic system at this time. When u > 0.118, the system remains asymptotically stable, which indicates that the government takes external intervention measures to accelerate the system to the stable state and ensures stable economic development.
Comparing Fig. 9 with Fig. 10, it is obvious that the control system in Fig. 10 enters the stable state earlier than that in Fig. 9, which means that the control effect of the variable feedback control method is better than that of the adjustment parameter control method.
Because of the extra cost of control, it is difficult for the distributor and retailer to actively control the chaotic vaccine supply chain. In addition, the safety of vaccines cannot be guaranteed in the chaotic vaccine supply chain, so external forces (for example, the regulation and policy of the government) usually are used to control the chaos of vaccine supply chain in time.

Conclusion
In this study, we mainly investigated the vaccine supply chain composed of a distributor and retailer. The distributor is responsible for the vaccine activity inspection of the vaccines before they are delivered to the retailer, whereas the retailer carries out the cold chain transportation of the vaccines. In addition, the decision-making of the distributor and retailer is not instantaneous, but rather time-delayed. We analyzed the complex dynamic characteristics of the system are by using the Neimark-Sacker bifurcation diagram, attraction domain, and entropy. We used the adjustment parameter control and variable feedback control methods to control the chaotic system. We obtained the following conclusions.
(1) When the time-delay parameter τ ∈ [0, τ 0 ), the vaccine transportation equilibrium is locally asymptotically stable. At this time, the distributor and retailer can cooperate well to ensure the vaccine activity. When the time-delay parameter τ ≥ τ 0 , the vaccines transportation system produces a Neimark-Sacker bifurcation and loses stability. (2) If the distributor and retailer adjust the decision variables too quickly, then the vaccine supply chain will bifurcate and even fall into chaos, which means that the distributor and retailer should choose appropriate adjustment speeds of the decision variables to prevent the vaccine supply chain from get into chaos. (3) The stability domain of the system shrinks as customers' sensitivity to cold chain transportation decreases, and by contrast it expends as customers' sensitivity to vaccine prices decreases. (4) Compared with the internal joint control of the distributor and retailer, the effect of external control, such as government intervention, will have a better control effect on the chaos of the system.