A Landscape-Level Assessment of Restoration Resource Allocation for the Eastern Monarch Butterfly

The Monarch butterfly eastern population (Danaus plexippus) is in decline primarily due to habitat loss. Current habitat restoration programs focus on re-establishing milkweed, the primary food resource for Monarch caterpillars, in the central United States of America. However, individual components of the Monarch life cycle function as part of an integrated whole. Here we develop the MOBU-SDyM, a migration-wide systems dynamics model of the Monarch butterfly migratory cycle to explore alternative management strategies’ impacts. Our model offers several advances over previous efforts, considering complex variables such as dynamic temperature-dependent developmental times, dynamic habitat availability, and weather-related mortality across the entire range. We first explored whether the predominant focus of milkweed restoration in the mid-range of the Monarch’s migration could be overestimating the Monarch’s actual habitat requirements. Second, we examined the robustness of using the recommended 1.2–1.6 billion milkweed stems as a policy objective when accounting for factors such as droughts, changes in temperature, and the stems’ effective usability by the Monarchs. Third, we used the model to estimate the number and distribution of stems across the northern, central, and southern regions of the breeding range needed to reach a self-sustainable long-term Monarch population of six overwintering hectares. Our analysis revealed that concentrating milkweed growth in the central region increases the size of the overwintering colonies more so than equivalent growth in the south region, with growth in the northern region having a negligible effect. However, even though simulating an increase in milkweed stems in the south did not play a key role in increasing the size of the overwintering colonies, it plays a paramount role in keeping the population above a critically small size. Abiotic factors considerably influenced the actual number of stems needed, but, in general, our estimates of required stems were 43–91% larger than the number of stems currently set as a restoration target: our optimal allocation efforts were 7.35, 92, and 0.15% to the south, central, and northern regions, respectively. Systems dynamics’ analytical and computational strengths provided us with new avenues to investigate the Monarch’s migration as a complex biological system and to contribute to more robust restoration policies for this unique species.


INTRODUCTION
There is a broad agreement that habitat loss is one of the most critical causes of species decline worldwide (Andrén, 1997), particularly with habitat species where life history is tightly bound to a limited number of habitat elements, such as specific host plants for breeding or feeding (Matthews et al., 2014). However, in many such species, knowledge of host-plant phenology is still not at the level needed for effective management. As a result, it is common practice to simplify the underlying dynamics driving host-plant phenology and assume that host availability is relatively constant. This assumption is potentially risk-prone, given the effects of climate change. This paper explores the implications of such assumptions for the eastern North American Monarch Butterfly, Danaus plexippus (referred from here on as "Monarch").
The eastern North American population of the Monarch Butterfly has the longest migration of any insect (Carlos Galindo-Leal, 2005). Every winter, individuals across the Eastern United States and Canada travel up to 4,000 km south to establish highly concentrated overwintering colonies within a few forest patches in Mexico, concentrated in the Monarch Butterfly Biosphere Reserve (Urquhart and Urquhart, 1976;Santiago et al., 2001). The high concentration of overwintering individuals allows the estimation of their population size by measuring the total area occupied by the overwintering colonies . Although this estimate is considered inaccurate (Hristov et al., 2019), it has been the primary tool for designing Monarch's conservation strategies. Critically, the total area of the Monarch overwintering colonies in Mexico has steeply decreased over the past two decades The Monarch's primary conservation challenge is habitat loss due to its specialized habitat needs, both during the breeding and the overwintering stages (Flockhart et al., 2015). On the overwintering grounds, canopy integrity is critical, and its loss, due primarily to illegal logging Vidal et al., 2013), exposes the Monarchs to sub-freezing temperatures  that reduce survival (Anderson and Brower, 1996). However, since the Monarch Biosphere Reserve's core zone has been virtually free from illegal logging for the last several years (WWF-Mexico, 2019), the breeding range currently receives the most attention and conservation efforts (Thogmartin et al., 2017a). While breeding and feeding in the United States and Canada, the Monarch requires patches of land with plants from the Asclepias genus, mainly Asclepias syriaca, commonly known as milkweed (Malcolm and Brower, 1986;Batalden and Oberhauser, 2015), linking their dynamics to milkweed habitat availability during this stage.
The leading hypothesis explaining the Monarch's decline is the loss of milkweed, driven by changing farming practices across the United States and Canada breeding areas. Milkweed grows naturally in disturbed areas such as ditches, abandoned land, side roads, and farm outskirts (Jeffery and Robison, 1971;Oberhauser et al., 2008). Evidence shows that milkweed abundance has decreased considerably over the last two decades due to changes in agricultural methods, including the introduction of genetically modified crops coupled with specialized herbicides such as glyphosate (Pleasants and Oberhauser, 2012). Researchers have cited this decrease as the primary cause of the Monarch's decline (Hartzler and Buhler, 2000;Hartzler, 2010;L. Brower et al., 2012). As a response, the "United States Cornbelt" encompassing the states of Iowa, Illinois, Indiana, southern Michigan, western Ohio, eastern Nebraska, and south Minnesota has received most of the attention for habitat restoration efforts. The main argument for this restoration approach is that the United States Cornbelt provides the natal habitat for most Monarchs overwintering in Mexico (Pleasants and Oberhauser, 2012;Flockhart et al., 2013;Green et al., 2018), is the region with the highest concentration of GMO farmland, and it has the least amount of remaining native plant composition (Holt, 2017).
In this paper, we examine two current assumptions considered in the design of Monarch conservation strategies. First, following classic exponential growth, we hypothesize that first-generation breeding success across the southern part of the range may be more critical to the demographics of subsequent generations leading to the final pre-overwintering population, compared with the mid-section of the breeding range as it is commonly assumed (Crewe et al., 2019). Mortality and lack of reproductive activity while overwintering (L. P. Brower and Missrie, 1999) translate into the lowest Monarch abundance during the early spring, i.e., when Monarchs arrive at the southern edge of their breeding range from the overwintering grounds. Not disregarding central North America's importance as a critical Monarch breeding region, that early-spring Monarch's low abundance suggests that the southern region could play a pivotal role in the cycle's success. Second, recovery strategies must consider that Monarchs will not use the entirety of the stems provided, either due to limiting factors such as droughts and temperature changes or by idiosyncratic factors related to the stem's appeal to the Monarch. In other words, it is crucial to understand the difference between available and total stems and that a successful strategy will need to provide a considerably larger number of total stems to reach the Monarch's required available stems. With that in mind, we hypothesize that the number of 1.3-1.6 billion new milkweed stems per annum needed to reach a minimum viable population size falls short of the actual number the Monarch needs to establish a self-sustaining population successfully (Pleasants, 2017;Thogmartin et al., 2017b).
To test our hypotheses, we developed the Monarch Butterfly Systems Dynamics Model (MOBU-SDyM). The MOBU-SDyM is a full-migration, temperature-dependent systems dynamics model of the Monarch butterfly. Through a Bayesian-driven estimation approach, we used this model to test the hypotheses presented above and estimate the yearly number of stems that would need to be added across the breeding region to reach a self-sustainable long-term population of six overwintering hectares. Throughout the analysis, we highlight the differences between considering total versus available milkweed stems. To our knowledge, three simulation models of the Monarch's full migration exist (Yakubu et al., 2004;Flockhart et al., 2015;Oberhauser et al., 2017) but none of these models explicitly incorporate tightly coupled biotic and abiotic components of the system. Innovatively, we investigate the effects of temperature and humidity on milkweed production and the temperaturedependent behavior of Monarch reproductive capacity and intergenerational time across the landscape, all critical factors affecting Monarch population dynamics.

