The Mayaro virus and its potential epidemiological consequences in Colombia: an exploratory biomathematics analysis

Mayaro virus (Togaviridae) is an endemic arbovirus of the Americas with epidemiological similarities with the agents of other more prominent diseases such as dengue (Flaviviridae), Zika (Flaviviridae), and chikungunya (Togaviridae). It is naturally transmitted in a sylvatic/rural cycle by Haemagogus spp., but, potentially, it could be incorporated and transmitted in an urban cycle by Aedes aegypti, a vector widely disseminated in the Americas. The Mayaro arbovirus dynamics was simulated mathematically in the colombian population in the eight biogeographical provinces, bearing in mind the vector’s population movement between provinces through passive transport via truck cargo. The parameters involved in the virus epidemiological dynamics, as well as the vital rates of Ae. aegypti in each of the biogeographical provinces were obtained from the literature. These data were included in a meta-population model in differential equations, represented by a model structured by age for the dynamic population of Ae. aegypti combined with an epidemiological SEI/SEIR-type model. In addition, the model was incorporated with a term of migration to represent the connectivity between the biogeographical provinces. The vital rates and the development cycle of Ae. aegypti varied between provinces, having greater biological potential between 23 °C and 28 °C in provinces of Imerí, biogeographical Chocó, and Magdalena, with respect to the North-Andean Moorland (9.33–21.38 °C). Magdalena and Maracaibo had the highest flow of land cargo. The results of the simulations indicate that Magdalena, Imerí, and biogeographical Chocó would be the most affected regarding the number of cases of people infected by Mayaro virus over time. The temperature in each of the provinces influences the local population dynamics of Ae. aegypti and passive migration via transport of land cargo plays an important role on how the Mayaro virus would be disseminated in the human population. Once this arbovirus begins an urban cycle, the most-affected departments would be Antioquia, Santander, Norte de Santander, Cesar (Provinces of Magdalena), and Valle del Cauca, and Chocó (biogeographical province of Chocó), which is why vector control programmes must aim their efforts at these departments and include some type of vector control to the transport of land cargo to avoid a future Mayaro epidemic.

