Modelling the transmission of pathogens by considering environmental- and direct transmission mechanisms

Hele tekst

(1)Modelling the transmission of pathogens by considering environmental- and direct transmission mechanisms. A. Scholtz 20109660. Dissertation submitted in partial fulfilment of the requirements for the degree Magister Scientiae in Applied Mathematics at the Potchefstroom Campus of the North-West University.. Supervisor: Prof. I.M. Schoeman Co-supervisor: Prof. J. Spoelstra October 2015.

(2) Contents Abstract. vi. Acknowledgement. vii. Nomenclature. ix. Introduction and preliminaries. 1. 1 Mathematical modelling. 9. 1.1. Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 1.2. Modelling methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10. 1.3. Model fitting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12. 1.4. Dimensional analysis and scaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13. 1.5. Population growth models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19. 2 Mathematical techniques 2.1. 9. 26. Systems of ordinary differential equations . . . . . . . . . . . . . . . . . . . . . . . . 26 2.1.1. Adaptive Runge-Kutta-Fehlberg method (ode45) . . . . . . . . . . . . . . . . 26. 2.1.2. Stability of systems of differential equations . . . . . . . . . . . . . . . . . . . 28. 2.2. Parameter Sensitivity Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31. 2.3. Mathematical optimisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32. 3 Modelling of pathogen transmission. 35. 3.1. Disease transmission . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36. 3.2. Compartmental models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39. 3.3. 3.2.1. SI Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39. 3.2.2. SIR Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41. 3.2.3. Other Compartmental Models . . . . . . . . . . . . . . . . . . . . . . . . . . . 43. Basic Reproductive Number . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46. 4 Model analyses 4.1. 50. Li model description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50. i.

(3) 4.2. 4.3. 4.1.1. Non-dimensionalisation and scaling . . . . . . . . . . . . . . . . . . . . . . . . 51. 4.1.2. Stability analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54. Breban model description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.2.1. Non-dimensionalisation and scaling . . . . . . . . . . . . . . . . . . . . . . . . 57. 4.2.2. Stability analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59. Proposed model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 4.3.1. Non-dimensionalisation and scaling . . . . . . . . . . . . . . . . . . . . . . . . 63. 4.3.2. Stability analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64. 5 Empirical case study. 68. 5.1. Data acquisition and methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68. 5.2. Simplified Kermack & McKendrick SIR model . . . . . . . . . . . . . . . . . . . . . . 69. 5.3. 5.2.1. Non-scaled and dimensional model . . . . . . . . . . . . . . . . . . . . . . . . 69. 5.2.2. Scaled and dimensionless model . . . . . . . . . . . . . . . . . . . . . . . . . . 71. Proposed Model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73. 5.3.1. Non-scaled and dimensional model . . . . . . . . . . . . . . . . . . . . . . . . 73. 5.3.2. Scaled and dimensionless model . . . . . . . . . . . . . . . . . . . . . . . . . . 75. 5.4. Parameter sensitivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77. 5.5. Summary and discussion of results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77. 6 Discussion and conclusion. 83. 6.1. Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83. 6.2. Conclusion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85. A Data. 87. B MATLAB code. 89. B.1 Simplified Kermack & McKendrick SIR-model . . . . . . . . . . . . . . . . . . . . . . 90 B.2 Simplified Kermack & McKendrick SIR-model (Scaled and dimensionless) . . . . . . 91 B.3 Proposed model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 B.4 Proposed model (Scaled and dimensionless) . . . . . . . . . . . . . . . . . . . . . . . 93. ii.

(4) List of Figures Figure 1:. Image illustrating the endocytosis and replication of influenza in a host cell.. 7. Figure 1.1: Graphs illustrating the similarities between a non-scaled and scaled system.. 19. Figure 1.2: Illustration of the Malthusian population growth model. . . . . . . . . . . . 21 Figure 1.3: Illustration of logistic population growth.. . . . . . . . . . . . . . . . . . . . 23. Figure 3.1: Illustration of compartments in an SI model. . . . . . . . . . . . . . . . . . . 39 Figure 3.2: Graph illustrating dynamics of the SI model. . . . . . . . . . . . . . . . . . 42 Figure 3.3: Illustration of compartments in an SIR model. . . . . . . . . . . . . . . . . . 42 Figure 3.4: Graph illustrating the dynamics of the SIR model. . . . . . . . . . . . . . . 43 Figure 3.5: Illustration of compartments in an SEIR model. . . . . . . . . . . . . . . . . 44 Figure 3.6: Illustration of compartments in an SITR model. . . . . . . . . . . . . . . . . 45 Figure 3.7: Illustration of compartments in an SLIAR model. . . . . . . . . . . . . . . . 45 Figure 4.1: Graph comparing the Li model to its dimensionless form. . . . . . . . . . . 54 Figure 4.2: Graph comparing the proposed model to its dimensionless form. . . . . . . . 65 Figure 5.1: SIR model fit to incidence data. . . . . . . . . . . . . . . . . . . . . . . . . . 69 Figure 5.2: SIR model showing dynamics of all compartments. . . . . . . . . . . . . . . 70 Figure 5.3: SIR model fit to scaled incidence (density) data.. . . . . . . . . . . . . . . . 72. Figure 5.4: SIR model (scaled) showing dynamics of all the compartments. . . . . . . . 72 Figure 5.5: Proposed model fit to incidence data. . . . . . . . . . . . . . . . . . . . . . . 74 Figure 5.6: Proposed model showing dynamics of population compartments. . . . . . . 74 Figure 5.7: Proposed model fit to scaled incidence (incidence density) data. . . . . . . . 76 Figure 5.8: Proposed model (scaled) showing dynamics of population compartments. . . 76 Figure 5.9: SIR model illustrating dynamics by varying the parameter β. . . . . . . . . 78 Figure 5.10: Proposed model illustrating dynamics by varying the parameter βE . . . . . 79 Figure 5.11: Proposed model illustrating dynamics by varying the parameter βD . . . . . 80 Figure 5.12: Proposed model illustrating dynamics by varying the parameter η. . . . . . 81 Figure B.1: Script used to extract data from file. . . . . . . . . . . . . . . . . . . . . . . 90 Figure B.2: Function script used for the Kermack and McKendrick (1932) ODE. . . . . 90 Figure B.3: Script for implementing estimated parameters for SIR-model. iii. . . . . . . . . 90.

(5) Figure B.4: Script used to extract and scale data from file.. . . . . . . . . . . . . . . . . 91. Figure B.5: Function script used for the scaled and dimensionless Kermack and McKendrick (1932) ODE. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 Figure B.6: Script for implementing estimated parameters.. . . . . . . . . . . . . . . . . 91. Figure B.7: Function script used for proposed model ODE. . . . . . . . . . . . . . . . . 92 Figure B.8: Script for implementing estimated parameters.. . . . . . . . . . . . . . . . . 92. Figure B.9: Function script used for proposed model ODE. . . . . . . . . . . . . . . . . 93 Figure B.10: Script for implementing parameters estimated by method of the optimisation algorithm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 Figure B.11: Implementing Levenberg-Marquardt via lsqnonlin. . . . . . . . . . . . . . 94. iv.

(6) List of Tables Table 1.1:. Dimensions of common physical quantities.. Table 2.1:. Coefficients used in the implementation of the Runge-Kutta-Fehlberg method. 27. Table 4.1:. Dimensions of parameters in the Li model. . . . . . . . . . . . . . . . . . . . 52. Table 4.2:. Dimensions of parameters in the Breban model. . . . . . . . . . . . . . . . . 58. Table 4.3:. Dimensions of parameters in the proposed model. . . . . . . . . . . . . . . . 62. Table 5.1:. Estimated parameters for non-scaled and dimensional SIR model. . . . . . . 70. Table 5.2:. Estimated parameters for scaled and dimensionless SIR model. . . . . . . . 71. Table 5.3:. Estimated parameters for non-scaled and dimensional proposed model. . . . 73. Table 5.4:. Estimated parameters for scaled and dimensionless proposed model. . . . . 75. Table 5.5:. Estimated parameters for scaled and dimensionless models. . . . . . . . . . 82. Table 5.6:. Varied parameter values for SIR- and proposed models.. Table A.1:. Influenza A data used in application model fitting. . . . . . . . . . . . . . . 88. v. . . . . . . . . . . . . . . . . . . 14. . . . . . . . . . . . 82.