METHODS
We first describe how the MOBU-SDyM captures the population dynamics of the Monarch's migration, including the model's initial setup and the main equations and assumptions. Next, we focus on performance metrics to assess the model's capacity to capture the real-life system. After that, we present a sensitivity analysis that tests the leverage 1 that variations in the number of milkweed stems have over the overwintering colonies' geographic area. Finally, we describe our optimization-driven process to estimate the most significant yearly milkweed allocation across the breeding regions. The last two steps underscore the importance of differentiating total and available milkweed stems. The TRACE documentation (Grimm et al., 2014) of  [Gossard and Jones (1977), Zalucki (1981)

Overview of the MOBU-SDyM
The MOBU-SDyM is a spatially explicit stage-structured model with four geographic regions (Breeding regions in the South, Central, and North and one Overwintering region in Mexico) and three Monarch stage classes (Immature, Breeding Adults, and Overwintering Adults) 2 . The breeding range is divided latitudinally in the model from 30°-35°, 35.1°-40°, and 40.1°-50°f or the South, Central and North regions, respectively and longitudinally from 75°-105°( Figure 1). A system of coupled partial differential equations described briefly in the following sections, and detailed in Supplementary Appendix A1, generates the model's dynamics. Monarch development time is temperature-dependent and typically measured in Degree-Days (°D). As such, spatio-temporal variation in temperature and milkweed habitat drive Monarch population dynamics in the MOBU-SDyM through a temperature variable based on historical records. That variable is converted into°D to calculate the expected total development time specific to each region at each time step. The model transitions individuals from the Immature to the Adult class using the Monarch's requirement of°D, bracketed by the Monarch's heat impairment threshold and developmental zero temperatures (Zalucki, 1982). Additionally, the adult's requirement of°D (Zalucki, 1981) dictates the average lifespan of Adults. Temperature, along with soil moisture, also drives milkweed habitat in the model, which, in turn, is essential for the reproductive success of the Monarch (Veihmeyer and Hendrickson, 1927;Jeffery and Robison, 1971;Flanagan and Johnson, 2005;Pleasants and Oberhauser, 2012).
The MOBU-SDyM initializes in January 1994, with all individuals concentrated in the Overwintering class and remaining there until mid-March. The initial value model represents the 6.23 ha of overwintering colonies recorded in 1994 multiplied by a density conversion factor (OWdens), defined later in this document, to obtain the actual number of Monarchs. The model then has individuals transition to the Breeding classes based on the estimated monthly proportion of individuals moving across their range along with their associated mortality (Flockhart et al., 2015). Once within the Breeding classes, a transition matrix (Flockhart et al., 2015) dictates individuals' transitioning across the three Breeding regions (South, Central, and North). While within the Breeding regions, Adults reproduce for approximately four generations, considerably increasing the size of the population.
A breeding season's success is dictated in the model mainly by the number of individuals arriving from the Overwintering class, density-dependent effects, temperature, and habitat availability. The south migration begins around mid-August with individuals transitioning from the North to Central regions and then mid-September to the South region. Finally, when the migration trigger is set off in the South region (around early October), individuals take an average of 15 simulation days to transition to Simulated Monarchs reproduce and disperse across the three breeding regions during the breeding months and mobilize to a non-breeding overwintering region during the overwintering season. The colored dots give the key to the system drivers that affect MOBU-SDyM dynamics. Numbers in parentheses provide sections in the paper where a description of that driver can be found.
the Overwintering class. Individuals that return to the Overwintering class do not reproduce and reduce in numbers due to background and weather-related mortality. By mid-December, all the individuals within the simulation concentrate in the Overwintering class, remaining there until mid-March to start a new cycle.

Population Dynamics
Here we detail how the model captures the Monarch's lifecycle dynamics by presenting the stock equations that define the number of Immatures and Adults across the different sections of the model. These descriptions are further explained as the flows that integrate across the stock equations 3 . The life cycle dynamics outlined here are repeated for each of the Breeding regions, with any differences in the Adult classes of different regions highlighted. The remaining section, Drivers, focuses on defining the factors that drive the life cycle dynamics and how they are explicitly incorporated into the model.

Immature Individuals
The number of Immature individuals (including eggs, larvae, and pupae) found in region x at the end of time step 4 t (Immature x, t ; Eq. 1) depends on four variables: the number of Immature individuals at the end of the previous time step (Immature x, t−1 ), the number of eggs laid in the current time step (E x,t ) and the Immature numbers transitioning out of the class because of death( M x,t ) or reaching maturity in the current time step (Id x,t ).
From Eq. 1, the number of eggs laid per Adult female multiplied by the total number of females in time t determines the number of new individuals joining the Immature class (E x,t ; Eq. 2). The total number of Adults (Adults x,t ) and the sex ratio (Sr t ) dictates the number of females in the Adult class. In turn, the number of eggs laid per female in a time step is dictated by (i) the average oviposition potential per female (Ef x,t ), which is dictated by the estimated mean female age; and (ii) environmental conditions such as the interpatch distance (Dis x,t ), depensatory effects (Al x,t ), and whether the time step is before or after the fall migratory trigger is set off (determined by the Boolean Ft x,t ).
Previous Monarch models have assumed a constant development rate for the Monarchs (Flockhart et al., 2015;Oberhauser et al., 2017). However, this does not capture development time's variation related to the number of Degree-Days (°D) to which individuals are exposed. The MOBU-SDyM captures this variation by calculating the number of individuals transitioning from the Immature to the Adult class at a time step (M x,t ; Eq. 3) as the number of Immature individuals (Immature x, t ) divided by the immature development time in region x, at time step t (see Section 2.2.4.1 below).
The number of new individuals that die daily within the Immature class (Id x, t ; Eq. 4) is a conditional equation that combines the mortality effect of freezing temperatures (-10°C) with background mortality. When the mean modeled temperature within a region is below -10°C, all the new Immature individuals in that region during that time step die. If the temperature is above that threshold, the sum of individuals joining the Immature class in that time step (E x, t ) multiplied by a basal death rate (IDd x, t ) and a density-dependent death rate (Densd x, t ) dictate the background mortality of Immature individuals. Since available literature only states overall survivorship of the immature individuals (rather than daily death rates), the model calculates the number of Immature individuals dying per time step (Id x, t ) as the number of eggs laid in that time step, i.e., the new individuals entering the class, multiplied by the proportion of all Immatures expected to die.

Adult Individuals
We estimated the number of individuals within the Adults classes (Adults x,t ; Eqs 5-7) for each region every time step using similar stock equations. The number of individuals within the Adults class increases with individuals transitioning from the same region's Immature class (M x,t ), and individuals transitioning from Adult classes of other regions (D x y,t where x is the current region, and y is any of the other breeding regions). Conversely, the number of individuals decreases with outgoing individuals toward Adult classes in other regions (D y x,t where x is the current region), the incoming individuals from other regions who die (Sd x,t ), the number of dying Adults within the region due to their age (Ad x,t ) and the number of Adults dying at the start of their fall migration (Fd y,t where y is any of the three breeding regions). The South region has an additional term that describes the movement of Adult individuals to the Overwintering region (Fs t ).
Adults North,t Adults North,t−1 + M North,t + D North Adults Central,t Adults Central,t−1 + M Central,t + D Central Adults South,t Adults South,t−1 + M South,t + D South Overwintering,t + D South The movement of any one individual across Breeding regions is assumed to be instantaneous (Eq. 8), and the model uses lookup tables to describe the monthly proportion of individuals moving across the four regions (D y x,t ). The Overwintering transition to the Breeding regions expands across 30 days before May, transitioning out any remaining individuals on May 1st (Eq. 9). Previous research shows that, for migrating Lepidoptera, mortality during non-migratory stages of their cycle is relatively low compared to migrating phases (Chapman et al., 2012). The MOBU-SDyM assumes that adults moving across regions (D y x,t ) have a higher probability of dying than the ones remaining within the same region (Dd y x, t ), which is reflected in the model by a specific death rate during movement (Sd y x,t ; Eq. 10).

Migration and Overwintering
A Sun angle of 52°at solar noon is assumed to be a cue for Monarch migratory behavior, although little peer-reviewed research exists in this regard (Perez and Taylor, 2004;Taylor et al., 2019). However, the migratory behavior does not appear precisely the same day every year as it would happen if the Sun angle were the only cue. Evidence suggests that temperature and host plant senescence additionally influence the reproductive diapause before migration (Goehring and Oberhauser, 2002), although the relationship between Sun angle, plant senescence, and weather variables is still unknown. The model explores this uncertainty by assuming a non-fixed migratory trigger (Ft x,t ), i.e., not dependent on a specific date. Instead, the combination of Sun angle at solar noon below 52°and a five-consecutive-day mean weather temperature below a posterior sample of a temperature migration threshold sets off the migration trigger in each model's Breeding region. As in the natural system, followed by migration trigger sets off earlier in the North region, following the Central and being last in the South. The priors used to estimate the temperature migration threshold's posterior distribution were the mean temperatures during the months that migration began within each region between 1994 and 2019.
Once migration is triggered within each region (Ft x,t 1), a southward transition of Adults starts, being affected by stochastic environmental factors that increase their death rate while migrating (FMigDR x,t ) and consequently, the number of dead individuals (Fd x,t ; Eq. 11) at every time step.
The Central region migrating Adults, already merged with individuals from the North, converge in the South until the migration trigger is activated there. Then, the South region transitions all the Adults toward the Overwintering region. Since the distance from the south to the overwintering region is considerably more than the distance among breeding regions, the MOBU-SDyM assumes that migration will take 15 days for an individual beginning its migration from the South region (Fs t ; Eq. 12) to its arrival at the Overwintering region (Fe t ; Eq. 13). Mid-Texas' distance to the overwintering grounds and the average flight speed of south-migrating Monarchs [65-80 km/day; (Downhower, 1988;Howard and Davis, 2015)] give the value to this time delay.
Fe t delay (Fs t , 15) − Fs t p Hd t 15 (13) To accommodate for the 15 days between departure from the south and arrival to the overwintering grounds, along with the mortality associated with it, we added a stock variable (Migrating t ; Eq. 14) holding migrating individuals until their arrival to the Overwintering region.
During that migration time, many factors not sufficiently studied but broadly assumed to be present (Hd t ) affect the Monarch's survival: road mortality, mosquito-controlling pesticides (Tracy et al., 2019), "sequestration of migrants" by tropical milkweed (Asclepias curassavica) in the southern region (Batalden and Oberhauser, 2015), and climatic adversities such Atlantic hurricanes (Ries et al., 2018), and droughts that could strongly alter the availability of nectar sources (L. P. Brower et al., 2006). The MOBU-SDyM considers these factors by including an extra term to the mortality definition (Fwd t ; Eq. 15).

Fwd t Migrating t p Hd t
The Overwintering region holds all the individuals within the system for approximately five months (Adults Overwintering,t ; Eq. 16); the exact duration depends on when the South region's migration trigger is activated. The Overwintering region increases its number of Monarchs with the arrival of individuals from the Breeding regions (Fe t ), and decreases its number through the season (Od t ; Eq. 17) due to predation and background mortality (Pd t ), and weather-related factors (Cd t ). When the overwintering season ends, the surviving Adults transition back north to the Breeding regions (D y Overwintering,t where y is any of the three breeding regions).  (Zalucki, 1982); in theory, at low or high temperatures, a Monarch could take as many as 30 days or as few as ten to reach maturity respectively. Moreover, research on Monarchs (Zalucki, 1981) and other butterfly species (Gossard and Jones, 1977) suggests that the adult lifespan is also negatively related to temperature. The MOBU-SDyM captures that variability in development time (Adult lifespan x,t ; Eq. 18 and Immature Development time x,t ; Eq. 19), number of deaths (Ad x,t ; Eq. 20), and number of Immature individuals maturing per time step (M x,t ; Eq. 21) by calculating the daily Degree-Days in a region (°D x,t ; Eq. 22) (Baskerville and Emin, 1969) using the current simulated temperature 5 (Temp x,t ), the Monarch's developmental zero 6 (Tl 11.5°C) and heat impairment threshold 7 (Tth 32°C). To calculate the proportion of individuals transitioning through classes, the Adults°D requirement to reach their estimated longevity (Ta x,t ) or the sum of°D required for immature Monarchs to overcome their specific developmental stage (Th x,t , Te x,t , and Tp x,t for an eggs, a larvae, and pupae, respectively) are divided by the current°D to give a current expected Adult lifespan (Eq. 18) or current expected Immature development time (Eq. 19), respectively.
Adult lifespan x,t Ta x,t°D DD°x ,t if Temp x,t > Tth : OR : Temp x,t < Tl then (0) else Temp x,t − Tl

Milkweed Availability
In addition to its effect on Monarch development time, milkweed growth and subsequent habitat availability are also affected by temperature. Most experts agree that habitat availability across the breeding regions is the primary driver of Monarch population size (Pleasants and Oberhauser, 2012;Pleasants, 2017;Thogmartin et al., 2017a;Stenoien et al., 2018). Previous studies used the estimated total number of milkweed stems that a region can support as a proxy for habitat availability (Flockhart et al., 2015;Thogmartin et al., 2017b). However, it is essential to consider that the actual number of stems available to the Monarch is highly dependent on weather conditions, e.g., a drought can reduce the number of stems across the landscape. Considering that butterflies are an R-selected species (Pianka, 1970) with highly fluctuating population sizes, models should capture the factors behind these fluctuations, such as dynamic milkweed availability across the landscape. However, no peerreviewed studies exist linking weather and milkweed growth to the necessary level of detail.
With no quantitative studies linking abiotic parameters to milkweed growth, the MOBU-SDyM includes a hypothesized relationship between temperature, soil moisture, and milkweed availability (Milkweed Availability x,t ; Eq. 23). Since those parameter values are unknown, we used a Bayesian parameter inference approach, using as prior early work on soil moisture and temperature for other similar plant species (Veihmeyer and Hendrickson, 1927). The interaction between the hypothesized ideal temperature (T), soil moisture ( H) and a scaling moisture estimate (Me) for milkweed growth defines such a curve. When the ideal temperature and soil moisture are met within a region, the milkweed availability is at its highest; if any of these parameters are below or above this value point, the availability decreases. This relationship is instantaneous, meaning there is no time between conditions and the subsequent milkweed availability. Further iterations of this model might try to implement such an element. Supplementary Appendix A1 explains and justifies the form of Eq. 23.

Interpatch Distance
Recently, the discussion of milkweed as the main driver of Monarch decline has expanded from considering absolute stem numbers to including their distribution across the landscape (Cutting and Tallamy, 2015;Kasten et al., 2016;Zalucki et al., 2016;Stenoien et al., 2018). We decided to account for interpatch distance since it is an element with considerable impact on the Monarch's egg-laying success (Zalucki et al., 2016). We included this effect by estimating the average interpatch distance across the landscape (Dis x, t ) with a logistic-shaped relation (Eq. 24) shaped by a midpoint (Ia) and slope (Ib) based on Zalucki et al. (2016) and experts' opinion of how egg-laying success decreases with distance needed to fly between patches.
Egg Laying Success x,t 1 e Ip x,t Ia Ib The MOBU-SDyM assumes that each patch is circular and estimates the interpatch distance (Dis x,t ; Eq. 25) of the average nearest neighbor (Gotoh et al., 1978) as the square root of the available milkweed patches on the current region (x) on timestep t, divided by the area within the region that is not a suitable habitat for a Monarch to lay its eggs (Aw x ; Eq. 26). The habitat not suitable for Monarch oviposition was considered to be the 5 Randomly generated using the temperature distribution backcasting provided by (Flato, 2007) for every region. 6 Developmental Zero is the temperature below no measurable development occurs (Zalucki, 1982). 7 Heat Impairment Threshold is the temperature above measurable development does not occur (Zalucki, 1982).
Frontiers in Environmental Science | www.frontiersin.org May 2021 | Volume 9 | Article 634096 total area of the region minus the overall area of available milkweed patches. The number of available milkweed patches (P x,t ; Eq. 27) depends on the available milkweed area (Am x,t ; Eq. 28) divided by the mean patch area [Pa x,t ; (Hartzler and Buhler, 2000)]. The available milkweed area (Am x ), in turn, is the result of dividing the number of available stems [MWa x,t ; Eq. 29; (Flockhart et al., 2015)] by the stem density [MWdens; (Ralph, 1977)]. Based on the authors' empirical observations, the MOBU-SDyM assumes that, a female will lay its eggs only if no other adult individuals are on that same stem. So, the available stems result from subtracting the number of Adults within a region from the total number of stems in that same region. Finally, the area not suitable for egg-laying is defined as the total area of region x (A x ) minus the area occupied by milkweed.

Egg Density
According to Flockhart et al. (2012), there is a marked decrease in egg survival when the egg density on a milkweed stem exceeds a threshold. Considering that the milkweed is, arguably, the most critical limiting factor for Monarch success, the densitydependent death rate (Densd x,t ; Eq. 30) of eggs is a likely key variable to consider. This parameter is estimated based on a laboratory experiment from Flockhart et al. (2012) and the egg density is calculated by dividing the number of eggs (Eggs x, t ; Eq. 31) by the total available stems in region x during timestep t (Milkweed Availability x,t ). Since the model does not consider a separate class for eggs (the Immature class includes eggs plus all larval states plus pupae), the number of eggs is back-estimated by dividing the current days to hatch (Th x, t ) by the total number of days to complete the Immature stage: egg (Th x, t ), larva (Tm x, t ), and pupae (Tp x, t ), and multiplying this by the total number of Immature individuals (Immature x, t ).

Adult Density
Problems that arise at small population size, e.g., inbreeding depression, difficulty in finding a mate, and protection against predators, are grouped within the broad category of depensatory or negative density-dependent effects (also known as Allee effects), all of which can hinder a population's potential to recover from a perturbation (Freckleton et al., 2005). Currently, no studies have linked the Monarch directly with any of these phenomena. However, these effects are present in other butterfly species (Murphy et al., 1990;Kuussaari et al., 1998). Here, the MOBU-SDyM assumes an individual's reproductive success logistically decreases once it reaches the lower threshold of approximately 0.5 individuals per hectare (Eq. 32). We employed a Powell hill-climbing optimization algorithm to estimate this parameter (Powell, 1971;Burns and Janamanchi, 2007), using as a start search parameter the estimated Monarch density across the breeding regions during the lowest overwintering season recorded (2014)

Fall Migration Mortality
Multiple demographic drivers exist during the fall migration that could strongly alter the mortality during this stage, e.g., availability of nectar sources (L. P. Brower et al., 2006), road mortality, mosquito-controlling pesticides (Tracy et al., 2019), "sequestration of migrants" by tropical milkweed (Asclepias curassavica) in the southern region (Batalden and Oberhauser, 2015), and climatic adversities such as Atlantic hurricanes (Ries et al., 2018) and droughts. Moreover, some authors propose that some of these drivers may not be solely detrimental but could improve migration success (Ries et al., 2018;Saunders et al., 2019); however, there is an ongoing debate about their magnitude and impact. The MOBU-SDyM indirectly incorporates these arguments by including a fall migration mortality parameter (named "Hurricanes Estimate" in the model) parameterized with its estimated posterior distribution. The MOBU-SDyM assumes that all the factors comprising the fall migration mortality are constant over time, except for the hurricane's season intensity, expected to increase in the future (Villarini and Vecchi, 2013). The model uses the historical records of the Power Dissipation Index 8 and its forecasted behavior (Villarini and Vecchi, 2013) as a proxy of mean hurricane intensity. The posterior distribution of the fall migration death rate gives a general sense of the possible magnitude of mortality during the fall migration but does not discern which of all those elements are the main contributors to that value. Eq. 13 in Section 2.2.3 section above describes the way the "Hurricanes estimate" (Hd t ) is incorporated into the model.

Weather-Related Overwintering Mortality
Most of the eastern Monarch population gather in just a few forest patches from early November to mid-March (Urquhart and Urquhart, 1978). This elevated concentration of individuals makes the overwintering sites extremely vulnerable to perturbations, where a localized threat can be devastating . Although predation by birds and rodents is a constant pressure, the greatest threat to the overwintering Monarchs is extreme weather events such as winter storms (Calvert et al., 1983). Some such events have killed up to 30% of the entire population at a time (L. P. Brower et al., 2004). Exposure to open air due to forest canopy degradation, driven mainly by illegal logging, can exacerbate such casualties (Anderson and Brower, 1996;Saunders et al., 2019).
In the model, the number of individuals dying within the Overwintering region (Od t ; Eq. 33) is a function of the proportion of individuals who die from background drivers, e.g., predation, tourism, exhaustion (Pd t ), and the proportion of individuals who die due to extreme weather events (Cd t ; Eq. 34). The MOBU-SDyM calculates the extreme weather-related death rate based on the interaction of precipitation (Ppt OW,t ) and temperature (Temp OW,t ), the probability of Monarchs' mortality to freezing temperatures when they are either wet (OWd wet ; Eq. 35) or dry (OWd dry ; Eq. 36), and the canopy's "blanket effect" on the Monarch's body temperature (BodyTemp t ; Eq. 37), which is the function of the parameter Exposure to open sky (Anderson and Brower, 1996). Since the exposure to the open sky is highly variable depending on every overwintering colony's specific location, which varies considerably every year, we estimated the posterior distribution for this parameter.

Od t
Overwintering t p Pd t + Overwintering t p Cd t (33)

Model Performance
We assessed the MOBU-SDyM's performance through posterior predictive checks (Meng, 1994) . We compared the standard deviates' means from both datasets using a One-Sample t-test and their variances with a chi-squared test. Complete cycle migratory models of the Monarch have used this same measure of performance in the past (Flockhart et al., 2015). A model's accuracy depends on how it can recreate the same mean and variance of the observed data and recreate the same type of oscillations. As a final test of model performance, we used a Receiver Operating Characteristic (ROC) curve (Hanley and McNeil, 1982) to evaluate the model's capacity to recreate the same oscillation pattern as the pattern of the observed data. The MOBU-SDyM performance metrics were estimated using R (R Core Team, 2019); the "bayesplot," "stats," and "pROC" packages were used to estimate the posterior predictive checks, standard deviates metrics, and ROC curve, respectively.

Model Calibration
There are not enough empirical data in most ecological systems for full parameterization without using highly-aggregated variables that ignore their underlying dynamics and impede emergent behavior within the system (Yates, 2012). Based on sound ecological theory, the structural hypothesis that guided our system suggested that seven parameters necessary to the model were either highly uncertain or had inexistent empirical data. The parameters whose posterior distributions were sampled to calibrate the model belong to one of three groups: 1) Milkweed growth: ideal temperature and humidity for milkweed growth and a moisture estimate scaling factor, 2) Fall migration: Temperature migration threshold and hurricanes estimate, and 3) Overwintering: Monarchs density, and exposure to open sky while overwintering.
We used a Bayesian inference approach (Tulsyan et al., 2016) through the VENSIM's built-in Markov Chain Monte Carlo & Simulated Annealing method to sample the posterior distributions of these parameters and estimate the most likely vector of values for all of them (Rahmandad et al., 2015). Supplementary Appendix A3 displays the convergence metrics. We also used the sampling of the parameters' posterior distributions to set confidence bounds around the model's output and establish confidence limits to our policy recommendations (Gelman et al., 2013).

