- About this Journal
- Abstracting and Indexing
- Aims and Scope
- Annual Issues
- Article Processing Charges
- Articles in Press
- Author Guidelines
- Bibliographic Information
- Citations to this Journal
- Contact Information
- Editorial Board
- Editorial Workflow
- Free eTOC Alerts
- Publication Ethics
- Reviewers Acknowledgment
- Submit a Manuscript
- Subscription Information
- Table of Contents
Journal of Applied Mathematics
Volume 2013 (2013), Article ID 927369, 11 pages
Modeling and Analysis of Bifurcation in a Delayed Worm Propagation Model
1College of Information Science and Engineering, Northeastern University, Shenyang 110819, China
2Key Laboratory of Medical Image Computing, Northeastern University, Ministry of Education, Shenyang 110819, China
Received 18 January 2013; Revised 4 August 2013; Accepted 26 August 2013
Academic Editor: Yannick De Decker
Copyright © 2013 Yu Yao et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
A delayed worm propagation model with birth and death rates is formulated. The stability of the positive equilibrium is studied. Through theoretical analysis, a critical value of Hopf bifurcation is derived. The worm propagation system is locally asymptotically stable when time delay is less than . However, Hopf bifurcation appears when time delay passes the threshold , which means that the worm propagation system is unstable and out of control. Consequently, time delay should be adjusted to be less than to ensure the stability of the system stable and better prediction of the scale and speed of Internet worm spreading. Finally, numerical and simulation experiments are presented to simulate the system, which fully support our analysis.
In recent years, Internet is undoubtedly one of the fastest increasing scientific technologies, which brings about convenience in people’s daily work and changes people’s life in variety of aspects. With rapid development of network applications and the increase of network complexity, security problems emerge progressively. Among them, the problem of Internet worms has become the focus with its wide infection range, fast spread speed, and tremendous destruction. Enlightened by the researches in epidemiology, plenty of models have been constructed to predict the spread of worms and some containment strategies have been taken into consideration. In addition, birth and death rates are widely applied in epidemiology because individuals in the ecological system may die during the spread of diseases. Meanwhile baby individuals are born everyday and join the ecological system [1–3]. In the computer science field, computers are like individuals in an ecological system. As a result of being infected by Internet worms or quarantined by intrusion detection systems (IDS), hosts will get unstable and unreliable which will result in system reinstallation by their owners. Besides, when many new computers are brought, most of them are preinstalled with operating systems without newest safety patches. Furthermore, old computers are discarded and recycled at the same time. These phenomena are quite similar to the death and birth in epidemiology. Thus, in order to imitate the real world, birth and death rates should be introduced to worm propagation model.
Quarantine strategies have been exploited and applied in the control of disease. The implementation of quarantine strategy in computer field relies on the IDS . The IDS include two major categories: misuse intrusion detection and anomaly intrusion detection system. The anomaly detection system is commonly used to detect malicious code such as computer virus and worms for its relatively better performance [5, 6]. Once a deviation from the normal behavior is detected, such behavior is recognized as an attack and appropriate response actions, such as quarantine or vaccination, are triggered.
Mechanism of time window was brought in IDS in order to balance the false negative rate and false positive rate . The introduction of time window is used to decide whether an alarm is a true or false positive based on the number of abnormal behaviors detected in a time window. It implies that the size of time window affects both the number of true positive and false positive rates. However, the import of the time window leads to time delay. Therefore, in order to accord with actual condition, time delay should be considered. In this paper, time delay is introduced in the worm propagation model along with birth and death rates, and its stability is analyzed. Moreover, it is indicated that overlarge time delay may result in the bifurcation which would do little help to eliminate the worms. Consequently, in order to guarantee the simplification and stability of the worm propagation system, time delay should be decreased appropriately by a decrease in the window size.
The rest of the paper is organized as follows. In the next section, related work on time delay and birth is death rates and introduced. Section 3 gives a brief introduction of the simple worm propagation model and quarantine strategy. Afterwards we present the delayed worm propagation model with birth and death rates and analyze the stability of the positive equilibrium. In Section 4, numerical and simulation experiments are presented to support our theoretical analysis. Finally, Section 5 draws the conclusions.
2. Related Work
Due to the high similarity between the spread of infectious biological viruses and computer worms, some scholars have used epidemic model to simulate and analyze the worm propagation [8–14]. For instance, Staniford first constructs the propagation of Internet worms by imitating epidemic propagation models, called simple epidemic model (SEM) model . susceptible-infected-removed (SIR) model plays a significant role in the research of worm propagation model . On the basis of SIR model, Zou et al. propose a worm propagation model with two factors on Code Red . Ren et al. give a novel computer virus propagation model and study its dynamic behaviors . Mishra and Pandey formulate an e-epidemic SIRS model for the fuzzy transmission of worms in computer network . L.-X. Yang and X. Yang investigate the propagation behavior of virus programs and propose that infected computers are connected to the Internet with positive probability . Mathematical analysis and simulation experiment of these models are conducted, which are helpful to predict the speed and scale of Internet worm propagation. In addition, to our knowledge, the use of quarantine strategies has produced a great effect on controlling disease. Enlightened by this, quarantine strategies are also widely used in worm containment [15–18]. Yao et al. construct a worm propagation model with time delay under quarantine, and its stability of the positive equilibrium is analyzed [15, 16]. Yao also proposes a pulse quarantine strategy to eliminate worms and obtains its stability condition . Wang et al. propose a novel epidemic model which combines both vaccinations and dynamic quarantine methods, referred to as SEIQV model .
Furthermore, some scholars have done some researches on time delay [19–21]. Dong et al. propose a computer virus model with time delay based on SEIR model . By proposing an SIRS model with stage structure and time delays, Zhang et al. perform some bifurcation analysis of this model . Zhang et al. also consider a delayed predator-prey epidemiological system with disease spreading in predator population . In addition, the direction of Hopf bifurcations and the stability of bifurcated periodic solutions are studied, which are our extension direction in the future research.
3. Worm Propagation Model
3.1. The Simple Worm Propagation Model
Realizing the similarities between Internet worms and biological viruses in propagation characteristics, classical epidemic models have been applied to the research of worm propagation models. Initially, we introduce a simple propagation model, the Kermack Mckendrick model (KM model) , as the basis of our research.
The KM model assumes that all Internet hosts are in one of three states: susceptible state (), infectious state (), and removed state (). And hosts can only maintain one state at any given moment. Infectious hosts are converted from susceptible hosts by a worm infection and can shift to removed state through killing worms by antivirus softwares and installing safety patches. Once patches are installed, the hosts can no longer be infected by worms. The state transition diagram of the KM model is given in Figure 1.
denotes the number of susceptible hosts at time , denotes the number of infectious hosts at time , and denotes the number of removed hosts at time . The total number of hosts in the network is and remains constant. is infection rate at which susceptible hosts are infected by infectious hosts, and is recovery rate at which infectious hosts get recovered. The KM model can be formulated by a set of differential equations from the state transition diagram as follows:
Although KM model adopts recovery feature and does generate some braking containment effect on the worm propagation, it only describes the initial stage of worm propagation and does not control the outbreaks of worms. More suppression strategies should be taken to further control the worm propagation.
3.2. The Worm Propagation Model with Quarantine Strategy
Quarantine strategy, which relies on the intrusion detection system, is an effective way to diminish the speed of worm propagation. On the basis of the KM model, quarantine strategy should be taken into consideration. Firstly, infectious hosts are detected by the systems and then get quarantined and patched. Moreover, considering that hosts could get patched whatever state the hosts stay, we add a new path from hosts in susceptible state to vaccinated state to accord with actual situation. The state transition diagram of the worm propagation model with quarantine strategy is given in Figure 2.
In this model, vaccinated state equals to removed state in KM model, and denotes the number of vaccinated hosts at time . Susceptible hosts can be infected by worms with an infection rate due to their vulnerabilities. After infection, the hosts become infectious hosts, which means the hosts can infect other susceptible hosts. The hosts can be directly patched into the vaccinated state before getting infected at rate . By applying misuse intrusion detection system for its relatively high accuracy, we add a new state called quarantine state, but only infectious hosts will be quarantined. denotes the number of quarantine hosts at time . Infectious hosts will be quarantined at rate which depends on the performance of intrusion detection systems and network devices. When infectious hosts are quarantined, the hosts will get rid of worms and get patched at rate . Meanwhile, infectious hosts can be manually patched into vaccinated state at recovery rate . It can be perceived that if a host is patched, it has gained immunity permanently.
The differetial equations of this model are given as
The total number of this model is set to , which means the sum of hosts in all four states is equal to one; that is . Then, the infection-free equilibrium of this model and its stability condition will be studied.
Due to the infection-free state which the system holds in the end, the number of infectious hosts is equal to 0. Obviously, the number of susceptible hosts and number of quarantined hosts are both equal to 0, because all of the hosts in the network will get patched at last, no matter which way the hosts take. Thus the infection-free equilibrium point of system (3) can be obtained. Since , the infection-free equilibrium of the model with quarantine strategy is .
Although all infectious hosts in quarantine model convert to vaccinated hosts, in other words, worms have been eliminated, it is hardly appropriate for real world situation. Actually, not all the hosts will convert to vaccinated hosts, and no-patch hosts are still existing in the network and suffering high risk of worm attack. In addition, considering that hosts are consumer electronic products, the recycling of old hosts happens frequently every day. In order to imitate the facts in the real world, birth and death rates must be taken into consideration.
Additionally, due to the time windows of intrusion detection system, which decreases the number of false positives, time delay should be considered to accord with actual situation. Therefore, in the next section, the new model is proposed.
3.3. The Delayed Worm Propagation Model with Birth and Death Rates
By adding time delay, along with birth and death rates, the delayed worm propagation model with birth and death rates is presented. Figure 3 shows the state transition diagram of the delayed worm propagation model with birth and death rates.
In this model, the total number of hosts denoted by is divided into five parts. Compared with quarantine model, a new state, delayed state, is added. The hosts do not have the ability to infect other susceptible hosts but are not quarantined yet. Therefore, a new state is needed to represent these hosts. denotes the hosts in delayed state at time . Furthermore, the hosts in susceptible, infected, quarantined, and removed states have the same rate to leave the network. If the hosts are quarantined by IDS, some of applications may not access to network because their ports are occupied by worms. At this moment, the users would like to reinstall OS and leave the network system. Thus, the birth and death rates are correlative to the number of infectious hosts and quarantined hosts in the network system. And the rate is
The hosts in delayed state are not likely to leave the network, accounting for the activeness of these hosts. The newborn hosts enter the system with the same rate , of which a portion is recovered by installing patches instantly at birth. This differs from the transition denoted by , because the hosts are vaccinated at the time of entering the network by installing patched operating system. The transition represents the rate by which those susceptible in the network get vaccinated by installing patches later. The delayed differential equations are written as
Other notations are identical to those in the previous model. To understand them more clearly, the notations in this model are shown in Table 1.
As mentioned in Table 1, the population size is set to , which is set to unity:
3.4. Stability of the Positive Equilibrium and Bifurcation Analysis
Theorem 1. The system has a unique positive equilibrium , when the condition is satisfied, where
Proof. From system (5), according to , if all the derivatives on the left of equal sign of the system are set to 0, which implies that the system becomes stationary, we can get Substituting the value of each variable in (8) for each of (6), then we can get: After solving this equation with respect to , a solution of is obtained Thus, if is satisfied, (10) has one unique positive root , and there is one unique positive equilibrium of system (5). The proof is completed.
Theorem 2. The positive equilibrium is locally asymptotically stable without time delay, if the following holds: ,where is satisfied.
Proof. If , then (14) reduces to
Because , equation (19) can be further reduced to
According to Routh-Hurwitz criterion, all the roots of (20) have negative real parts. Therefore, it can be deduced that the positive equilibrium is locally asymptotically stable without time delay. The proof is completed.
Obviously, () is a root of (14). After separating the real and imaginary parts, it can be written as
Lemma 3. Suppose that , is satisfied.(i)If one of followings holds: (a) , ; (b) , , and , then all roots of (14) have negative real parts when , and is a certain positive constant.(ii)If the conditions (a) and (b) are not satisfied, then all roots of (14) have negative real parts for all .
Proof. When , equation (14) can be reduced to
By the Routh-Hurwitz criterion, all roots of (20) have negative real parts if and only if , , and .
Considering (26), it is easy to see from the characters of cubic algebra equation that is a strictly monotonically increasing function if . If and or , but , then has no positive root. Thus, under these conditions, equation (14) has no purely imaginary roots for any . In addition, this also implies that the positive equilibrium of system (5) is absolutely stable. Therefore, the following theorem on the stability of positive equilibrium can be easily obtained.
Theorem 4. Assume that and are satisfied and and or , , and . Then the positive equilibrium of system (5) is absolutely stable; that is to say is asymptotically stable for any time delay .
It is assumed that the coefficients in satisfy the condition as follows: .
Then, according to lemma in , it is known that (26) has at least a positive root ; namely, the characteristic equation (14) has a pair of purely imaginary roots .
In view of the fact that (14) has a pair of purely imaginary roots , the corresponding is given by eliminating in (21):
Let be the root of (14) so that and are satisfied when .
Lemma 5. Suppose that . If , then is a pair of simple purely imaginary roots of (14). In addition, if the conditions in Lemma 3 are satisfied, then
This signifies that there exists at least one eigenvalue with positive real part for . Differentiating both sides of (14) with respect to τ, it can be written as Therefore where . Then it can be derived that
It follows the hypothesis that . Therefore, transversality condition holds and the conditions for Hopf bifurcation are satisfied when , according to Hopf bifurcation theorem  for functional differential equations. Then the following results can be obtained.
Theorem 6. Suppose that the conditions and are satisfied.(i)For system (5), the equilibrium is locally asymptotically stable when but unstable when .(ii)When condition is satisfied, system (5) will undergo a Hopf bifurcation at the positive equilibrium when (), where is defined by (28).
Theorem 2 implies to us that when birth and death rates and the mechanism of time window in intrusion detection system are considered, the length of time window is crucial for the stability of the propagation process. When time delay is less than the threshold , the system will stabilize at its infection equilibrium point, which is beneficial to further implement containment strategies to eliminate worms. However, when time delay is greater than the threshold, the system will be unstable and undergo a bifurcation. Although the propagation of worm presents a periodical phenomenon, in real world, complicated environment may make the propagation pass the critical state and reach to an unpredictable chaos.
4. Numerical and Simulation Experiments
4.1. Numerical Experiments
The parameters of the delayed model will be chosen properly, according to the stability of the positive equilibrium, bifurcation analysis, and the practical environment. 500,000 hosts are picked as the population size, and the worm’s average scan rate is per second. The worm infection rate can be calculated as , which means that average 0.4657 hosts of all the hosts can be scanned by one host. The infection ratio is . The rest of the parameters are , , , , , and . Initially, , which means that there are 50 infectious hosts while the rest of the hosts are susceptible at the beginning.
Figure 4 shows the curves of five kinds of hosts when . All of the five kinds of hosts get stable within 4 minutes, which illustrates that is asymptotically stable. It implies that the number of infectious hosts maintain a constant level and thus can be predicted. Further strategies can be developed and utilized to eliminate worms. But when time delay gets increased and then reach the threshold , will lose its stability and a bifurcation will occur. Figure 5 shows the susceptible, infected, and recovered hosts when . In this firgure, we can clearly see that the number of infectious hosts will outburst after a short period of peace and repeat again and again but not in the same period.
In order to see the influence of time delay, is set to a different value each time with other parameters remaining the same. Figure 6 shows the number of infected hosts in the same coordinate with time delays , , , and , respectively.
Figure 7 gives the four curves in four coordinates. In these two figures, the impact of time delay on infected hosts is evident. Initially, the four curves are overlapped which means that time delay has little effect in the initial stage of worm propagation. With the increase of time delay, the curve begins to oscillate. When time delay passes through the threshold , the infecting process gets unstable. Then simultaneously, it can be discovered that the amplitude and period of the number of infected hosts get increased.
Figure 8(a) shows the projection of the phase portrait of system (5) when in -space. The curve converges to a fixed point which suggests that the system is stable. Figure 8(b) shows the projection of the phase portrait of system (5) when in -space. The curve converges to a limit circle which implies that the system is unstable. Figures 8(c) and 8(d) show the phase portraits of susceptible and infected as and , respectively.
4.2. Simulation Experiments
The discrete-time simulation is an expanded version of Zou’s program simulating Code Red worm propagation and has been modified to run on a Linux server. The system in our simulation experiment consists of 500,000 hosts that can reach each other directly, which is consistent with the numerical experiments, and there is no topology issue in our simulation. At the beginning of simulation, 50 hosts are randomly chosen to be infectious and the others are all susceptible.
Identical with the results of numerical experiments, Figure 9 shows all five kinds of hosts when . We find that the model will reach a stable state after a short period of time, which suggests that the number of infetious hosts can be predicted and we can develop further strategy to eliminate worms.
When time delay increases but is less than the threshold derived from theoretical analysis, the number of infectious hosts and other hosts present a damped oscillation and finally reach an approximately stable state. Figure 10 shows the number of five kinds of hosts when . An interesting result is that the number of each host maintain a slight random vibration. However, the overall amplitude is only 0.273% to the population size, and the maximum amplitude is only 0.953% to the population size. The reason why this phenomenon occurs is due to the random events in simulation experiments. The implement of transition rates of the model is based on probability. That is to say, when the removal rate of infectious hosts is set to 0.03, it means that the probability of transition from infectious hosts to vaccinated state is 3%; it also means that infectious hosts get patched with probability of 3%. To investigate this phenomenon, phase portrait of susceptible and infected is depicted when as shown Figure 11. It shows that the curve is converged to a fixed point, which means that the system gets stable.
When time delay passes the threshold, a bifurcation appears. Figure 12 shows the number of susceptible, infectious and vaccinated hosts when . Figure 13 shows the difference between simulation and numerical results of susceptible and infectious hosts, respectively. In this figure, the periods of these two curves are well matched, which suggests that the results of numerical and simulation experiments are identical. But there is a difference of 9.6% of total population size in amplitudes. This mainly results from the precision of experiments. The number of hosts in numerical experiments can be either integers or decimals because it is the solution of differential equations. However, in simulation experiments, the number of hosts must be integers. Furthermore, when the number of infectious hosts are very little, this tiny difference will result in a big gap afterwards.
In order to accord with actual facts in the real world, time delay generated by time window in IDS is introduced to construct the worm propagation model. Dynamic birth and death rates are considered in this paper, accounting for the reinstallation of OS which users are more likely to do after when the hosts suffer worm’s destruction or quarantined to the Internet. In addition, combined with birth and death rates, time delay may lead to bifurcation and make the worm propagation system unstable.
In this paper, a delayed worm propagation model with birth and death rates is studied. Next, the stability of the positive equilibrium is analyzed. Through theoretical analysis and numerical and simulation experiments, the following conclusions can be derived.(i)The introduction of time window in IDS may lead to time delay in worm propagation model, which results in bifurcation. The critical time delay where the bifurcation appears is obtained.(ii)The worm propagation system is stable when time delay . In this situation, the worm propagation can be easily predicted and eliminated.(iii)A bifurcation emerges when time delay , which means the worm propagation is unstable and may be out of control.
Consequently, time delay should remain less than by decreasing the time window in IDS, which is helpful to predict the worm propagation and even eliminate the worm.
In this paper, we have only discussed the cases which satisfy conditions , , and . But for cases that satisfy , , and , the dynamical behavior of this model has not been studied. These issues will be studied and discussed in our future work.
The authors are very grateful to the anonymous referee for many valuable comments and suggestions which led to a significant improvement of the original paper. This paper is supported by the National Natural Science Foundation of China under Grant no. 60803132; the Natural Science Foundation of Liaoning Province of China under Grant no. 201202059; Program for Liaoning Excellent Talents in University under LR2013011; the Fundamental Research Funds of the Central Universities under Grant nos. N120504006 and N100704001; and MOE-Intel Special Fund of Information Technology (MOE-INTEL-2012-06).
- N. Yoshida and T. Hara, “Global stability of a delayed SIR epidemic model with density dependent birth and death rates,” Journal of Computational and Applied Mathematics, vol. 201, no. 2, pp. 339–347, 2007.
- M. Song, W. Ma, and Y. Takeuchi, “Permanence of a delayed SIR epidemic model with density dependent birth rate,” Journal of Computational and Applied Mathematics, vol. 201, no. 2, pp. 389–394, 2007.
- J. Liu and T. Zhang, “Bifurcation analysis of an SIS epidemic model with nonlinear birth rate,” Chaos, Solitons and Fractals, vol. 40, no. 3, pp. 1091–1099, 2009.
- V. Yegneswaran, P. Barford, and D. Plonka, “On the design and use of internet sinks for network abuse monitoring,” Computer Science, vol. 3224, pp. 146–165, 2004.
- D. Yeung and Y. Ding, “Host-based intrusion detection using dynamic and static behavioral models,” Pattern Recognition, vol. 36, no. 1, pp. 229–243, 2003.
- N. Stakhanova, S. Basu, and J. Wong, “On the symbiosis of specification-based and anomaly-based detection,” Computers and Security, vol. 29, no. 2, pp. 253–268, 2010.
- G. P. Spathoulas and S. K. Katsikas, “Reducing false positives in intrusion detection systems,” Computers and Security, vol. 29, no. 1, pp. 35–44, 2010.
- C. Lin, B. Liu, H. Hu, F. Xiao, and J. Zhang, “Detecting hidden Malware method based on “In-VM” model,” China Communications, vol. 8, no. 4, pp. 99–108, 2011.
- S. Staniford, V. Paxson, and N. Weaver, “How to own the internet in your spare time,” in Proceedings of the 11th USENIX Security Symposium, pp. 149–167, San Francisco, Calif, USA, 2002.
- S. Qing and W. Wen, “A survey and trends on Internet worms,” Computers and Security, vol. 24, no. 4, pp. 334–346, 2005.
- C. C. Zou, W. Gong, and D. Towsley, “Worm propagation modeling and analysis under dynamic quarantine defense,” in Proceedings of the 2003 ACM Workshop on Rapid Malcode (WORM '03), pp. 51–60, Washington, DC, USA, October 2003.
- J. Ren, X. Yang, Q. Zhu, L. Yang, and C. Zhang, “A novel computer virus model and its dynamics,” Nonlinear Analysis: Real World Applications, vol. 13, no. 1, pp. 376–384, 2012.
- B. K. Mishra and S. K. Pandey, “Fuzzy epidemic model for the transmission of worms in computer network,” Nonlinear Analysis: Real World Applications, vol. 11, no. 5, pp. 4335–4341, 2010.
- L.-X. Yang and X. Yang, “Propagation behavior of virus codes in the situation that infected computers are connected to the internet with positive probability,” Discrete Dynamics in Nature and Society, vol. 2012, Article ID 693695, 13 pages, 2012.
- Y. Yao, X. Xie, H. Guo, G. Yu, F. Gao, and X. Tong, “Hopf bifurcation in an Internet worm propagation model with time delay in quarantine,” Mathematical and Computer Modelling, 2011.
- Y. Yao, W. Xiang, A. Qu, G. Yu, and F. Gao, “Hopf bifurcation in an SEIDQV worm propagation model with quarantine strategy,” Discrete Dynamics in Nature and Society, vol. 2012, Article ID 304868, 18 pages, 2012.
- Y. Yao, L. Guo, H. Guo, G. Yu, F. Gao, and X. Tong, “Pulse quarantine strategy of internet worm propagation: modeling and analysis,” Computers and Electrical Engineering, 2011.
- F. Wang, Y. Zhang, C. Wang, J. Ma, and S. Moon, “Stability analysis of a SEIQV epidemic model for rapid spreading worms,” Computers and Security, vol. 29, no. 4, pp. 410–418, 2010.
- T. Dong, X. Liao, and H. Li, “Stability and Hopf bifurcation in a computer virus model with multistate antivirus,” Abstract and Applied Analysis, vol. 2012, Article ID 841987, 16 pages, 2012.
- T. Zhang, J. Liu, and Z. Teng, “Stability of Hopf bifurcation of a delayed SIRS epidemic model with stage structure,” Nonlinear Analysis: Real World Applications, vol. 11, no. 1, pp. 293–306, 2010.
- J. Zhang, W. Li, and X. Yan, “Hopf bifurcation and stability of periodic solutions in a delayed eco-epidemiological system,” Applied Mathematics and Computation, vol. 198, no. 2, pp. 865–876, 2008.
- W. Kermack and A. McKendrick, “A contribution to the mathematical theory of epidemics,” Proceedings of the Royal Society of London A, vol. 115, no. 772, pp. 700–721, 1927.
- S. Wang, Q. Liu, X. Yu, and Y. Ma, “Bifurcation analysis of a model for network worm propagation with time delay,” Mathematical and Computer Modelling, vol. 52, no. 3-4, pp. 435–447, 2010.
- B. Hassard, D. Kazarino, and Y. Wan, Theory and Application of Hopf Bifurcation, Cambridge University Press, Cambridge, UK, 1981.