(7) Abstract. Environmental transmission of pathogens is a concern in society. With the rise of emerging diseases such as influenza type A, Ebola and extensively drug-resistant tuberculosis, there is a demand for understanding the transmission of these diseases in order to devise control strategies. A thorough literature study is done on the methodology used in epidemiological modelling. Models incorporating environmental transmission mechanisms from literature are analysed and a model is proposed which incorporates both direct- and environmental transmission mechanisms. An empirical case study is done on influenza type A using acquired data. It is shown that the mechanism of environmental transmission is viable in human-to-human transmission. Keywords: Mathematical modelling, Mathematical epidemiology, Environmental transmission, Direct transmission, Basic reproduction number. vi.

(8) Acknowledgements. First and foremost, thank you to the South African Centre for Epidemiological Modelling and Analysis, the National Research Foundation, the North-West University and the School for Computer, Statistical and Mathematical Sciences for the funding I received. To my supervisors, Prof. I.M. Schoeman, Prof. J. Spoelstra and Mr. L.M. Viljoen for your constant support and motivation which made this study possible. Secondly, thank you to my friends and family. The people who stand by me and support me unconditionally. Every day. A special thanks to my fiancé Dawie van Nieuwenhuizen, and best friend Cobi Hayes, for your infinite, unconditional support, friendship and love. My lovely and loving mother Marie Scholtz, without whom my academic career would not have been what it is today. Also, my father André Scholtz, my brother Jacques Scholtz and my sister Michelle MacDonald. Lastly, thank you to all my friends, those who are new, those who are old, and those who have come and gone.. vii.

(9) “It is not for him to pride himself who loveth his own country, but rather for him who loveth the whole world. The earth is but one country, and mankind its citizens.” - Baha’u’llah.

(10) Nomenclature β-Coronavirus:. A class of coronavirus from the subfamily Coronavirinae from the family Coronaviridae.. Basic reproductive. The expected number of secondary infections that may arise due to. number:. the introduction of one typically infected individual in a totally susceptible population. Abbreviated: R0. Emerging infectious. A disease that is new to a host population or changing in some manner. disease:. within the host population.. Endemic disease:. A disease that is always present in a population.. Epidemic disease:. A disease that a large fraction of a population may acquire over a short period of time.. Fomite:. An inanimate object which may be contaminated by micro-organisms and facilitate their transmission.. Haemagglutinin:. A glycoprotein present as an envelope spike on influenza viruses. Used in identification of different influenza viruses. Abbreviated: HA, or H when describing influenza virus subtypes.. Human. A lentivirus of the family Retroviridae which causes Acquired Immune. immunodeficiency. Deficiency Syndrome (AIDS).. virus:. Abbreviated: HIV. Incidence:. The fraction of a population that contracts a certain disease during a given period of time.. MERS-coronavirus:. Middle East respiratory syndrome caused by a variant β-coronavirus. Abbreviated: MERS-CoV. Neuraminidase:. A glycoside hydrolase enzyme present as an envelope spike on influenza viruses. Used in identification of different influenza viruses. Abbreviated: NA, or N when describing influenza virus subtypes.. Pandemic disease:. A worldwide epidemic.. Pathogen:. A specific causative agent of disease.. Prevalence:. The fraction of the population that currently, or at a certain point in time, has a given disease.. ix.

(11) Ribonucleic acid:. A. polynucleotide. composed. of. ribonucleotides. joined. by. phosphodiethers. RNA-virus:. A virus of which the genetic material consists solely of ribonucleic acid.. Sporadic disease:. A disease that occurs only occasionally in a population.. x.

(12) Introduction and preliminaries Describing worldly phenomena from a mathematical point of view is fascinating; numerous procedures exist for finding approximate solutions to real-world problems. The usual approach to be followed when solving real-world problems is to initially analyse a problem based on known methods and existing research and then to develop the simplest, most elegant model which describes the physical phenomena accurately. Einstein made this clear, although he stated it broadly, not limiting this philosophy to physics or mathematics: “It can scarcely be denied that the supreme goal of all theory is to make the irreducible basic elements as simple and as few as possible without having to surrender the adequate representation of a single datum of experience.” Einstein (1934) In modelling physical systems, we attempt to develop elegant, simple models which are complex enough to describe physical systems or phenomena as accurately as possible. In particular, this modelling paradigm is also followed in mathematical modelling. Mathematical modelling attempts to answer problems of practical concern (Brauer, 2009). The modelling of biological phenomena is no exception. The modelling of infectious disease has become an area of interest in the mathematical sciences as much as it is of importance to the biological sciences. Mathematical application in infectious disease modelling was first introduced by Daniel Bernoulli in 1760 when he attempted to calculate the mortality of smallpox, and what the effects of inoculation would be on a population (Brauer and Castillo-Chavez, 2012:345). This opened a window of opportunity where scientists would be able to use mathematical tools to describe disease dynamics and develop control strategies in the case of threatening epidemics (Brauer and Castillo-Chavez, 2012:345). Important contributions to the modelling of infectious disease include, among many others, Ross’ malaria model (1911), Kermack and McKendrick’s compartmental model (1927), and the definition of the basic reproductive ratio by Diekmann et al. (1990). The first model of infectious disease modelling was provided by Kermack and McKendrick (1927). This compartmental SIR model (where SIR abbreviates Susceptible, Infected, and Removed) was an innovative idea, which created new opportunities and tools for scientists to model the spread of contagious diseases. This would be the novel approach in infectious disease modelling over the next few decades. Compartmental models have simplified the modelling of the propagation of infectious diseases within a population. These models are simple enough to construct, yet the question is whether such models deliver adequate results concerning their description of disease transmission between hosts, as well as their predictive power. The simplified deterministic model proposed by Kermack and McKendrick describes the between-host dynamics of an infectious disease in a closed population. In its simplicity and elegance, this model has laid the foundation for compartmental modelling of infectious disease dynamics.. 1.

(13) Introduction and preliminaries The Italian physician Girolamo Fracastoro proposed that infectious diseases are transmitted by either direct or indirect contact, and was the first to use the word “fomes”1 (Brock, 1961). Certain pathogens persist environmentally (Li et al., 2009) on fomites, in water and aerosols (Breban, 2013). Aerosols are defined as a solid or liquid suspension in a gas, of which the stability may vary between a few seconds to more than a year (Hinds, 2012:3-4). An assumption that suggests the inclusion of an environmental transmission mechanism into the SIR model owes to the fact that certain pathogens are environmentally persistent, and therefore environmental transmission of the pathogen is a contributing factor to overall transmission dynamics. This has already been modelled by Codeço (2001), Breban et al. (2009), Li et al. (2009), Breban et al. (2010) and Kisdi and Boldin (2013). More recently, models reflecting SIR-dynamics have been implemented to a great extent to describe the prevalence of modern-day infections such as human immunodeficiency virus (HIV), human coronavirus (HCV), influenza, etc. According to Li et al. (2009) and Breban (2013), the classical approach to modelling the propagation of infectious disease in a population has mostly considered direct transmission mechanisms. These models imply that each susceptible individual must come into physical contact (defined as direct contact) with an infected individual before the susceptible individual has been exposed and is possibly infected. In this study, focus will be placed on the effect of pathogens transmitted by means of an environmental transmission mechanism, given that the pathogen is persistent in the environment for a reasonable measure of time. Typical of mathematical modelling, there will always be a quest to develop a way to describe a system more efficiently, in order to infer important information from this system (Velten, 2009; Chubb and Jacobsen, 2010). A thorough literature review of epidemiological models is presented (see Chap. 3). Modelling of infectious diseases is approached from a mathematical background, and thorough mathematical analyses and verification on models from literature are done, including those proposed by Kermack and McKendrick (1927), Li et al. (2009) and Breban (2013) (see Chap. 4). These analyses include dimensional analysis and scaling (see Chap. 1), parameter sensitivity analysis, stability analysis and numerical analysis (see Chap. 2). In Chapter 4 (Section 4.3), we propose a direct and environmental transmission model for the modelling of influenza type A, which is analysed and verified based on the models by Kermack and McKendrick 1927, Li et al. (2009) and Breban (2013). Data have been obtained on the prevalence of influenza type A in citizens from the United States. In Chapter 5, the data set is fit to a model containing only a direct transmission mechanism and then to a model with both direct and environmental transmission mechanisms. The model is also shown to accurately describe the disease dynamics with an included environmental transmission mechanism. R0 is calculated for influenza A based on both the original SIR model and the proposed model. Furthermore, the calculated value of R0 is compared to R0 from literature. Also, the issue of whether environmental transmission is a viable modelling tool to be included in deterministic models is discussed. The question of whether the work done in this study is viable is answered by an empirical case study, 1. From Latin, touchwood, tinder. pl. fomites. Inanimate objects which may be contaminated by infectious agents. (Merriam-Webster, 2006). 2.