Region's Leverage and Total vs. Available Milkweed Stems
One of our research objectives was to examine how varying the number of milkweed stems affects the overwintering colonies size. Furthermore, we wanted to test the difference in the resulting area of the overwintering colonies when the milkweed restoration strategy's objective assumes, or not, that every stem planted will be available and used by Monarchs, as currently is accepted. We accomplished these two objectives by setting up an experimental design to "artificially" increase and decrease the total number of stems by 3 × 10 6 , 141 × 10 6 , 723 × 10 6 , 2.06 × 10 9 , and 4.5 × 10 9 and then allocating them either to one region, equally divided over two regions only, or equally divided across the three regions. This design, consisting of 700 individual runs, ran once per 100 vectors of parameter values drawn from the posterior distribution previously obtained (7,000 runs in total). Finally, we ran the same experiment but increasing or decreasing the available stems instead of total stems, essentially bypassing any assumption of the resource's incomplete use. For this set of experiments, we used a stable version of the model with all the dynamic variables affecting the model over time fixed to their 2019 values, e.g., temperature, precipitation, sex ratio, PDI. For the analysis of these experiments, we performed a simple linear regression in R (R Core Team, 2019) of the mean simulated overwintering area, dependent variable, versus the number of milkweeds supplemented to each region (South, Central, and North), and their interactions. We averaged the best and most parsimonious candidate models and used them to explain the differences between using total and available stems.