Brazil, Venezuela, Peru, Bolivia, Surinam, Costa Rica, Guatemala, Panama, Mexico, Haiti, and, sporadically, in Colombia [3,4]. Since its discovery, three genotypes, D, L and N, have been identified; one of them (L) is exclusive of the Amazon region of the State of Pará, Brazil. The second genotype (D), has been isolated from the Amazon region of Peru, Bolivia, Venezuela, Colombia, Argentina and Brazil, and the third genotype (N) was isolated from a strain from Peru [5,6].
The disease caused by the Mayaro virus is characterized by high fever and intense joint pain, with duration from 3 to 7 days, causing medical disability [7]. Most infections with Mayaro virus are sporadic and occur in people with a recent history of activity in forest areas (rural or semirural) or their surroundings [8]. This disease is not fatal, but it is unknown whether it generates other symptoms, like Zika, which is related with microcephaly in neonates and Guillain-Barré syndrome [9]. The Mayaro virus shares characteristics such as the special tropism through osteomuscular tissue with more prominent arbovirus like dengue (Flaviviridae), Zika (Flaviviridae) and chikungunya (Togaviridae) [10]. In addition, it shares epidemiological, ecological and biogeographic traits, including dependence on vectors at the forest and/or urban level [11][12][13]. However, unlike dengue (Asian), chikungunya and Zika (African) viruses, it is native to the American continent [8].
The Mayaro virus is endemic in rural areas in South America where it maintains an enzootic sylvan cycle (similar to chikungunya) that involves wild vertebrates (birds and reptiles), including primates and hematophagous mosquitoes [14]. In sylvatic and rural areas, Mayaro is transmitted by the bite of females of the Haemagogus spp. infected with this virus [15]. Furthermore, in natural populations of Aedes (Stegomyia) aegypti (Linnaeus, 1762) from South America and North America it has been observed experimentally their susceptibility to being infected and transmitting the genotype D of the Mayaro virus [16,17]. Additionally, in natural populations of Ae. aegypti from South America, natural infection with the genotype L of the Mayaro virus has also been observed [18].
Aedes aegypti has the human population as preferential source of food, and urban areas as main reproduction areas [19], although it can also be found in rural areas [20]. Consequently, this species has the potential of serving as vector bridge between rural areas close to natural foci of the Mayaro virus and urban areas [8,[21][22][23]. Due to the aforementioned, it cannot be discarded that this arbovirus can emerge as a global pathogen through Ae. aegypti, as recently observed in chikungunya and Zika arboviruses [24,25].
The importance of Ae. aegypti lies in that it is present in large densities in all continents, with a high degree of anthropophilia, thus, being able to introduce the Mayaro virus into an urban cycle and, thereby, rapidly spread the etiological agent [21,26]. In the Americas, Ae. aegypti is the primary vector of dengue, chikungunya and Zika [27][28][29], which are distributed in the tropical regions of the continent [21]; in Colombia, this species is present up to 2300 m [30], representing 80% of the territory [9] that could be at risk of contracting the Mayaro virus.
Once the Mayaro virus starts to circulate in the urban area through Ae. aegypti and, given that currently no vaccine is available to prevent Mayaro infection in humans (with good results in mice) [31], the only alternative will be the control of this vector. For this purpose, elimination or closing of locations where larvae develop and use of insecticides to control larval and adult forms are ways currently used globally [28]. Nevertheless, frequent use of insecticides has caused resistance to practically all classes of insecticides used globally [32], which is why substitution options by other classes of insecticides are scarce due to the poor offer of environmentally safe products with distinct action sites thus hindering actions to control of this species [33]. Hence, understanding the Mayaro epidemiological cycle in the urban area is of vital importance to elaborate control programmes more in keeping with these new problems.
Given that the Mayaro virus is genetically related to the chikungunya virus [34], and chikungunya is transmitted by Ae. aegypti in Colombia, an epidemiological behavior would be expected for Mayaro virus similar to that observed in chikungunya. Recently, it has been observed experimentally that in Ae. aegypti individuals infected initially with chikungunya and then with Mayaro virus, both arboviruses can be transmitted with similar efficiency. Further, if it is initially infected with Mayaro virus and then with chikungunya virus, there is an exclusion due to over-infection to chikungunya; in any of the situations described, Mayaro virus will always be transmitted [35].
Since the natural distribution of the Mayaro virus is unknown in sylvatic areas and how it could be potentially transmitted by another vector different from Haemagogus spp., an option to know its epidemiological impact is to formulate a mathematical model that permits simulating its behavior and expansion in the human population, through the Ae. aegypti potential vector, given that it is distributed in peri-urban and urban areas.
For Ae. aegypti, the mathematical models have proven to be useful tools to understand dengue transmission [14,36] and, recently, chikungunya [37] and Zika [38], as well as to help in planning control strategies [39]. Among the different mathematical approaches to study infectious diseases through vectors, as with the Mayaro virus, the SEIR-type (susceptible, exposed, infected and recovered) epidemiological models [40,41] and in metapopulations, have been used widely [42][43][44][45][46]. Nonetheless, none of them considers the passive vector transport (dispersion process of invasive species associated with human activities) or its dynamic populations, under different biogeographical conditions; therefore, an option is to include a meta-population model in differential equations, where each local population includes a structured model on age for the dynamic population of Ae. aegypti and an SEI/SEIR-type epidemiological model for the population of humans, as well as the passive transport of the mosquito through land cargo, which, has also been used recently to study the dynamics of Rift Valley fever transmission in the human population, with Ae. aegypti as the main vector [47]. Among the causes for this heterogeneity, there is the level of connectivity among the places where the species is found and the set of biogeographical conditions where the individuals inhabit [48]. For Colombia, according to Morrone [49], there are eight biogeographical provinces (North-Andean Moorland, Cauca, Napo, Imerí, Chocó, Venezuelan Plains, Maracaibo and western Ecuador), each characterized by having different geographical, climatic, altitudinal and natural components. Among these characteristics, temperature plays an important role in the life-cycle and vector transmission of mosquitoes of medical importance [50,51]. For instance, optimal temperatures for development, longevity, and fecundity are between 22 °C and 32 °C [52], while at higher temperatures survival range of Ae. aegypti decreases; however, egg-laying time increases causing a growth in egg number [53]. Moreover, the extrinsic incubation period of the dengue virus is reduced, resulting in higher rates of viral transmission [54,55]. On the other hand, at lower temperatures, the transmission of arbovirus may be impeded in some cases [22]. Considering that temperature and passive transport are important factors in arbovirus transmission, here we evaluated how passive transport and the different vital parameters of Ae. aegypti affect the dynamics of Mayaro virus transmission in the Colombian human population, from which the study only quantified transport via land cargo due to the lack of data from other means of transport.

Methods
The parameters involved in the epidemiological dynamics of the Mayaro virus were obtained from the literature; data used here belong to Mayaro genotype D (Table 1). In a parallel manner, based on the biogeographical regions for Latin America proposed by Morrone [49], Colombia was divided into eight biogeographical provinces (Fig. 1). For each biogeographical province, from historical data available in the Wolrdclim-2 library on temperature from 1970 to 2000 [56], the minimum and maximum median historical temperature was estimated. These variables are recognized in the literature as the most influential in the life-cycle of Ae. aegypti [50,51]. The areas of the biogeographical provinces were defined for Colombia, as proposed by Morrone [49] in R software version 3.5.1, using the raster, rgdal [57], sp [58], dplyr [59] and st [60] libraries. Thereafter, the layers of the biogeographical provinces were overlapped with each of the raster layers for temperature. Then, for each biogeographical province and each raster layer of temperature, the median, minimum, and maximum historical average temperatures were extracted. At the same time, the vital rates of Ae. aegypti corresponding to different temperatures were obtained from the literature and adjusted by average for each region, taking into account its temperature range. Additionally, in 2018 and for each province, the size of the human population was estimated through the projections of the Departamento Administrativo Nacional de Estadística (DANE; Table 1 These data were included in a meta-population model in differential equations, where each province was considered a local population. For each province, a SEI-type model structured on age for the population dynamics of Ae. aegypti was combined with a SEIR-type epidemiological model for the human population. In addition, the model included a term that represents the passive transport of the mosquitoes through land cargo in Colombia, using data from the Colombian Ministry of Transportation for 2017 [62]. Figure 2 represents the compartmental diagram of the model. In the model structured on age for the population dynamics of Ae. aegypti, four stages were considered: egg ( H ); larva ( L ); pupa ( P ); and adult (A). The adult stage is divided, in turn, into three classes, based on the SEI epidemiological model: susceptible ( As ), exposed ( Ae ), and infected ( Ai ). The vital rates of Ae. aegypti in each stage vary according to the biogeographical province ( r ). Time (t) is measured in days.
The variables and parameters used are defined by: H r (t) : number of eggs over time t , in province r ; L r (t) : number of larvae over time t , in province r ; P r (t) : number of pupas over time t , in province r ; As r (t)r : number of adults susceptible over time t , in province r ; Ae r (t) : number of adults exposed over timet , in province r ; Ai r (t) : number of infectious adults over time t , in province r ; θ : proportion of females/males; K : carrying capacity for larvae; χ : intraspecific competition coefficient, proportional to K ; µ r x : death rate of individuals of stage x ( x = H , L, P and A ); r H : rate of transition from egg to larva, in province r ; r L : rate of transition from larva  [52,78,87] to pupa, in province r ; r P : rate of transition from pupa to adult, in province r ; r A : rate of oviposition of female mosquitoes, in province r ; Ae : rate of transition from exposed to infected adult mosquitoes; m r : net number of adult migrant mosquitoes per day in province r ; β H : rate of effective contacts between susceptible mosquitoes and infected humans; β M : rate of effective contacts between infected mosquitoes and susceptible humans.
The system of differential equations that interprets the dynamics of Ae. aegypti is the following: The human population was divided into four classes: susceptible ( S ); exposed ( E ); infected ( I ); and recovered dP r dt = r L L r − r P + µ r P P r dAs r dt = r P P r − β H (I r /NH r )As r − µ r A As r + m r As r dAe r dt = β H (I r /NH r )As r − Ae Ae r − µ r A Ae r + m r Ae r dAi r dt = Ae Ae r − µ r A Ai r + m r Ai r NM r = As r + Ae r + Ai r ( R ) over time t , (in days). The vital rates of the human population were considered the same for all the biogeographical provinces. The variables and the parameters used in this model are: S r (t) : number of humans susceptible over time t , in province r ; E r (t) : number of humans exposed over time t , in province r ; I r (t) : number of humans infected over time t , in province r ; R r (t) : number of humans recovered over time t , in province r ; E : rate of transition from humans exposed to humans infected; I : rate of transition from infected to recovered humans; µ r F : death rate of humans in class z(F = S, E, R, I) ; N r : total human population in province r ( N r = S r + E r + I r + R r ).
The system in differential equations that interprets the SEIR dynamics for humans is: The Ae. egypti passive transport through the movement of cargo trucks is represented with the term m r A r , where m r is the net daily rate of migration of adult mosquitoes between provinces. Assuming that for every 1000 trucks traveling from one province to another, the number of effective migrants (Mig ) is equal to 3 mosquitoes, as already observed in the field [63]. These numbers are consigned in Table 4, obtained in turn from official data on cargo transport in Colombia from the National Registry of Land Cargo Dispatches (RNDC, for the term in Spanish) (Table 3), and comprise the Mig matrix of the 8 × 8 order, whose ij component represents the number of individuals going from province i to province j.
The entry rate of adult individuals in state X (As, Ae, Ai) to province r from other provinces is given by: The exit rate of adult individuals is state X from province r to other provinces is given by: where PropX r = X r PTA r is the proportion of adult individuals in state X in province r , i.e. the ratio of the number of adult individuals in state X r and the total population of adult mosquitoes (PTA r ).
Mig(r, i) * PropX ir The numerical solution of the model was obtained by using the Runge-Kutta 4th order method implemented in MATLAB ® environment [64]. From this, and using the parameters obtained from the literature for the life-cycle dynamics of Ae. aegypti and the Mayaro virus in each biogeographical province, the behavior of the vector and human populations was simulated in each of the biogeographical provinces, with and without passive transport (migrations), for 300 days. The province of Imerí was taken as center of origin of the infection with Mayaro virus of adult individuals of Ae. aegypti, given that previous reports have been registered in vertebrates for this arbovirus [65]. Additionally, scenarios were simulated with the beginning of the infection, not in Imerí, but in Magdalena, taking into account that it is the province most connected by land [62]. In both simulations, initial proportions of human populations were taken from DANE data (Table 1), and the initial mosquito populations were hypothetical values because in Colombia and the rest of the world there are no precise estimates of these, there are some mathematical estimates that cannot be extrapolated for our case [66,67]. In all provinces, the initial populations were H: 300, L: 200, P: 200, As: 200 and Ae: 0. Relating to Ai, an infected mosquito was initially placed in the infection focus region and zero in the others.
A global sensitivity and uncertainty analysis was performed using the partial coefficient correlation index (PRCC). The parameters considered for this analysis were the infection rate from the mosquito to humans ( β M ), from humans to the mosquito ( β H ) the migration rate ( Mig ), and carrying capacity ( K ) assuming a uniform distribution of the parameters and statistical independence between them. For the uncertainty analysis, the multidimensional parametric space of the parameter values was explored using the Latin hypercube sampling (LHS) [68] method in the LHS package [69] using R software. The effect of the variation was analyzed of the parameters on the model response variable, infected adult mosquitoes (Ai), through a sensitivity test [70], by means of a partial correlation coefficient analysis (PRCC) implemented with the sensibility package [71] in R. Additionally, to determine the population dynamics of Ai for those parameters significantly correlated with it, from the data reported, they were simulated taking into account its respective interval scale.

Results
From the data of the life-cycle dynamics of Ae. aegypti, associated with data for temperature from the biogeographical provinces obtained from the literature, it is possible to observe that the provinces with the best conditions for the development of this mosquito species presenting the ranges of most optimal temperatures for the development and life-cycle of Ae. aegypti, are Imerí, Chocó and Magdalena (Table 1). Table 2 shows the parameters involved in the epidemiology of the Mayaro virus, both in the vector such as in humans, obtained also from the literature. These parameters are not related directly with the biogeographical provinces, but with the conditions of the hosts.
Data on cargo transport movement available in the RNDC 2017 showed that the provinces with the highest flow of land cargo are Magdalena, Maracaibo and North Andean Moorlands, in contrast with Imerí and Napo which have the lowest connectivity (Table 3). Consequently, the same pattern was observed regarding the daily migration of mosquitoes between biogeographical provinces (Table 4). Figures 3, 4, 5 correspond to some numerical tests obtained from the model proposed in this study, considering (I) or not (II), passive transport through land cargo and when the focus of infection is in the Imerí region (Figs. 3a, b, 4a, b, 5a, b), which is the province with the lowest connectivity, or in the province of Magdalena (Figs. 3c, d, 4c, d; 5c, d), which is the province with the highest connectivity.
To appreciate the population dynamics corresponding to the different vital rates of each province, Fig. 3 shows the evolution of the proportion of adult mosquitoes in relation to the total population over time in each of the provinces. When passive transport (II) is not considered, the highest values were those of the Venezuelan Plains where the highest transfer rates of immature stages occur, and the lowest values corresponded to the Northern Moorland Province, the province with the highest mortality rate; these dynamics are independent of the infection focus province (Fig. 3b, d). When considering passive transport (I) there were small changes in the dynamics within these regions due to the arrival of individuals from other provinces; similar results were obtained for both scenarios.
To demonstrate the influence of passive transport on the spread of the virus, some simulations are presented here corresponding to scenarios where initially there are only infected mosquitoes in a region; furthermore, the only source of infection considered is due to migrations. Two infection focus regions, Imerí and Magdalena, were chosen, as they had, respectively, the lowest and highest road connectivity with the remaining regions. Figure 4 corresponds to proportions of infected mosquito populations when the focus of infection is Imerí (Fig. 4a, b) or Magdalena (Fig. 4c, d) and when there is passive transport (I) or it is not considered (II). When migrations are not considered, Fig. 4b and d correspond to the local dynamics in Imerí and Magdalena, respectively; in these cases, the rest of the regions keep the infected population at zero. The most notable difference between the two provinces is the time it takes populations to start growing and the size of the peaks, which is due to the difference between the vital rates of each region. On the other hand, if there is passive transport, infected mosquitoes begin to migrate from the infection focus region, Imerí (Fig. 4a) or Magdalena (Fig. 4c), to neighboring regions and spread according to connectivity between provinces (Table 3) and they are developed following local vital rates (Table 1). Figure 5 shows the dynamics of infected human populations when the focus of infection is Imerí (Fig. 5a, b) and Magdalena (Fig. 5c, d) and when there is passive transport (I) or not (II). Initially there are no infected humans, so, the dynamics depends on infected mosquitoes. Figure 6 shows the results of the sensitivity analysis made with the parameters human-to-mosquito contagion rate ( β H ), mosquito-to-human contagion rate ( β M ), number of mosquitoes transported per 1000 trucks ( Mig ) and larvae carrying capacity ( K ). The parameters β H (positive) and β M . (negative) were the only statistically significant and therefore are the parameters that influence this study.
Knowing that β H and β M are sensitive and significant parameters in the results of the simulations, Fig. 7 shows in the province of Imerí (Fig. 7a, b, e and f ) and Magdalena (Fig. 7c, d, g and h), the population dynamic of number of infected humans and the proportion of infected  [31,107] Intrinsic time of virus incubation 3-11 (days) [108] Extrinsic time of virus incubation 7 (days) [16,35] Rate of vector-host effective contact β M 0.5-0.8 [109] Rate of host-vector effective contact β H 0.29-0.39 [109] Larvae carrying capacity K. 50,000 [41,75] (Fig. 7a, c, e and h) and β M (Fig. 7b, d, f and g) parameters. In Imerí and Magdalena providences, the population dynamics for infected humans (Fig. 7a, c) and infected mosquitoes (Fig. 7e, h) in the simulations with respect to β H parameter showed that with β H less than 0.5, both the number of infected humans (Fig. 7a, c) and the proportion of infected mosquitoes (Fig. 7e, h) decrease. In contrast, simulations with different values o. f showed no difference in the behavior of infected humans at Imerí (Fig. 7b) or Magdalena (Fig. 7d) as well as infected mosquitoes at Imerí (Fig. 7f ) or Magdalena (Fig. 7g).

Discussion
Our principal findings indicate that once the Mayaro virus starts to disseminate through Ae. aegypti in the vector and human populations, regardless of the center of origin of the infection, Imeri or Magdalena, the provinces of Colombia most affected by the virus will be Magdalena, biogeographical Chocó, and Imerí, having as a basic difference the time in which the infection begins to grow.  and to the colonization of mosquitoes of epidemiological importance such as Ae. aegypti, the main vector in the Americas for dengue, chikungunya and Zika [75] and Ae. albopictus, the main vector of these arboviruses in Asia and Europe [27,76]. In Ae. aegypti, it has been noted that passive migration occurs principally through human activities, such as movement of land, air, sea and riverine cargo at local, regional and international levels [77,78], which evidences the urgency to quantify the effect of passive transport on migration [78][79][80]. Recently, for Ae. albopictus, it has been shown that adult individuals are transported passively in private cars and that, for every 1000 vehicles, between two and 11 carry the mosquito [63]. This agrees with the findings in the present study, where the number of individuals for migration to be effective between biogeographical provinces was between one and three mosquitoes for every 1000 trucks. Migration of Ae. aegypti, through land cargo transport, influenced directly the migration of this vector and the potential dissemination of the Mayaro virus among biogeographical provinces, principally in Magdalena, with the particularity that when the infection origin center is Magdalena, the Mayaro would spread faster to the other biogeographical provinces compared to when the infection origin center is Imeri. In Colombia, the departments of Antioquia, Cundinamarca, Santander, Norte de Santander and Huila from the province of Magdalena, are the most interconnected and, hence, with the greater probability of migration of mosquitoes infected with the Mayaro virus or not, through cargo transport [62]. Furthermore, changing climate conditions, like global warming, are generating new niches (in latitude, longitude and altitude) and optimal for subsequent habitat expansion not only in Ae. aegypti (temperature and precipitation play an important role in the abundance of this species [81]), but also in other species of medical importance [82,83]. In tropical countries, it has been found that the duration of the development cycle and the vital rates of Ae. aegypti vary with respect to temperature [84][85][86], where temperatures between 22-33 °C are optimal for its vital development [52], in contrast to temperatures ≤ 16 °C [87]. Consequently, places with optimal temperatures have population growth and shorter life development for Ae. aegypti, as well as increased probability of arbovirus transmission [88], in contrast to places with lower temperatures [87]. Our results show   ). Additionally, the carrying capacity (K) did affect the population size of living beings; however here, K did not affect the population size of Ae. aegypti, in these or in the other biogeographical regions, like the size [89] and amount of available water in the larva breeding sites [90], or climatic factors such as humidity and precipitation [91], could be affecting our results; hypotheses that were not tested here. Meanwhile, it should be stressed that these areas in Colombia have the highest levels of precipitation [92]. According to García et al. [93], areas with high levels of precipitation are characterized by higher population densities of Ae. aegypti than areas with less rainfall. According to Mohammed & Chadee [94] temperatures ≥ 34 °C cause a gradual drop in the vital rates of Ae. aegypti, with temperatures > 39 °C being lethal for this species [95,96], suppressing embryo development and diminishing survival of larval stages [97], having less individuals of the vector in the long term. The only province with temperatures over 34 °C was the Venezuelan Plains, where the vital rates of egg development time, larva, pupa and days of survival of Ae. aegypti adults are lower with respect to the other provinces.
The model presented here has limitations, which originate from the assumptions considered, among which it was assumed that Mayaro disease is caused by the virus genotype D, the most frequent, and it is dispersed from where there are records of disease outbreaks (Imerí Province [5,6]), by mosquitoes that are transported by land cargo. However, an additional limitation is the lack of data on Mayaro outbreaks in the other local populations. Despite this, the epidemic dynamics obtained in Colombia is consistent with that reported in other studies available in the literature of dengue [65,66], chikungunya [67] and Zika [66] for the country.
The HLS/PRCC analysis shows that the contact rate of infected humans and susceptible mosquitoes ( β H ) and the contact rate of infected mosquitoes and susceptible humans ( β M ) have an effect on the variation in the number of infected mosquitoes over time; this behavior has been previously evidenced by Requena & Juárez [40] in their model for chikungunya suggesting that the most effective way to control or reduce the rate ( β H ) is through the use of repellents, mosquito nets, or fumigations at Ae. aegypti hatcheries. In this scenario, population control of this vector is essential to prevent a possible Mayaro epidemic; however, currently there are difficulties in controlling Ae. aegypti populations due to the spread of resistance to the insecticides used to control it [98,99]. Therefore, in Colombia, new alternatives for the control of Ae. aegypti populations are being evaluated, such as the use and release of mosquitoes infected with Wolbachia pipientis-WMel lineage (mosquitoes with refraction to arboviruses transmitted by Ae. aegypti) [100]. Additionally, the potential vectorial role of other species such as Ae. albopictus, a species that shares ecological niches in Colombia and America with Ae. aegypti and that in other latitudes is the main vector of dengue, chikungunya and Zika is unknown [101][102][103].

Conclusions
Our results indicate that temperature influences local vector proliferation, since vital rates depend on it. That is, in each biogeographical province the dynamics of the vector is different, while the passive migration of Ae. aegypti plays an important role on how the Mayaro virus would be disseminated in human populations. Once this arbovirus begins an urban cycle from the Imerí or Magdalena provinces, the most-affected departments would be Antioquia, Santander, Norte de Santander, Cesar (Provinces of Magdalena), and Valle del Cauca and Chocó (Province of biogeographical Chocó) as already occurred with dengue, chikungunya and Zika. Therefore, vector control programmes must aim their efforts at these departments and include some type of vector control to the transport of land cargo to avoid a future Mayaro epidemic. Abbreviations A: Adult; Ae: Exposed adult mosquito; Ai: Infected adult mosquito; As: Susceptible adult mosquito; E: Egg; L: Larva; P: Pupa; RNDC: National Registry of Land Cargo Dispatches; SEI: Susceptible, exposed and infected; SEIR: Susceptible, exposed, infected and recovered.