(14) Introduction and preliminaries where the new model is compared to work done by Li et al. (2009) and Breban (2013). The study concludes with a discussion of the results, points of concern, and key elements which should be addressed through further investigation.. Problem statement Does inclusion of environmental transmission, by means of adapting an existing SIR model, deliver better results in modelling disease transmission dynamics in a population, and what information can be inferred from estimated parameters?. Aim • We discuss and analyse the consequence of environmental- and direct transmission mechanisms on an adapted SIR model. • Investigate whether inclusion of an environmental transmission mechanism into an SIR model is essential in obtaining accurate model descriptions for certain classes of infectious disease. • An empirical case study is used to compare the classic SIR model by Kermack and McKendrick (1932) and a model incorporating direct- and environmental transmission proposed in Chapter 4.. Outline The outline of this study is summarised as follows. Introduction An introduction to the topic of the study, and introduces the reader to previous literature regarding between-host dynamics; mathematical epidemiology, specifically referring to studies on environmental transmission of pathogens; the history of infectious diseases in society; and the biology of influenza. The research problem is identified and the aim of the study is defined, followed by a detailed outline of the study. Chapter 1: Mathematical modelling An introduction and discussion on the preliminaries and basic principles, as well as the process, of mathematical modelling. This chapter introduces models of population dynamics including population growth and population interaction in terms of predator-prey models. Chapter 2: Mathematical techniques An overview of mathematical techniques used in analyses of systems of differential equations used in mathematical modelling, solving of systems of differential equations, stability of a system of differential equations, sensitivity analysis, and mathematical optimisation. 3.

(15) Introduction and preliminaries Chapter 3: Modelling of pathogen transmission A discussion on the modelling of pathogen transmission, a review of selected models in disease dynamics from literature, and an overview of the basic reproduction number. Chapter 4: Model analyses Thorough analyses of the models proposed by Li et al. (2009) and Breban (2013) are done and an adapted SIR model with an incorporated environmental transmission mechanism is proposed and analysed. Chapter 5: Empirical case study Influenza type A data were acquired from the Centre for Disease Control and Prevention in the United States of America. These data are used to fit the classic SIR model and the model proposed in Chapter 4; in both cases for unscaled and dimensional, and scaled and dimensionless versions of these models. The parameters in the models are estimated by means of the Levenberg-Marquardt optimisation algorithm. Furthermore, the value of R0 is calculated from the estimated parameters. Chapter 6: Discussion and conclusion Discussion of the results obtained in Chapter 5 and conclusion of this study.. 4.

(16) Introduction and preliminaries. Literature review of contagious diseases A short history of contagious diseases in society Contagious diseases are diseases which are transmitted between individuals through direct contact with an infected individual, the discharges of an infected individual or by indirect means (MerriamWebster, 2006). The adjectives “contagious” and “communicable” will be used interchangeably, as both refer to the same set of diseases in this study, i.e. diseases which may be transmitted from one individual to another. The following discussion is based on the text by Brauer and Castillo-Chavez (2012:345–349). Diseases caused by pathogens such as influenza, tuberculosis, measles, etc. have threatened society for years. These diseases are a well-known fact of life in the western world, giving rise to numerous studies being conducted in order to devise methods of control and intervention. Even though many of the known pathogens are easily treatable, they still threaten less developed countries. The effects and prevalence in these countries are mostly unknown, but are of great importance. Communicable disease outbreaks and epidemics with high mortality in less developed countries have devastating effects on the mean life expectancy of individuals in a population, and this may ultimately have dire consequences for the economies of these countries. Emerging diseases such as the MERS-coronavirus (MERS-CoV), a variant β-CoV, which recently emerged, can wreak havoc in the world (Lu and Liu, 2012). It is therefore important that disease transmission is modelled effectively in order to make inferences from the proposed models and develop control procedures. The idea of microscopic life being the cause of disease goes back to as early as Aristotle (384 – 322 BC), but it only developed as a theory in the 16th century. The existence of micro-organisms was demonstrated by Van Leeuwenhoek using a predecessor of the microscope. In 1840 Jacob Henle introduced the first expression of the germ theory of disease. This theory was then further developed by Robert Koch, Joseph Lister and Louis Pasteur from the late 19th to the early 20th century. The mechanisms of transmission of many diseases are generally well known. Diseases caused by viruses such as CoV, influenza, measles, etc. usually confer immunity to reinfection. Diseases transmitted by bacteria such as tuberculosis, bacterial meningitis, pneumonia, etc. do not confer immunity in hosts, who are then open to reinfection. Diseases like malaria are parasitic infections which usually spread by means of a vector, such as an insect, from person to person. The sudden outbreak of a disease which affects a substantial proportion of a population and then disappears, is referred to as an epidemic. Epidemics, however, usually leave a certain proportion of the population untouched. It does happen that recurring outbreaks occur in the same population, but with lesser effect on the total population due to acquired immunity within the host population. This attribute of recurrent epidemics draws an important connection between epidemics and pathogen evolution.. 5.

(17) Introduction and preliminaries The western historian W.H. McNeill 1998 argues that communicable diseases have had an important influence on society in history. In the 18th century, there was a sudden sharp increase in population throughout the world. Many factors may have contributed to this increase, but the single factor that has been isolated was the decrease in mortality due to recurring epidemic infections. Advances in medicine and increased travel contributed to this decline in mortality, because of the increased circulation and co-circulation of diseases. Historical accounts of epidemic diseases can even be found in biblical references as far back as in the time of Moses. The fall of several great empires were due to epidemic outbreaks. The fall of the Roman Empire during the Antonine plagues, which may have been a measles epidemic, contributed to the collapse of the Roman economy, and eventually caused the empire to disintegrate into disorganisation. Communicable disease was the precursor of the fall of the Han empire in China in the 3rd century AD. The defeat of millions of Aztecs by Hernán Cotrés and his company of only 600 men was aided by smallpox, which travelled with them all the way from Europe. The disease spread as far as Peru, devastating the Incas whom Cortés later defeated as well. Other illnesses, such as measles and diphtheria, also made their way from Europe to North America accompanying the early European explorers. In 1346, the first wave of the Black Death (or bubonic plague) made its way from Asia to Europe and was followed by several other recurrences. It is estimated that this epidemic caused mortalities of as much as a third of the entire European population after only four years since its introduction, and it swept through Europe for the next 300 years. This recurring epidemic had an immense impact on the political and socio-economic development in Europe during the medieval era.. Orthomyxoviridae: Influenza Type A The following discussion on influenza was adopted from Willey et al. (2014:856–859): Avian influenza, commonly referred to as ‘bird flu’, is an infectious disease of birds. It does happen, however, that some strains can cause serious infections in humans. Influenza was first described by Hippocrates in 412 BC. The first documented influenza epidemic in recorded history occurred in 1580; approximately 31 possible epidemics have occurred since, and four of these occurred in the 21st century. A new influenza strain was documented in 2011, an H3N2 variant which has its origin in pigs. Influenza is a negative-strand RNA-virus belonging to the family Orthomyxoviridae. It is a respiratory disease transmitted mainly through the ingestion or inhalation of respiratory secretions. The virus incubates in the body for about 1 to 2 days, adhering to epithelial cells facilitated by haemagglutinin proteins on the virus’s envelope. Another enzyme found on the virus envelope, neuraminidase, may assist in hydrolysing the mucus layer of the respiratory tract, enabling adherence to the epithelial layer. The attached virus then enters the cell by means of receptor-mediated endocytosis, whereby the virus is enclosed in an endosome and enters the epithelial cell. When the endosomal pH starts to decrease, the haemagglutinin molecules on the virus envelope undergo a dynamic conformational change in structure: the hydrophobic tail-ends of the haemagglutinin turn outward and extend toward the membrane of the endosome. Contact between the haemagglutinin 6.