Milkweed Improvement
Finally, we wanted to calculate the minimum yearly increase of total milkweed stems needed across the three predefined regions to produce a long-term mean size of the Overwintering colonies of 6 ha 9 . To meet this goal, we set an optimization objective that would reduce the residuals between the expected long-term mean area of 6 ha until the year 2050 and the simulated data by modifying the yearly increase in the number of total milkweed stems across the three breeding regions, while also penalizing the total number of stems needed to achieve that goal. We used a Powell optimization algorithm (hill-climbing optimization), which is the built optimizer in VENSIM. To reduce the computational burden, we set the optimization step to 1 million stems, i.e., the optimizer modified the search values in increments of 1 million. We ran four optimizations for each of a 10% random draw taken from the posterior distribution sample (∼5.3*10 6 simulations) that we estimated earlier and used the best payoff as the "optimal" solution, given the parameters. Since the actual policy objective set by policymakers is the total number of stems, this was the number we used as the independent variable in this analysis while, at the same time, we recorded the resulting available stems on each run.
We obtained and used climate projections from Villarini and Vecchi (2013) and the B1 scenario from the SRES experiment done by the Canadian Center for Climate Modeling Analysis model with T47 resolution (Flato, 2007) to parameterize the weather variables (PDI, soil moisture, and mean temperature) during the forecasting period to ensure our analysis was relevant under climate change conditions.