(18) Introduction and preliminaries. Figure 1. Image illustrating the endocytosis and replication of influenza in a host cell. Source: Shi et al. (2014). and the endosome membrane results in fusion, and this in turn results in the release of the RNAnucleocapsid into the cellular cytoplasm. The process of replication is then initiated. The viral RNA enters the epithelial cell nucleus where it is then replicated and transcribed. Newly formed ribonucleo-proteins are then transported into the cellular cytoplasm, and assembly into virions occurs at the cell membrane. Neuraminidase then facilitates the release of virions from the infected cells. The host cell ruptures, and mature virions take the plasma membrane of the host cell as their envelope. The classification of influenza into different subtypes is attributed to the virus’s surface glycoproteins located on the envelope, namely haemagglutinin (HA) and neuraminidase (NA). At this time there are 18 HA and 11 NA antigenic forms known (Tong et al., 2013). Combinations of these two glycoproteins are used to classify subtypes. The efficiency of these proteins defines the virulence of influenza type A – haemagglutinin facilitates the attachment and entrance into the host cell, and neuraminidase facilitates the process of viral shedding. The host specificity of influenza type A is associated with differences in receptor binding, therefore the influenza virus prefers binding to the cell membrane of human tracheal epithelial cells. Humans are not infected by avian influenza viruses as easily because of the virus’s preference to adhere to the intestinal epithelial cells of birds, restricting transmission (for the most part) to birds. The genome of the influenza virus, however, is very plastic. This genomic plasticity can contribute to changes in the viral proteins which are synthesised, modifying the virus’s antigenicity. Small changes in the genome due to an accumulation of mutations are referred to as antigenic drift. When two different strains of virus infect the same host cell, however, genetic reassortment of the. 7.

(19) Introduction and preliminaries genomes (and incorporation into a single new nucleocapsid) results in a large antigenic change and possible a new strain of virus. This process is called antigenic shift. There is a greater chance at any given time that antigenic shift will occur, with the potential to lead to major epidemics and even pandemics – as was the case with the 2009 H1N1 pandemic. Antigenic variations occur almost on a yearly basis in influenza A. The question now arises of how people contract influenza A infections. In this respect, animal reservoirs play a crucial role in the epidemiology of influenza A. Influenza A is common among birds, and although bird-to-human transmission rarely happens, bird-to-pig transmission happens on a regular basis. Pig-to-human transmission and vice versa happen easily. Thus, the recombinant reassortment of bird and human strains occur in pigs, leading to new strains as a result of major antigenic shift.. 8.

(20) Chapter 1. Mathematical modelling This chapter is an introductory discussion on the ideas and concepts that form an important part of the mathematical modelling process. In Section 1.1, theoretical concepts used to distinguish between variables and parameters are defined. Thereafter, in Section 1.2, a general framework for the methodology of mathematical modelling is discussed, followed by the procedure of fitting data to a model in Section 1.3. Dimensional analysis and scaling of systems (Section 1.4), and how this assists in the verification of terms in a model are discussed. Lastly, the concepts of population dynamics are discussed in Section 1.5. A mathematical model is a tool used by a modeller to describe a physical system and then infer information from the model in order to make decisions that will influence the physical system (Velten, 2009). This is described by Minsky (1965) simply as follows: “To an observer B, an object A∗ is a model of an object A to the extent that B can use A∗ to answer questions that interest him about A.” This definition of a model contains a valuable philosophy: by modelling a phenomenon, the modeller solves real-world problems and make decisions based on inferences made from the model Velten (2009). Moreover, according to Velten (2009), the best model to describe real-world phenomena is the simplest model that still serves its purpose (much like Einstein’s 1934 description), that is, the model should be complex enough to describe the system for the relevant question.. 1.1. Preliminaries. Mattheij and Molenaar (1991) explain that when a mathematical modeller encounters a real-world problem that needs to be solved, the first task of the modeller is to translate the subject matter into a mathematically tractable form. Mattheij and Molenaar (1991) define a mathematical model as “...the description of an experimentally verifiable phenomenon by means of the mathematical language”. Definition 1.1.1: (Mathematical model) Velten (2009) A mathematical model is a triplet (S, Q, M ) where S is a system, Q is a question relating to S and M is a set of mathematical statements which can be used to answer Q. . Blank text. Variables The first class of quantities contained in a mathematical model is variables, and within this class one must make a clear distinction between dependent and independent variables Mattheij and Molenaar. 9.

(21) 1. Mathematical modelling (1991). The dependent variables are differentiable with respect to the independent variables, and therefore the change in a system may be described by means of differential equations. In this study the only independent variable is time, and therefore attention is focussed on ordinary differential equations as opposed to partial differential equations, where more than one independent variable exist (Mattheij and Molenaar, 1991; Velten, 2009; Giordano et al., 2003; Kreyszig, 1999). Definition 1.1.2: (State variables) Velten (2009) Let (S, Q, M ) be a mathematical model. Mathematical quantities s1 , s2 , ..., sn which describe the state of the system S in terms of M and which are required to answer Q are called state variables of (S, Q, M ). . Blank text. Parameters The second class of quantities used in mathematical models is parameters. It is of importance to distinguish between two types of parameters, i.e. constant parameters and adjustable parameters. Constant parameters are those which are fixed according to universal laws, such as gravitational acceleration and radioactive decay time. Adjustable parameters are those that may be controlled by the experimentalist Mattheij and Molenaar (1991); Velten (2009). It is said that a mathematical model is solved when the dependent variables, in principle, are known as functions of the parameters and independent variables describing the system. Definition 1.1.3: (Reduced system and system parameters) (Velten, 2009) Let s1 , s2 , ..., sn be the state variables of a mathematical model (S, Q, M ). Let β1 , β2 , ...βm be mathematical quantities that describe properties of the system S in terms of M and which are needed to compute the state variables. Then Sr = {β1 , β2 , ..., βm } is known as the reduced system and β1 , β2 , ...βm are the system parameters of (S, Q, M ). . Blank text. 1.2. Modelling methodology. We now summarise the contents of the majority of this section into a framework for mathematical modelling. Mattheij and Molenaar (1991) give a general outline for the process of, or procedure to follow in, modelling: 1. Orientation: In this first stage the modeller gets acquainted with the system to be modelled through observations and by gathering information from experts and literature. 2. Formulation of relations between variables and parameters: The next step is formulating relations between variables and parameters that describe the system. Certain classes of systems have established relations which may be applied from the literature. There is, however, an important condition that must be adhered to, i.e. consistency: the conditions proposed by the modeller should not contain contradictions. 10.

(22) 1. Mathematical modelling 3. Non-dimensionalisation: The relevant quantities in the model description often have physical dimensions, e.g. time. The specific measuring system should not have an influence on the actual behaviour of these quantities. 4. Reduction: When modelling physical systems, initially most models tend to be too complex. It is therefore necessary to reduce the model to a simpler form, which is still a reliable model. 5. Asymptotic analysis: Models describing physical systems are mostly too complicated to solve analytically. By applying asymptotic analysis to the model, the modeller may try to solve the system analytically in order to find an approximation within the range of the parameter values. These results may act as an initial check of the reliability of the proposed model. 6. Numerical analysis: In this step the modeller presents a quantitative solution to the problem by means of numerical methods. 7. Verification: It is of importance that the modeller investigates the solution as a function of the independent variables and adjustable parameters, and tests these features against data. The modeller may also propose experiments in order to validate the model. A model which is too complex to be verified experimentally or contains too many parameters is useless. 8. Reiteration of steps 2–7: After one has compared the data to the model solution, one may reiterate the steps as stated, in order to improve on the proposed model. It is therefore important to follow a modular approach when modelling: first try the simplest relevant model, and add information as one progresses. 9. Implementation: The final step is where the model is implemented and the final results are used by nonmathematicians. These steps are merely a guideline in the modelling process, and not absolute. The modeller ultimately decides the process he follows by trial and error.. 11.