RESULTS
The results section first presents the obtained posterior distributions of the parameters described in Section 2.3 for the calibration of the MOBU-SDyM, then transitions to describe the model's performance, followed by an examination between total and available milkweed stems and the different leverage that the three Breeding regions have to changes in milkweed on the Overwintering colonies' modeled size. Lastly, we provide results on the milkweed stems' number needed to reach restoration goals.

Model Performance
The model showed a remarkably good performance in estimating the observed overwintering population's size, trend, and oscillations under the metrics we used. The posterior predictive checks of the models' data, also known as "Bayesian p-value" (Meng, 1994), showed 96% accuracy. In other words, the size of the overwintering colonies posterior distribution generated by the model included 96% of the times the data points from the observed colonies' size. The empirical cumulative distribution function's visual overlay, another step of the posterior predictive checks, also accurately resembled the observed data's cumulative distribution function (Supplementary Appendix A1). The ROC curve evaluated the model's accuracy to predict the direction of observed data's year-to-year growth, i.e., if the population grows or shrinks. This metric yielded a level of accuracy of 67.91% ). The One-Sample t-test to the standard deviates mean, which tested differences between the observed and simulated growth rates, failed to reject the null hypothesis that it was different from 0 (95% CI: [-3.87, 0.54], t -1.56). The standard deviates variance, evaluated via a Chi-squared test, yielded a variance of 0.1215, rejecting the null hypothesis that it was equal to one, suggesting that the patterns obtained from the MOBU-SDyM tended to vary more than the observed data. The system's variance has been a problematic element to capture, and previous models have not been able to do so, either (Flockhart et al., 2015).

Model Calibration
The posterior distribution of all seven parameters subject to Bayesian inference, i.e., ideal temperature and humidity for milkweed growth, moisture estimate scaling factor, temperature migration threshold, hurricanes estimate, Monarchs density, and exposure to open sky while overwintering, fell within credible ranges (Figure 2). The joint posterior distribution for the first group of three variables representing milkweed growth reveals that, even at optimal conditions, stem availability is not greater than 31.37% (Figure 3). The parameter posteriors belonging to the fall migration group showed a considerably narrow distribution, which denotes the system's high sensitivity to temperature changes during the fall, migration delays, and mortality events during the migration. Finally, with the overwintering group's parameters, the exposure parameter's posterior was considerably skewed toward zero while the posterior for the density per hectare of overwintering Monarchs had a wide distribution ranging from 20 to 140 million Monarchs/hectare.

Region's Leverage and Total Vs. Available Milkweed Stems
We compared the differences between modifying the number of total and available stems, i.e., between assuming a perfect use of all milkweed or not due to weather conditions and Monarch's idiosyncrasy. Under the conditions of the model, an increase of available stems affected the size of the overwintering colonies, almost one order of magnitude more than the same addition of total stems. We also tested the relative contribution (leverage) to the Overwintering colonies' size of changes of milkweed stems on each of the three Breeding regions. The top regression models concur that one stem of milkweed in the Central region translates into a greater number of Monarchs for the next overwintering season and more than one stem in South, and considerably more than one stem the North region (Table 2). However, each region's relative contribution substantially changed if the measure used was either total or available stems. More specifically, while adding total milkweed stems in the Central region showed an effect 1.61 times larger than adding this same number of stems in the South, the difference dampened to 0.96 times more if considering available stems instead (Figure 4). This dampening effect is most likely due to climate change affecting milkweed availability and consequent Monarch productivity more in the Central region than in the South (Supplementary Appendix A2).

Milkweed Improvement
We explored different configurations of milkweed enhancement across the landscape in a search for its optimal allocation. In agreement with the previous step of the analysis, the payoff surface's slope was considerably steeper for the Central breeding region, followed by South, and the shallowest for North region. We did not find any evidence of multiple optima. Under this optimization, the MOBU-SDyM generated posterior distributions for the best habitat rehabilitation strategy entailing mean yearly increases of 1.86*10 8 , 2.41*10 9 , and 5.49*10 6 total stems for the South, Central and North regions, respectively. These numbers translate into 3.09*10 7 , 2.88*10 8 , and 6.36*10 5 available stems. Figure 5 depicts the population's long-term modeled trajectory under this optimization, and Figure 6, the optimal distribution of total and available number of stems given the posterior samples.