(23) 1. Mathematical modelling. 1.3. Model fitting. When the relations between variables and parameters have been determined and the modeller has data at his disposal, the question arises of whether the model describes the physical phenomenon accurately. The methods described here fit the proposed mathematical description (function) of the system to the data which have been collected, or which are available. The procedure is to find the parameters β that minimise the error function for a given set of n data points (xi , yi ), i = 1, 2, ..., n, with a class of functions of the form f (x; β), which describe the proposed model (Giordano et al., 2003; Velten, 2009). The components of the error function are given by . Ei (xi , β) = yi − f (xi ; β),. i = 1, ..., n. and. β1.   β2  β= .  .. .     ,  . βm where β1 , β2 , . . . , βm are the parameters in the model. Since the functions in this study are elements of the space of continuous functions with continuous n-th derivatives in R, i.e. C n (R), we may define any norm on the space, that is (C n (R), k · kp ) with p > 0 (Kreyszig, 1989). The following approximation criteria are based on those presented in Giordano et al. (2003). 1. l∞ -approximation criterion Given a function y = f (x) with unknown parameters and a collection of n data points (xi , yi ), minimize the largest absolute deviation over the entire collection. That is, determine the parameters of the function y = f (x) such that (. min E∞ (xi , β) = min. )

(24)

(25) max

(26) Ei (xi ; β)

(27). i∈[1,n]. (. = min. )

(28)

(29) max

(30) yi − f (xi ; β)

(31) .. i∈[1,n]. (1.1). 2. l1 -approximation criterion Given a function y = f (x) with unknown parameters and a collection of n data points (xi , yi ), minimize the sum of absolute deviations over the entire collection. That is, determine the parameters of the function y = f (x) such that min E1 (xi , β) = min. n

(32) n

(33) X X

(34)

(35)

(36) Ei (xi ; β)

(37) = min

(38) yi − f (xi ; β)

(39) . i=1. (1.2). i=1. 3. l2 -approximation criterion Given a function y = f (x) with unknown parameters and a collection of n data points (xi , yi ), minimize the sum of squared deviations over the entire collection. That is,. 12.

(40) 1. Mathematical modelling determine the parameters of the function y = f (x) such that min E2 (xi , β) = min. ( n X. )1/2 2. Ei (xi ; β). ( n )1/2 X

(41)

(42) 2

(43)

(44) = min yi − f (xi ; β) .. i=1. 1.4. (1.3). i=1. Dimensional analysis and scaling. When attempting to model physical phenomena, it is of importance that one includes information that is available, but no more than what is necessary. Dimensional analysis ensures that all physical information is used in modelling a system. Dimensional analysis is a technique used to restructure the original dimensional variables in a model or system into a set of dimensionless products with the result being dimensional homogeneity between the variables (Vignaux, 1992; Giordano et al., 2003). The method of nondimensionalising a mathematical model is a constructive way of formulating the model in terms of only dimensionless quantities. Nondimensionalising also yields insight into the scaling relations of the model (Van Groesen and Molenaar, 2007). The following definitions are fundamental in understanding the implications that operations performed on variables have on their dimensions. Definition 1.4.1:(Drobot, 1953) Two elements X and Y of a linear space Π have the same dimension if there exists a number ξ such that Y = ξX. We write this symbolically as [X] = [Y ] .. (1.4). The space Π can be partitioned into disjoint classes, such that each element of each class has the same dimension (Drobot, 1953). It is useful to introduce the concept of products and powers of dimensions. Definition 1.4.2:(Drobot, 1953) The product of dimensions of variables [X] and [Y ] is given by [X][Y ] = [XY ],. (1.5). and the power of the dimension of a variable [X] to a real number a is given by [X]a = [X a ].. (1.6). In particular, it follows that if a is a number, then we have [aX] = [a][X] = [X], such that [a] = 1. Therefore the dimension of any number is 1. The fundamental dimensions of physical systems are summarised in Table 1.1, adapted from the 13.

(45) 1. Mathematical modelling Physical quantity. Dimension. Physical quantity. Dimension. Mass Length Time Velocity Acceleration Specific weight Force Frequency Angular velocity Angular acceleration Angular momentum Energy. M L T LT −1 LT −2 M L−2 T −2 M LT −2 T −1 T −1 T −2 M L2 T −1 M L2 T −2. Momentum Work Density Viscosity Pressure Surface tension Power Rotational inertia Torque Entropy Heat. M LT −1 M L2 T −2 M L−3 M L−1 T −1 M L−1 T −2 M T −2 M L2 T −3 M L2 M L2 T −2 M L2 T −2 M L2 T −2. Table 1.1. Dimensions of common physical quantities. Source: Giordano et al. (2003). text by Giordano et al. (2003). Denote the dimension of time by T , the dimension of length by L and the dimension of mass by M . The following discussion was adopted from literature by Van Groesen and Molenaar (2007) and Mattheij and Molenaar (1991): Consider a system with scalar variables x1 , x2 , ..., xk and scalar parameters β1 , β2 , ..., βl . The total number of quantities involved in the system is N = k + l. Now form the product r. r. xr11 xr22 ...xrkk , β1k+1 β2k+2 ...βlrN .. (1.7). The question of interest here is for which choices of the exponents rj these products will be dimensionless. This is accomplished by replacing each variable and parameter with its respective fundamental dimensions. Say we have m fundamental dimensions involved in this system, then the replacement delivers another product ds11 ds22 ...dsmm ,. (1.8). where the numbers si , i = 1, 2, ..., m are linear functions of the exponents rj , j = 1, 2, .., N . Now by equating si = 0,. for every. i = 1, ..., m,. (1.9). we obtain a system of m linear equations for the N unknowns r1 , r2 , ...rN . It follows then that there are N − m linearly independent solutions, which correspond to N − m dimensionless quantities, say πi with i = 1, 2, ..., (N − m). This is formalised in the following theorem:. 14.

(46) 1. Mathematical modelling Theorem 1.4.1: Buckingham’s π-theorem: Van Groesen and Molenaar (2007); Mattheij and Molenaar (1991) Consider a system with variables x1 , x2 , ..., xk and parameters β1 , β2 , ..., βl in which there are m fundamental dimensions involved. Then, k + l − m dimensionless quantities πi can be identified, which are products of the above-mentioned variables and parameters. A scalar model equation f (x1 , ..., xk ; β1 , ..., βl ) = 0,. (1.10). can then be replaced by a function H(π1 , ..., πk+l−m ) = 0,. (1.11) . of dimensionless products.. This theorem enables us to develop a scheme for nondimensionalising a model. We apply the technique discussed in Van Groesen and Molenaar (2007) and Giordano et al. (2003) to a typical system of three differential equations.. Example Consider the system of equations, dX = −βXY, dt dY = βXY − νY, dt dZ = νY. dt The variables in the equations are X, Y , Z and t and the parameters are β and ν. Therefore k = 4, l = 2 and k + l = N = 6, and the fundamental dimensions of these equations are (time) and (person), such that m = 2. The parameters β and ν are both rates, but β has dimension (person)−1 (time)−1 and the parameter ν has dimension (time)−1 . Take note that the variables X, Y and Z represent the number of people within that compartment at a certain time and therefore have dimension (person). We write the product, tr1 X r2 Y r3 Z r4 β r5 ν r6 ,. 15. (1.12).

(47) 1. Mathematical modelling Replacing the variables and parameters by their respective fundamental dimensions, we find, . (time)r1 (person)r2 (person)r3 (person)r4 (time)−1 (person)−1. r5 . (time)−1. r6. ,. ⇒(time)r1 (person)r2 (person)r3 (person)r4 (time)−r5 (person)−r5 (time)−r6 , ⇒(time)r1 −r5 −r6 (person)r2 +r3 +r4 −r5 . From this we construct the equation For (time) : r1 − r5 − r6 = 0, For (person) : r2 + r3 + r4 − r5 = 0.. (1.13) (1.14). This delivers two equations in six variables, which means that there exist only two linearly independent solutions, such that four variables are free variables, i.e. r1 = r5 + r6 ,. r5 , r6 ∈ R,. r2 = r5 − r3 − r4 ,. r3 , r4 , r5 ∈ R.. (1.15) (1.16). Let us consider the different cases: • Choose r3 = 1, r4 = 0, r5 = 0 and r6 = 0. Then r2 = −1 and r1 = 0, which delivers the dimensionless product π1 =. Y . X. (1.17). • Choose r3 = 0, r4 = 1, r5 = 0 and r6 = 0. Then r2 = −1 and r1 = 0, which delivers the dimensionless product π2 =. Z . X. (1.18). • Choose r3 = 0, r4 = 0, r5 = 1 and r6 = 0. Then r2 = 1 and r1 = 1, which delivers the dimensionless product π3 = tβX.. (1.19). • Choose r3 = 0, r4 = 0, r5 = 0 and r6 = 1. Then r2 = 0 and r1 = 1, which delivers the dimensionless product π4 = tν.. 16. (1.20).

(48) 1. Mathematical modelling Using equation (1.20), proceed by substituting the value t∗ = tν, which delivers dX = −βXY, d(t∗ /ν) dY = βXY − νY, d(t∗ /ν) dZ = νY. d(t∗ /ν). (1.21) (1.22) (1.23). Simple algebra then delivers the new system, dX β = − XY, ∗ dt ν dY β = XY − Y, dt∗ ν dZ = Y, dt∗. (1.24) (1.25) (1.26). β . Note that this does not deliver a dimensionless system. The choice for using ν the dimensionless parameter π4 was done as an “informed choice”. Proceeding with scaling by a with the parameter. quantity that has dimension (person), will ensure that this system will be completely dimensionless.. Scaling Scaling of a model enables us to scale the variables in the model to values in the interval [0, 1] (Van Groesen and Molenaar, 2007; Segel, 1972). We proceed to apply a scaling technique to the same system as discussed above. Consider. dX = −βXY, dt dY = βXY − νY, dt dZ = νY. dt Scale the variables t, X, Y and Z by t∗ =. t , a. X∗ =. X , b. Y∗ =. Y , c. Z∗ =. Z , d. where the scaling factors a, b, c and d are not specified. Substitution into the system yields. 17. (1.27).

(49) 1. Mathematical modelling. dX ∗ = −acβX ∗ Y ∗ , dt∗ dY ∗ = abβX ∗ Y ∗ − aνY ∗ , dt∗ dZ ∗ ac = νY ∗ . dt∗ d Now choose the scaling factors by means of an informed guess. Scaling does not necessarily guarantee that the model will be dimensionless, but one can devise a strategy to achieve this. Firstly, based on what was done in the part on dimensional analysis, and in order to illustrate how to apply 1 the process of nondimensionalising a system, coupled with scaling of the system, choose a = . ν Secondly, by choosing the scaling factors b, c and d wisely, this ensures that the variables in the system are always in the interval [0, 1] (Segel, 1972). Choose these values as. b = max |X(t)| = N, t∈I. c = max |Y (t)| = N, t∈I. d = max |Z(t)| = N. t∈I. Note that b = c = d = N , where N denotes some saturation level for the system and X +Y +Z = N . This choice is made based on the true maximum value that each of these variables may take for all t in an unspecified interval I. Implementing these choices for the scaling factors into the system delivers the dimensionless and scaled system dX ∗ = −p1 X ∗ Y ∗ , dt∗ dY ∗ = p1 X ∗ Y ∗ − Y ∗ , dt∗ dZ ∗ = Y ∗, dt∗ X0 βN with dimensionless parameter p1 = . The scaled initial values are (X0∗ , Y0∗ , Z0∗ ) where X0∗ = , ν N Y0 Z0 Y0∗ = , Z0∗ = . N N Figure 1.1 compares the non-scaled and scaled (and dimensionless) systems. These graphs serve as an illustration as to what we achieve by scaling and nondimensionalisation of a system. Here the time interval is taken as 100 days, the value of N as 500, initial values (X0 , Y0 , Z0 ) = (499, 1, 0), and parameter values β = 1 × 10−3 and ν = 0.1.. 18.

(50) 1. Mathematical modelling. 500. 1. X (t) Y (t) Z (t). 450. 0.8. State variables (scaled). 400 350. State variables. X ∗( t ) Y ∗( t ) Z ∗( t ). 0.9. 300 250 200 150. 0.7 0.6 0.5 0.4 0.3. 100. 0.2. 50. 0.1. 0. 0 0. 20. 40. 60. 80. 100. 0. 20. t. 40. 60. 80. 100. t. Figure 1.1. Graphs illustrating the similarities between a non-scaled and scaled system. Source: Own construction. 1.5. Population growth models. In this section, the discussion is based on population growth models proposed by Thomas Malthus and Pierre François Verhulst, and the predator-prey model proposed by Alfred J. Lotka and Vito Volterra. The fundamental forms of these models play an important role in understanding population dynamics in epidemiology (Brauer and Castillo-Chavez, 2012).. Malthusian model Malthus presented an exponential growth model for the human population, from which he concluded that the population would, at some stage, exceed its capacity to grow adequate food resources (Giordano et al., 2003). The problem can be described as follows, based on the text by Giordano et al. (2003): Suppose the total population is known at an initial time, say P0 = P (t0 ), where P (t) is a function of population change over time. We are interested in describing the population size at a time t = t1 by means of the function P (t). Assume that the population may only change due to births and deaths. Factors pertaining to births are described simply by a rate β and deaths by a rate µ. The new population after a change in time is then given by P (t + δt), governed by the equation P (t + δt) = P (t) + βP (t)δt − µP (t)δt,. (1.28). or, after manipulating the equation P (t + δt) − P (t) = βP (t) − µP (t) = (β − µ)P (t) := κP (t). δt. (1.29). This implies that the change in the population over a period of time is proportional to the population. 19.

(51) 1. Mathematical modelling size. Taking the limit as δt → 0, the result is the differential equation, dP = κP, dt. P (t0 ) = P0 ,. t0 ≤ t ≤ t1 ,. (1.30). where κ > 0 for growth. By using the method of separation of variables we solve the model: Z. dP = κ dt, P ln P = κt + c0 . Z. Substituting initial values, we find ln P0 = κt0 + c0 , such that. c0 = ln P0 − κt0 .. Then, ln P = κt + ln P0 − κt0 , P = P0 eκt−κt0 , delivering the equation P (t) = P0 eκ(t−t0 ) .. (1.31). Equation (1.31) is known as the Malthusian model of population growth and predicts that a population grows exponentially. Figure 1.2 illustrates the Malthusian population growth model for an initial population of P0 = 200 with an average growth rate of κ = 0.5 over a period of 10 years.. Verhulst model A refinement of the Malthusian model is proposed by reconsidering the initial assumptions made about the population in the previous section. Consider the proportionality factor κ, which is a measure of the rate of growth of the population. Verhulst introduced the idea of refining the model to reflect limited growth in the population (Giordano et al., 2003). Assume the population increases until it approaches a maximum size, say M . In order to ensure this, the rate of growth, κ should decrease over time. A simple submodel is described by the linear function κ = ρ(M − P (t)),. 20. ρ > 0,. (1.32).

(52) 1. Mathematical modelling 4. x 10 3 2.5. Population. 2 1.5 1 0.5 0. 0. 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. Time Figure 1.2. Illustration of the Malthusian population growth model. Source: Own construction. with ρ a constant. Substitute this into equation (1.30) and find dP = ρ(M − P )P, dt. P (t0 ) = P0 ,. t0 ≤ t ≤ t1 .. Solving equation (1.33) delivers the function for the logistic growth of the population:. dP = ρ dt, (M − P )P 1 P ln = ρt + c0 . M M −P Z. Z. Substituting initial values and solving for the constant c0 delivers 1 P0 c0 = ln M M − P0 . 21. . − ρt0 .. (1.33).