DISCUSSION
We developed the most comprehensive full life cycle model of the Monarch Butterfly biology to date, the MOBU-SDyM, which, despite its complexity, showed remarkable performance under several metrics. By including complex dynamic variables, along with the systems dynamics' analytical and computational strengths, we created new avenues to investigate the complex biological system of the Monarch's migration. Our model allowed us to test our original hypotheses, i.e., the South region's importance, the relevance of considering incomplete use of milkweed in setting policy objectives, and the consequent underestimation of the number of currently recommended milkweed stems needed to bring back the Monarch population to a minimum self-sustaining size. In the following section, we discuss our results, the relevance and limitations of including novel dynamic variables, and the implications of our contribution to the overall Monarch's conservation efforts moving forward.
We aimed to identify the Breeding region with the highest leverage (Senge and Forrester, 1980) on the following year's Overwintering colonies' size. Our sensitivity analysis showed substantial leverage of the milkweed availability across the Central region, slightly lower in the South region and almost negligible leverage in the North. Our results agree with previous researchers that identified the United States's Midwest region (represented by the Central region in our model) as the restoration priority for Monarch conservation (Thogmartin et al., 2017a). This difference in leverage between the South and Central region refutes our initial hypothesis that the South should have more leverage than the Central Region, even when such leverage is substantially dampened when considering available instead of total milkweed stems. One possible explanation behind this result is that high temperatures surpassing the heat impairment threshold of 32°C in the South during several months in mid-summer halt Immature development and considerably decrease the available milkweed stems. If this presumption is correct, the high summer temperatures would have generated a yearly bimodal peak of butterflies in the South, which, indeed, happens in the wild (Inamine et al., 2016) and also occurred in the model (Supplementary Appendix A2). However, this result does not mean that the Central region is the only critical one. We showed through our sensitivity analysis a non-linear interaction between the South and Central regions in which, regardless of considering total or available stems, the reduction (not increasing) of milkweed stems in both regions simultaneously brought about catastrophic consequences ( Figure 4). This result suggests that, even though the South may not play a key role in increasing the Overwintering colonies' size to the desired area, it seems vital to keep the population above a critically small number. Next, we estimated the mix of habitat restoration efforts across the three regions that would meet the safe minimum threshold of a longterm average of 6 ha of Overwintering colonies (Semmens et al., 2016). We showed that a significant focus should be habitat restoration in the Central breeding region, followed by the South and, to a lesser extent, in the North. The milkweed stem number recommended by our analysis brings the stem count considerably above the estimated number of stems across the breeding regions before 1995 (Flockhart et al., 2015), and it suggests a strategy with restoration efforts of one order of magnitude larger than previously determined (Pleasants, 2017). Our higher than previously suggested stem recommendation results from the imperfect use of the resource, meaning that not every single stem across the landscape will be visited by a Monarch throughout the season or be entirely usable when the Monarch visits. Interestingly, although our stems recommendation is considerably higher than previously suggested (Pleasants, 2017;Thogmartin et al., 2017b), when accounting only for available simulated stems over a 7-years simulation period 10 , the resulting number is appreciably close to the current recommendation but with a higher resolution in terms of those efforts being focused primarily on the Central region ( Figure 5). In other words, previous research is correct regarding the number of stems used by the Monarch but failed to make the distinction that such number is the number used by the Monarch and not the number required to be planted.
Combining our continent-wide milkweed allocation recommendations with land-type specific research can be a powerful guide for policymaking. For example, there is a call for using roadside and railroad right-of-ways for planting milkweed and other nectar flowers to support the Monarch's breeding and FIGURE 4 | Sensitivity analysis to changes in total and available stems across the three breeding regions-confidence intervals obtained from the posterior distribution of the uncertain parameters. Note the difference in magnitude (y-axis) between the effect of changing total (top row) and available (bottom row) stems.
10 Average simulation years it took the model to reach a stable minimum safe population.
Frontiers in Environmental Science | www.frontiersin.org May 2021 | Volume 9 | Article 634096 migration (Kasten et al., 2016). According to Pleasants (2017), there are 3.2 × 10 6 hectares of roadsides across the Midwest. To achieve the restoration requirement that we propose, using the milkweed stem density of 33,030 stems/ha (Ralph, 1977), the total roadside area required to be covered with milkweed would be 72,963 ha, i.e., 2.26% of roadsides across the Midwest. This is a feasible objective considering that roadsides are a relatively idle land where a mix of planting and strategic mowing (Knight et al., 2019) can be easily achieved, compared to other land types. However, Pitman et al. (2018) warns of the potentially increased mortality due to vehicle collisions. We recommend focusing more research on differential roadside mortality based on vehicle throughput and, probably, aiming roadside rehabilitation efforts on low throughput roads. Conservation strategies need to consider metapopulation dynamics such as patch configuration and heterogeneity while applying our recommendations too (Hanski and Ovaskainen, 2003;Haddad et al., 2017). Even though spatially explicit models would be more suitable for including this kind of element and that analyzing the effect of interpatch distance was not the model's primary objective, we decided to account for interpatch distance since it is an element with considerable impact over the Monarch's egg-laying success (Zalucki et al., 2016). Regarding metapopulation dynamics, we also included depensatory effects as a potential population driver; however, at least at the current population size, that effect proved negligible (Supplementary Appendix A1).
We added a dynamic temperature variable affecting the whole system by simulating its influence on Monarch biology via direct effects such as the heat impairment threshold, developmental zero, temperature migration threshold, and subsequent impacts on Monarch's developmental rate and longevity. Including this element becomes particularly relevant in the light of climate change, where other butterfly species' first and peak flight times have shifted by as much as ten days per 1°C climate warming (Roy and Sparks, 2000), leading to potentially critical phenology mismatches with their host plants (Posledovich et al., 2015). Our model produces such phenology shifts and shows that temperature affects Monarch population dynamics directly by delaying the migration onset and reducing the intergenerational time when temperatures rise (Supplementary Appendix A1). For example, constant warm temperatures could lead to one whole extra generation of butterflies within the yearly cycle (Supplementary Appendix A2), which, for a species with such reproductive potential as the Monarch, would entail population changes of multiple orders of magnitude (Altermatt, 2010). This result parallels actual biological events, where development time and Monarchs' longevity (and other insect species) strongly correlate with temperature (Gossard and Jones, 1977;Zalucki, 1981;Zalucki, 1981).
We also included indirect temperature effects as other variables' drivers interacting with the Monarch's biological FIGURE 5 | Expected recovery of Monarch overwintering counts following the amount and allocation of milkweed stems across the three Breeding regions prescribed by the optimization results. Each blue line represents one simulation using the specific optimum obtained given a specific parameter vector drawn from the posterior distribution. Orange bars represent the observed size of the overwintering colonies.
6 Developmental Zero is the temperature below no measurable development occurs (Zalucki, 1982). system, namely its interaction with the soil layer moisture content and its impact on milkweed availability. Common knowledge assumes that plants strongly depend on temperature and soil moisture, and a deviation from an ideal mix of these parameters could have consequences on the quality and availability of milkweed at a landscape level. However, despite soil moisture and temperature's importance on milkweed growth and, consequently, Monarch's population success, literature documenting such dynamics is scarce. We generated a joint posterior distribution of this relation and its effect on the "ideal" milkweed availability revealing that optimal milkweed availability peaks with a temperature of 28.1°C and 0.65 kg/m 2 of soil moisture (0-10 cm). However, even at this optimum, milkweed availability only reaches 31.37% because, being nongeographically explicit, the model assumes that an individual can utilize every stem available in the system, which is not the case. The soil moisture-temperature joint posterior captures the resource's incomplete usage by making only a portion of the total number of stems available at optimal conditions. Indeed, data obtained from the Monarch Larvae Monitoring Program showed a highly variable weekly density of 0.091 immature individuals/stem during May-September of 2010-2018 (Monarch Larva Monitoring Project, 2019). This imperfect total availability highlights the importance of differentiating between total and available milkweed stems. Our model, as others (Flockhart et al., 2015;Oberhauser et al., 2017), considers only the availability of common milkweed (Asclepias syriaca), basing all our calculations regarding milkweed on rough approximations of the available number of stems in the ground (Flockhart et al., 2015) and leaving aside several other Asclepias species that Monarchs also use (Malcolm and Brower, 1986).
Previous research in other species has highlighted the importance of considering habitat availability as a dynamic variable subject to temporally explicit weather variables. Peterson (1997) analyzed the relationship between host-plant phenology and dispersal of the Dotted Blue butterfly (Euphylotes enoptes) in which the egg-laying behavior was uneven and correlated with their host plant's artificially staggered phenology. Dennis et al. (2003) extended this idea by proposing a "Functional Resource-Based Concept" for habitat, drawing examples from several butterfly species. He argues that a patch's size and internal habitat dynamics, such as resource synchronization with stadia, resource disturbance, complementary resources required, and abundance frequency of resources, drive butterflies' meta-population dynamics.
We also accounted for weather-related mortality during the fall migration, another novel element of the MOBU-SDyM. The fall migration is, arguably, the Monarch's migratory cycle stage that carries most uncertainty about the nature and intensity of drivers. The migratory Monarch life cycle's importance has recently been heavily debated in the literature (Agrawal, 2019;Saunders et al., 2019). Unlike previous modeling efforts that, due to lack of data, opted to omit factors influencing the Monarch's fall migration mortality, we included a mortality variable grouping all of those variables (and other unknown ones) into a single Bayesian posterior distribution. The posterior FIGURE 6 | Required total number and allocation of milkweed restoration efforts given the posterior sampling as prescribed by the optimization results. The plot shows the difference between the numbers in the simulation as "prescribed by science," considering perfect use and availability of the resource (Red) and the simulated total number of stems a policy objective should consider (Cyan).
Frontiers in Environmental Science | www.frontiersin.org May 2021 | Volume 9 | Article 634096 distribution showed that its effect in the overwintering colonies area is historically small, though slightly increasing towards the present, and substantially increasing into the future (Supplementary Appendixes A1, A2) due to an expected Atlantic hurricane intensity increase (Villarini and Vecchi, 2013). The lack of reproductive activity during migration and overwintering explains this parameter's small effect over the overwintering colonies. More specifically, one fall-migrating Monarch's death represents the loss of only one overwintering Monarch; in contrast, one individual's death earlier in the season removes its expected potential demographic contribution to the next overwintering population. However, a factor not included in our model that could influence the effect of conditions during the fall migration is the loss of fat reserves during migration; current research is undergoing to bridge this knowledge gap.