(53) 1. Mathematical modelling Consequently, solving for P , P 1 P0 1 ln = ρt + ln − ρt0 , M M −P M M − P0 P P0 ln = ρM t + ln − ρM t0 , M −P M − P0 P P0 eρM t−ρM t0 , = M −P M − P0 P0 P +M −M eρM (t−t0 ) , = M −P M − P0 M P0 −1 + eρM (t−t0 ) , = M −P M − P0 M P0 eρM (t−t0 ) + 1, = M −P M − P0 P0 ρM (t−t0 ) M = (M − P ) e +1 , M − P0 P0 P0 ρM (t−t0 ) ρM (t−t0 ) e +1 −P e +1 , M =M M − P0 M − P0 P0 P0 P eρM (t−t0 ) + 1 = M eρM (t−t0 ) + 1 − M, M − P0 M − P0 P0 P0 ρM (t−t0 ) ρM (t−t0 ) P e +1 =M e , M − P0 M − P0 P0 M eρM (t−t0 ) M − P 0 , P = P0 eρM (t−t0 ) + 1 M − P0 P0 M eρM (t−t0 ) M − P0 M −P0 P = · P M − P0 0 eρM (t−t0 ) + 1 , M − P0 . . . . which delivers the function, P (t) =. P0 M eρM (t−t0 ) . M − P0 + P0 eρM (t−t0 ). (1.34). Figure 1.3 illustrates the logistic growth of a population over a period of 10 years, with an initial population of P0 = 200, a maximum population of M = 10000, at the rate ρ = 1 × 10−4 . The curve shown in Figure 1.3 is known as a logistic curve. By rewriting equation (1.34) in the form M P0 , P (t) = P0 + (M − P0 )e−ρM (t−t0 ). (1.35). it can easily be seen that as t → ∞, then P → M . It is important to realise that, when modelling population growth, the maximum population, M , is not always known. When working with data, a modeller may decide that the growth of the population is essentially represented by logistic growth. The maximum population size may then 22.

(54) 1. Mathematical modelling. 10000 9000 8000. Population. 7000 6000 5000 4000 3000 2000 1000 0. 0. 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. Time Figure 1.3. Illustration of the logistic population growth model. Source: Own construction. be estimated by considering the inflection point given by the second derivative of the function P (t), d2 P = ρP 0 (M − 2P ) . dt2. (1.36). Setting the second derivative equal to zero delivers the inflection point of a function. Simplifying d2 P = 0 delivers P = M/2. This means that when the growth dP/dt is most rapid, the population dt2 P reaches half the limiting population M before the growth rate diminishes and tends to zero.. 23.

(55) 1. Mathematical modelling. Predator-Prey model Introducing population growth models, a strong theme in mathematical epidemiology, leads to a discussion on predator-prey models. These models describe the population growth of two species, where one species is the primary food source of the other (Giordano et al., 2003). Epidemiological models describe parasitic relationships between hosts (or a number of hosts) and pathogens, revealing similarities to the modelling assumptions made in predator-prey models (Brauer and Castillo-Chavez, 2012). The following discussion is an adaptation from the text by Giordano et al. (2003): Consider two species, X and Y , X being the prey species of which the population is denoted by x(t) at time t, and Y the predator species of which the population is denoted by y(t) in the model. If we assume that the environment can support an unlimited population of species X, a rough model of exponential population growth is given by the equation dx = αx, dt. α > 0,. where α is a constant. Now, assume that the main predator species of X is Y . The growth rate of X is diminished proportionally to the total population of the predator Y . Assume that this interaction is represented by the equation dx = αx − βxy = (α − βy)x, dt with β a constant. This intrinsic growth rate, φ = α − βy of the population of prey is a linear function of the population of predators Y , and the rate decreases as the population of predators increases. The constants α and β are the degree of self-regulation of the prey X and the predation of the predators Y , respectively. Now consider the population of predators, y(t). In the absence of prey, the predators have no food. Assume that the population of Y declines at a rate proportional to the population at time t, which gives the exponential decay equation, dy = −µy, dt. µ > 0,. with µ a constant. However, in the presence of the prey species, the predator population may grow based on the. 24.

(56) 1. Mathematical modelling interactions between the two species, and therefore, dy = −µy + ηxy = (−µ + ηx)y, dt with η a constant. In this equation the intrinsic growth rate, ψ = −µ + ηx, of the population of Y is a linear function of the population of prey, X, in which the predators increase as the population of prey increases. This delivers the autonomous system of differential equations modelling the population interactions between X and Y , i.e. dx = (α − βy)x dt dy = (ηx − µ)y dt. (1.37) (1.38). The equations (1.37) - (1.38) are known as the Lotka-Volterra equations. This chapter has given insight into the process of mathematical modelling by discussing the broad framework proposed by Mattheij and Molenaar (1991) in Sections 1.1 and 1.2. It is, however, of importance that modelling a dynamic system is approached with the correct mathematical techniques, and that the described models are mathematically sound. The importance of this fact is emphasised in point 2 of Section 1.2, “Formulation of relations between variables and parameters”. In the succeeding chapter, the reader is introduced to the fundamental mathematical concepts for modelling the types of dynamical systems investigated in this study, and which are necessary to analyse these models. These systems are modelled by use of ordinary differential equations, which are by no means absolute. These systems may also be modelled by means of partial differential equations and difference equations.. 25.

(57) Chapter 2. Mathematical techniques In order to analyse mathematical models and to infer useful information from such models, mathematical instruments are needed to perform these analyses. This chapter discusses some of these mathematical techniques. Note that their discussion is by no means absolute; this merely serves to highlight the aspects which are important in our analyses.. 2.1. Systems of ordinary differential equations. Consider a system of first-order differential equations of the form.   x˙ 1 = f1 (t, x1 , x2 , ..., xn )      x˙ 2 = f2 (t, x1 , x2 , ..., xn ) . .   ..    . (2.1). x˙ n = fn (t, x1 , x2 , ..., xn ). The aim is to determine the n unknown functions x1 (t), x2 (t), ..., xn (t) which are dependent on time. The notation x˙ k denotes the derivative of the function xk (t) to time. A convenient way for writing a system of ordinary differential equations is by introducing vector notation. That is, . x˙ 1. .    x˙ =   . x˙ 2 .. ..    ,  . x˙ n. . f1 (·). . . x1. .   f2 (·)  F (·) =  .  .. .      .    x=  . x2 .. ..    ,  . and. fn (·). xn. where x, ˙ F and x are all (n × 1)-vectors, allowing us to write the system of differential equations (eq. 2.1) as x˙ = F (t, x).. 2.1.1. (2.2). Adaptive Runge-Kutta-Fehlberg method (ode45). In this study the method ode45 is used to numerically solve for systems of non-linear differential equations, since there is a lack of analytic methods for doing so. The ode45 algorithm in Matlab® implements the adaptive fourth- and fifth-order Runge-Kutta-Fehlberg method for systems of firstorder differential equations. Fehlberg derived this method from the classical Runge-Kutta equations using a fourth-order method, with five function evaluations and a fifth-order method using six function evaluations. Initially, the adapted method required a total of eleven function evaluations. Fehlberg, however, simplified the process by selecting parameters in such a way that it had the 26.

(58) 2. Mathematical techniques i. ai. ai − bi. ci. dij. 1. 16 135. 1 360. 0. 0. 2. 0. 0. 1 4. 1 4. 3. 6656 12825. 128 − 4275. 3 8. 3 9 32 , 32. 4. 28561 56430. 2197 − 75240. 12 13. 1932 7200 7296 2197 , − 2197 , 2197. 5. 9 − 50. 1 50. 1. 439 3680 845 216 , −8, 513 , − 4104. 6. 2 55. 2 55. 1 2. 8 1859 11 − 27 , 2, − 3544 2565 , 4104 , − 40. Table 2.1. Summary of the coefficients used in the implementation the Runge-Kutta-Fehlberg method. Source: Kincaid and Cheney (2002). advantage of calculating only six function evaluations for both equations (2.3) and (2.4), (Kincaid and Cheney, 2002). This method has the added advantage of automatic step size control (adaptive step size), which contributes to better accuracy. Applying the notation as stated above for the system, the two formulae are given by,. xn+1 = xn + ¯ n+1 = x ¯n + x. 6 X i=1 6 X. ai Fi ,. (2.3). bi Fi .. (2.4). i=1. The Fi s are calculated from. Fi = hf (t + ci h, x +. i−1 X. dij Fj ),. (i = 1, 2, ..., 6).. (2.5). j=1. The error in the values calculated by the two formulae, equations (2.3) and (2.4), is given by. ¯ n+1 = e = xn+1 − x. 6 X. (ai − bi )Fi .. (2.6). i=1. This equation enables us to monitor the step size. Equation (2.3) is of the fifth order and equation (2.4) is of the fourth order. The values of the coefficients are summarised in Table 2.1. While the size of the error e is monitored, the step size is adjusted as follows: If kek < is doubled, i.e. use 2h. If kek ≥ δ the step size is halved, i.e.. 27. h 2.. δ 128. then the step size.