CONCLUSION
The drivers behind the Monarch's decline have been a matter of controversy for at least 20 years. Initially, that search pointed to illegal logging at the overwintering sites as the main reason , later shifting to the well-supported milkweed limiting hypothesis (Pleasants and Oberhauser, 2012). Lately, alternative hypotheses have joined the discussion, such as the nectar-limiting (L. P. Brower et al., 2006) and the tropical milkweed-sequestration  hypotheses.
Our results bring further evidence to support the milkweedlimiting hypothesis as the Monarch's main population driver. Our results also acknowledge the importance of including alternative hypotheses and lesser drivers (e.g., depensatory effects, migration mortality, patch configuration, open sky exposure while overwintering) to provide more realistic modeling through a systemic view. Furthermore, our study advances the knowledge of the different breeding regions' relative contributions to the next overwintering season's success and, consequently, better informs the geographical allocation of milkweed in large-scale restoration efforts. We also highlight the importance of considering the incomplete use of milkweed in setting policy objectives, currently resulting in underestimating the milkweed stems needed to bring back the Monarch population to a minimum self-sustaining size.
Our results also underscore the need for a deeper understanding of fundamental variables to reduce uncertainty in our Monarch migratory population predictions. More specifically, there is a need for a better population metric at the overwintering sites and to expand, between Texas and the overwintering sites, the long-time standardized monitoring and reporting efforts found throughout the rest of the migratory flyway such as the Integrated Monarch Monitoring Program (Cariveau et al., 2019), Mission Monarch (mission-Monarch. org), The Monarch Larvae Monitoring Program (Kountoupes and Oberhauser, 2008), The TriNational Monarch Knowledge Network (birdscanada.org/birdmon/tmkn/main), and Monarch Watch (Taylor et al., 2019), to name a few. Also, in the light of climate change, modeling efforts not accounting for the dynamic and complex interdependencies between organisms and the changing climate are bound for failure. As such, research should focus intensely on getting empirical evidence to describe the soil moisture-temperature-milkweed growth system that we estimated here. Another important source of uncertainty that we found is the behavioral rules dictating how Monarchs interact with milkweed, i.e., all the cues at the micro and macro scale that determine individuals' movement and behavior. Ethological research must elicit those rules to design habitat restoration strategies that maximize Monarch's milkweed usage while reducing egg dumping and immature mortality.
Finally, the Monarch butterfly is a complex social-ecological system in which actions and perceptions from different stakeholders may play a determinant role in conservation strategies' support and success. Research is undergoing to couple the MOBU-SDyM and the lessons we learned here with previous social research ( Solis-Sosa et al., 2019) to give policymakers the necessary tools to conserve the Monarch Butterly's eastern migratory phenomenon for generations to come.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
RS conceived the original idea. RS designed and coded the model with input from CS. RS analyzed the model's data and implemented a Bayesian approach with guidance from SC. RS wrote the manuscript, and it was revised by AM, CS, SC, with contributions of ML. TRACE documentation was written by RS and edited by CS and AM.

FUNDING
The RenewZoo program provided funding (Grant No. 481954

ACKNOWLEDGMENTS
We want to thank the Montreal Insectarium for hosting RS during a considerable part of the modeling process giving him access to the incredible expertize amassed there. Special thanks to all the Monarch Butterfly experts who gave shape to this multi-year effort through their example or guidance.