(59) 2. Mathematical techniques. 2.1.2. Stability of systems of differential equations. In the analysis of ordinary differential equations, it is important to understand certain underlying properties of the system, i.e. when a system will be stable or unstable. Consider a system of differential equations, x˙ = f (t, x),. (2.7). where x˙ is an n×1 vector and f (t, x) is a vector with components fi (t, x1 , x2 , ..., xn ), for i = 1, 2, ...n. Assume that the fi are continuous functions having at least continuous first order partial derivatives. This implies, given initial conditions, that the system’s solution exists and is unique (Barnett and Cameron, 1975:169). If the functions fi are not explicitly dependent on t, the system (2.7) is called autonomous. Definition 2.2.1: (Mattheij and Molenaar, 1991): An equilibrium state x = x∗ is (i) stable, if for any ε > 0, there exists a δ > 0, such that kx(t0 ) − x∗ k < δ implies that kx(t) − x∗ k < ε for t ≥ t0 , (ii) asymptotically stable if it is stable, and if limt→∞ kx(t) − x∗ k = 0, (iii) unstable if it is not stable, that is, if there exists an ε > 0, such that for every δ > 0, there exists an x(t0 ) with kx(t0 ) − x∗ k < δ, kx(t1 ) − x∗ k ≥ ε for some t1 > t0 . If this is true for all x(t0 ) in kx(t0 ) − x∗ k < δ then the blank text equilibrium is completely unstable.. . When the system (2.7) is unstable, small perturbations may cause the system to follow a completely different trajectory. Before we discuss the concept of stability for nonlinear systems, first consider a linear system given by x˙ = A(t)x,. (2.8). where A(t) is an n × n matrix. If A is a constant matrix, then the system (2.8) has the solution x(t) = x0 etA , where the matrix exponent is defined as, etA = I + tA +. t2 2 t3 3 A + A + ... . 2! 3!. (2.9). Definition 2.2.2: (Strang, 1976:422-425), (Ayres, 1962:206): If an n × n matrix A has m distinct eigenvalues λi , with i = 1, 2, ..., m, of which the eigenvalues. 28.

(60) 2. Mathematical techniques each has an algebraic multiplicity θi , then the Jordan canonical form of A is given by . λ1. 1. 0. ···. 0. .              J =             . 0. λ1. 1. ···. 0. 0 .. .. 0 .. .. λ1 · · · .. . . . .. 0. 0. 0. 0.              ,             . 1. · · · λ1 ... . λm. 1. 0. ···. 0. 0. λm. 1. ···. 0. 0 .. .. 0 .. .. λm · · · .. . . . .. 0. 0. 0. 0. 1. · · · λm. . where the dimension of each Jordan block is θi .. For an n × n matrix A with n distinct eigenvalues, the algebraic multiplicity of each eigenvalue is one, and the corresponding Jordan canonical form of A is given by . λ1. 0.    J =  . 0 .. .. λ2 · · · 0 .. . . . . . . .. 0. 0. ···. 0.     .  . · · · λn. The criteria for stability of a system of differential equations are summarised in the following theorem. Theorem 2.2.1: (Barnett and Cameron, 1975:177-178) The stability of the linear system (2.8) with A a constant n × n matrix is determined by the eigenvalues, λi of the matrix A. The following cases are possible. The system is (i) asymptotically stable when Re(λi ) < 0 for i = 1, 2, ..., n , (ii) stable when Re(λi ) ≤ 0 for i = 1, 2, ..., n, provided that the dimension of the Jordan-block for λi with Re(λi ) = 0 equals 1, (iii) unstable if any Re(λi ) > 0, and (iv) completely unstable if Re(λi ) > 0 for i = 1, 2, ..., n . Proof. (Barnett and Cameron, 1975:177-178) The solution to the system (2.8) is given by x(t) = x0 eAt . The Sylvester-formula for determining matrix functions for repeated characteristic roots is given by, f (A) =. r h X. i. f (λk )ζk,1 + f (1) (λk )ζk,2 + ... + f (αk −1) (λk )ζk,αk ,. k=1. 29. (2.10).

(61) 2. Mathematical techniques where the indices αi are the multiplicity of the eigenvalues of A. The terms f (n) (λk ) denote the n-th derivative of f with respect to λ evaluated at λ = λk . The matrices ζkj are constant and determined by A. From f (λk ) = eλk t , we have At. e. =. r X. [ζk,1 + ζk,2 t + ... + ζk,αk tαk−1 ] eλk t .. (2.11). k=1. Using properties of norms we obtain keAt k ≤. αk r X X ł. k. =. tł−1 keλk t kkζk,l k,. αk r X X. tł−1 eRe(λk )t kζk,l k −→ 0, as t → ∞,. (2.12). ł. k. provided that Re(λk ) < 0 for all k. From the solution of (2.8) we see that kx(t)k ≤ keAt kkx0 k,. (2.13). so the system is asymptotically stable. If any real part of an eigenvalue of A is positive, then it is clear that kx(t)k → ∞ as t → ∞, . so the origin is unstable.. This theorem enables us to test for different stability criteria in a linear system. As a result of this theorem, Liapunov introduced the idea of linearising a system of non-linear differential equations, in order to deduce stability criteria for non-linear systems. Suppose the components of f of the system (2.7) are such that we may apply Taylor’s theorem. We obtain a new system x˙ = f (x, t) = A0 x + g(x, t),. (2.14). by using any equilibrium state of the system f (x∗ , t) = 0. The matrix A0 denotes the constant matrix. . ∂f1 ∗  ∂x (x )  1  ∂f2 ∗   ∂x (x ) 0  1 A = ..  .    ∂fn. ∂x1. (x∗ ). ∂f1 ∗ (x ) · · · ∂x2 ∂f2 ∗ (x ) · · · ∂x2 .. . ··· ∂fn ∗ (x ) · · · ∂x2. . ∂f1 ∗ (x ) ∂xn  ∂f2 ∗  (x )  ∂xn .  ..  .. (2.15).  . ∂fn ∗  (x ) ∂xn. The component functions of g(x, t) have power series expansions in (x1 , x2 , ..., xn ), beginning with terms of second degree, and the system x˙ = A0 x,. 30. (2.16).

(62) 2. Mathematical techniques is referred to as the first (linear) approximation of the system (2.7). Lemma 2.2.1: (Barnett and Cameron, 1975:214) The real matrix A is a stability matrix if and only if for any given real symmetric positive definite matrix Q the solution P of the continuous Liapunov matrix equation AT P + P A = −Q,. (2.17) . is also positive definite. Theorem 2.2.2: (Liapunov’s linearisation theorem), (Barnett and Cameron, 1975:216). If the system (2.16) is asymptotically stable, or unstable, then the equilibrium point x∗ has the same property for the non-linear system in (2.7) Proof. (Barnett and Cameron, 1975:217) Consider the quadratic form V = xT P x,. (2.18). (A0 )T P + P A0 = −Q,. (2.19). where P satisfies. Q being an arbitrary real, symmetric, positive-definite matrix. If the system (2.16) is asymptotically stable then by Lemma 3.2.1 the matrix P is positive definite. The derivative of V with respect to the right hand side of (2.14) is V˙ = −xT Qx + 2g T P x.. (2.20). Because of the nature of g, the second term in the above equation has degree three at the very least, we have for x sufficiently close to x∗ that V˙ < 0. Hence the equilibrium point x∗ . is asymptotically stable.. 2.2. Parameter Sensitivity Analysis. In modelling physical systems it is of importance to determine which parameters in the model are most influential on state variable outputs (Hamby, 1994). Hamby (1994) highlights five reasons for conducting parameter sensitivity analysis, i.e., 1. Which of the parameters in the model require additional research in order to broaden the knowledge base. 2. Which of these parameters play an insignificant role and may be excluded in the final model. 31.

No results found