Inappropriate human activities contribute to the degradation of ecosystems in arid or semi-arid regions. Therefore, emphasizing the importance of strategies for restoring vegetation in these areas cannot be overstated. However, there has been insufficient research on how to develop effective restoration strategies at minimal cost. This paper addresses this gap by studying how optimizing the spatiotemporal distribution of human activities through local and boundary controls can reduce the level of desertification in vegetation pattern structures, thereby facilitating the recovery of arid land vegetation. The results indicate that vegetation restoration depends on the proportion and number of human activity areas, with a trade-off between them. Furthermore, consistent conclusions were obtained on circular regions, demonstrating the robustness of the approach to boundary shapes. This paper aims to offer new insights into the restoration of arid land vegetation and the prevention of catastrophic ecosystem changes from the perspective of optimal control.

Topics
Control theory, Water model, Phase transitions, Desertification, Reaction-diffusion system, Population ecology

Research on strategies for restoring vegetation in arid regions has long been a focal point of attention. While much attention has been paid to the impact of human activities on the distribution of dryland vegetation, the focus has primarily been on statistical methods. This paper based on a dynamic model employs local and boundary controls to reveal the influence of the proportion and distribution of human activity areas on the structure of vegetation patterns. The results show that regardless of whether the areas are regular or irregular, as the proportion of human activity zones increases, the effectiveness of vegetation restoration initially rises before reaching a plateau. When the proportion of human activity areas remains constant, fragmented human activities are more conducive to transitioning vegetation toward pattern structures with lower levels of desertification. Moreover, there exists a mutual dependence between the proportion and distribution of human activity zones, with a discernible threshold delineating their relationship.

Desertification, a complex phenomenon resulting from the interplay of various factors, such as climate change and human activity, has become increasingly severe, particularly in recent years characterized by rapid economic development and a surge in extreme climate events.1 Actions, such as deforestation, overgrazing, and land development, have further exacerbated this issue. The Third Edition of the World Atlas of Desertification, published by the European Joint Research Centre, highlights the alarming extent of the problem, indicating that over 75% of the Earth’s land surface is currently experiencing degradation, with projections suggesting this figure may exceed 90% by 2050. The far-reaching impacts of desertification should not be underestimated for all nations worldwide. Most notably, desertification leads to biodiversity loss and exacerbates the emergence of extreme climatic conditions, culminating in a self-perpetuating cycle.2 Consequently, the imperative to address desertification has gained prominence as a pressing global concern.

Vegetation serves as a key ecosystem indicator, offering a direct means of assessing desertification levels within a region. Influenced by climatic conditions and human factors, vegetation exhibits a spatially heterogeneous yet somewhat structured distribution, encompassing diverse patterns. These patterns have been extensively documented through field observations and satellite imagery, particularly in arid and semi-arid regions, such as Niger,3,4 Israel,5 and Australia.6 Previous researchers have conducted a series of studies from various perspectives to unravel the intrinsic connection between these vegetation patterns and ecosystem desertification.7–12 These endeavors primarily involve the formulation of reaction–diffusion equations rooted in vegetation growth mechanisms, shedding light on the underlying mechanisms governing pattern formation. Among these studies, Klausmeier’s pioneering work introduced the vegetation–water reaction–diffusion equation, elucidating that these intricate vegetation patterns emerge as a result of Turing instability, a form of spatial self-organization.7 This seminal contribution has provided a solid foundation for a plethora of subsequent research.13–16 

Furthermore, related studies have unveiled the role of resource scarcity in arid landscapes, leading to phase transitions within vegetation patterns. These transitions encompass shifts between uniform vegetation and patterned states, transitions from one pattern state to another, and shifts from patterned states to desertification.17 Notably, these transitions are not merely gradual but can culminate in catastrophic shifts toward desertification.18 Therefore, the imperative of monitoring and preventing such catastrophic transformations is paramount.

The resilience of ecosystems is intimately connected to their capacity to withstand catastrophic transformations. Traditionally, resilience is defined as an ecosystem’s ability to maintain its core functions in the face of external disturbances.19,20 Holling quantified resilience by characterizing an ecosystem’s tolerance for changes;19 ecosystems can be deemed controllable when exposed to minor changes, but they may exhibit catastrophic shifts in response to substantial alterations.21,22 In spatially extended systems, such as drylands and savannas, critical transitions are intricately linked to the formation of self-organizing spatial vegetation patterns. These patterns serve dual functions, acting as early warning signals for critical transitions22,23 and as indicators of resilience, aiding in the avoidance of critical thresholds that might trigger catastrophic shifts.24–27 Remarkably, Rietkert et al. provide a comprehensive overview of practical strategies for ecosystems to circumvent catastrophic tipping points through a variety of spatial pattern formation processes.27 Although substantial efforts have been dedicated to researching ecosystem resilience, continuous and in-depth investigations remain imperative, given the intricate nature of ecosystems.

Additionally, the impact of human activities on ecosystems, recognized as one of the most significant factors, is multifaceted and dynamic. Human actions can either induce ecosystem collapse, as observed in cases of overgrazing and land reclamation or promote ecosystem restoration through afforestation and ecological engineering.28–30 Therefore, the strategic implementation of human interventions to avert catastrophic ecosystem changes is a question that warrants careful consideration (Fig. 1). Optimal control theory, a valuable tool for formulating optimal solutions for human activities, has found extensive application in fields, such as epidemiology and ecology, typically focusing on problems of resource and time optimization.31–34 Nevertheless, studies examining optimal control problems from the perspective of a pattern structure remain relatively scarce.35–39 Building on our previous work in optimal control, which primarily investigated the impact of heterogeneous human activities on vegetation pattern structures, we have chosen to refine our approach by considering localized optimization, encompassing both local regions and boundaries.40 This paper embarks on an exploration from the standpoint of local and boundary control, seeking to unravel the intricate influence of human activities on vegetation through theoretical analysis and numerical simulations, with the ultimate goal of providing novel insights into averting catastrophic shifts within ecosystems.

FIG. 1.
Drought and external disturbances lead to a catastrophic transition of vegetation from a pattern state to bare soil.21,23 This paper aims to reverse the pattern structure by regulating human activities to avoid desertification.
View largeDownload slide

Drought and external disturbances lead to a catastrophic transition of vegetation from a pattern state to bare soil.21,23 This paper aims to reverse the pattern structure by regulating human activities to avoid desertification.

FIG. 1.
Drought and external disturbances lead to a catastrophic transition of vegetation from a pattern state to bare soil.21,23 This paper aims to reverse the pattern structure by regulating human activities to avoid desertification.
View largeDownload slide

Drought and external disturbances lead to a catastrophic transition of vegetation from a pattern state to bare soil.21,23 This paper aims to reverse the pattern structure by regulating human activities to avoid desertification.

Close modal

The main framework of this paper is as follows. First, we present the dynamic model under study and propose an indicator to measure the trend of regional desertification levels with different pattern structures. Second, we introduce the optimal control problem and conduct a series of theoretical analyses, including the proof of the existence of the optimal solution and the first-order necessity condition. Finally, we validate the theoretical analysis through numerical simulations, offering an intuitive demonstration of how regulating the spatial distribution of human activities can help avoid catastrophic changes, thus facilitating the prevention and control of desertification. Additionally, we perform a robustness analysis to demonstrate the validity of our algorithm.

Water scarcity presents a pressing challenge for plant survival in semiarid or arid regions. The relationship between vegetation and water resources is complex and intricate. Von Hardenberg et al. developed a realistic model that captures the nonlinear interactions and feedback mechanisms between these two vital components of the ecosystem.8 Based on this classic model, this paper couples the local impacts of human activities,
(1)
where all variables are in a dimensionless form. The specific meanings of each parameter are presented in Table I. n and w represent vegetation biomass and soil–water density, respectively. First, introduce the meaning of each term in the first equation of Eq. (1), the term w n 2 represents the biomass resulting from vegetation converting water absorption into growth. The term − σ n accounts for the reduction in vegetation due to herbivory and mortality. The diffusion term Δ n represents the propagation of plants, which mainly includes clonal propagation and seed dispersal. The newly introduced term χ ω u characterizes the influence of human activities on vegetation biomass, encompassing actions, such as revegetation, woodcutting, and grazing, with ω ⊂ Ω denoting a region of human activities and χ ω serving as the characteristic function of ω⁠,
χ ω = { 1 , x ∈ ω , 0 , x ∉ ω .
In the second equation of Eq. (1), water primarily arises from rainfall, denoted as a⁠, and water consumption mainly comprises evaporation and absorption via vegetation growth, represented by − w and − w n 2⁠, respectively. The term μ ( Δ w − δ Δ n ) accounts for the self-diffusion of water and the impact of vegetation diffusion on water. Additional details regarding the model and its dimensional parameters can be found in the literature.41 
TABLE I.

The implications of the relevant parameters in Eq. (1).

Parameter Biological meaning Value Source
σ  Vegetation loss rate  1.8  41  
a  Precipitation  4.4  41  
μ  Soil–water diffusion coefficient  10  41  
δ  Soil–water diffusion feedback intensity  [0, 0.17]  41  
Parameter Biological meaning Value Source
σ  Vegetation loss rate  1.8  41  
a  Precipitation  4.4  41  
μ  Soil–water diffusion coefficient  10  41  
δ  Soil–water diffusion feedback intensity  [0, 0.17]  41  
View Large

Evidently, in the absence of human activity impact, this model aligns with the one investigated by Sun et al., whose research results indicate that the water absorption capacity of vegetation roots can induce pattern phase transition.41 We also briefly studied the spatial distribution of vegetation under different water absorption capacities and found that vegetation primarily exhibited uniform, white-eye, striped, and green-eye distributions across different interval segments of water absorption capacity. It is worth noting that our study differs from the research of Sun et al. in the absence of mixed patterns. This distinction may be attributed to the insufficient setting of the terminal time during the simulation process of Sun et al., preventing the pattern structure from reaching a steady state (also see Fig. 2).

FIG. 2.
Spatial distribution of vegetation under different cross-diffusion coefficients δ. As the soil–water diffusion strength increases, the vegetation pattern undergoes a phase transition and the structure type of the soil–water diffusion strength remains the same for a certain interval. The relevant parameter values can be found in Table I.
View largeDownload slide

Spatial distribution of vegetation under different cross-diffusion coefficients δ⁠. As the soil–water diffusion strength increases, the vegetation pattern undergoes a phase transition and the structure type of the soil–water diffusion strength remains the same for a certain interval. The relevant parameter values can be found in Table I.

FIG. 2.
Spatial distribution of vegetation under different cross-diffusion coefficients δ. As the soil–water diffusion strength increases, the vegetation pattern undergoes a phase transition and the structure type of the soil–water diffusion strength remains the same for a certain interval. The relevant parameter values can be found in Table I.
View largeDownload slide

Spatial distribution of vegetation under different cross-diffusion coefficients δ⁠. As the soil–water diffusion strength increases, the vegetation pattern undergoes a phase transition and the structure type of the soil–water diffusion strength remains the same for a certain interval. The relevant parameter values can be found in Table I.

Close modal
Furthermore, to quantitatively assess potential desertification trends corresponding to different vegetation pattern structures, we have introduced the Regional Vegetation Desertification Level Index ( R V D L I )⁠,
R V D L I = 1 | Ω | ∫ n ( x , T ) ≤ n ∗ n ( x , T ) d x ,
where n ( x , T ) represents the steady-state vegetation pattern, and n ∗ is the vegetation biomass under a uniform distribution. We present the R V D L I under various vegetation pattern structures, and the results indicate significant disparities in R V D L I among different pattern structures, characterized by a marked increase in R V D L I when the pattern structure type shifts (Fig. 3, left). Among these, the R V D L I is the smallest for the homogeneous state, while the R V D L I for the white-eye pattern, the striped pattern, and the green-eye pattern progressively increase, signifying the lowest likelihood of desertification in the homogeneous vegetation state and the highest probability in the green-eye distribution pattern structure. Hence, the green-eye pattern can serve as an indicator for early desertification warning. Furthermore, the R V D L I for all pattern structures, except for the uniform vegetation state, exhibits a declining trend with increasing cross-diffusion intensity δ⁠. As δ represents the parameter that adjusts the effect of vegetation diffusion on water, this decline indicates that higher values of δ lead to a stronger impact of vegetation diffusion on water, causing greater water depletion, thus contributing to an increased severity in desertification trends (Fig. 3, right).
FIG. 3.
Left image: R V D L I corresponding to different types of pattern structure. The R V D L I jumps up in steps as the pattern type shifts sequentially between homogeneous, white-eye, stripe, and green-eye. Right image: partial enlarged view of R V D L I corresponding to each pattern type. R V D L I tends to increase with increasing the soil–water diffusion intensity δ. Other parameter values can be found in Table I.
View largeDownload slide

Left image: R V D L I corresponding to different types of pattern structure. The R V D L I jumps up in steps as the pattern type shifts sequentially between homogeneous, white-eye, stripe, and green-eye. Right image: partial enlarged view of R V D L I corresponding to each pattern type. R V D L I tends to increase with increasing the soil–water diffusion intensity δ⁠. Other parameter values can be found in Table I.

FIG. 3.
Left image: R V D L I corresponding to different types of pattern structure. The R V D L I jumps up in steps as the pattern type shifts sequentially between homogeneous, white-eye, stripe, and green-eye. Right image: partial enlarged view of R V D L I corresponding to each pattern type. R V D L I tends to increase with increasing the soil–water diffusion intensity δ. Other parameter values can be found in Table I.
View largeDownload slide

Left image: R V D L I corresponding to different types of pattern structure. The R V D L I jumps up in steps as the pattern type shifts sequentially between homogeneous, white-eye, stripe, and green-eye. Right image: partial enlarged view of R V D L I corresponding to each pattern type. R V D L I tends to increase with increasing the soil–water diffusion intensity δ⁠. Other parameter values can be found in Table I.

Close modal

The results of Sec. II showed that as the soil–water diffusion intensity δ becomes stronger, vegetation assumes a spotted distribution resembling a “green eye” pattern. Clearly, the distance between vegetation patches significantly expands, and the spatial distribution of vegetation takes on a more isolated character. This phenomenon could potentially lead to a drastic reduction in vegetation’s water storage capacity, ultimately culminating in desertification. To avert this scenario, we consider appropriate regulation of human activity, mainly including intra-regional and boundary regulation.

Intra-regional regulation considers practical considerations; accordingly, we opt to implement human activity restrictions within local areas. We aim to investigate how variations in the number and size of these human activity zones impact desertification prevention and control. On the other hand, boundary regulation addresses the rapid urban expansion and escalating land occupation rates. Given the complexities of real-world factors, we favor the use of the Robin boundary as the preferred boundary control method.

However, while striving to mitigate desertification through the regulation of human activities, we must also account for the associated costs. To address this concern, we formulate the following optimal control problem:
(2)
where κ i > 0 ( i = 1 , 2 , … , 5 ) represents weight coefficients; n t a r and w t a r denote the selected target patterns; and u⁠, b n⁠, and b w are used to describe the spatiotemporal distribution of human activities within the region and on the region boundary. They are subject to the constraints outlined in the feasible control set K f e⁠,
K f e = { ( u , b n , b w ) ∈ L ∞ ( Ω T ) × L ∞ ( Σ T ) × L ∞ ( Σ T ) | u ( x , t ) ∈ [ u 1 , u 2 ] f o r a . a . ( x , t ) ∈ Ω T , b n ( x , t ) ∈ [ b n 1 , b n 2 ] , b w ( x , t ) ∈ [ b w 1 , b w 1 ] f o r a . a . ( x , t ) ∈ Σ T = ∂ Ω × ( 0 , T ) } ,
where Ω ⊂ R 2 is an open bounded set, with boundary ∂ Ω ∈ C 2 + α ( α > 0 )⁠. n and w satisfy the following state system:
(3)
where g 1 ( n , w ) = w n 2 − σ n⁠, and g 2 ( n , w ) = a − w − w n 2⁠.

Based on this optimal control problem, we provide first-order and second-order optimality conditions (see  Appendix A) and then use the variable step projection gradient method (see  Appendix B) for numerical simulation to reveal the impact of human activities on vegetation patterns.

The proportion of human activity areas undeniably exerts a significant influence on the spatial distribution of vegetation. As depicted in Fig. 4, the figure illustrates the impact of varying the proportion of human activity areas on a vegetation pattern shift, with the white-eye pattern as the target configuration. The outcomes of this analysis reveal a gradual convergence between the spatial distribution of vegetation and the predefined target pattern as the proportion of human activity areas increases. At a human activity area proportion of 23.04 %⁠, the spatial distribution of vegetation predominantly comprises a mixture of green-eye and white-eye patterns, signifying suboptimal control effectiveness. However, as the size of human activity areas is expanded to 36 %⁠, the vegetation exhibits a white-eye pattern, aligning closely with the desired target configuration. Further augmentation of the human activity area, exceeding 50 %⁠, results in a spatial distribution of vegetation that perfectly conforms to the target pattern. This observation underscores the critical role of the proportion of human activity areas in achieving the desired vegetation distribution, with larger proportions contributing to a perfect transition between vegetation pattern structures.

FIG. 4.
The shift between pattern structures with increasing proportion of human activity areas under a constant number of human activity areas. The percentage at the top represents the corresponding proportion of human activity areas, increasing gradually from left to right. (a1)–(a4) The spatial distribution of vegetation. (b1)–(b4) The spatial distribution of human activity intensity within the area. Each small box represents a human activity area, and there are 25 in total. (c1)–(c4) The distribution of human activity intensity on the boundary (top, bottom, left, and right). Related parameter values: δ = 0.1, δ t a r g e t = 0.041, κ 1 = 1, κ 2 = 1, κ 3 = 1.0 e − 5, κ 4 = 1.0 e − 10, κ 5 = 1.0 e − 5, u 1 = b n 1 = − 2, u 2 = b n 2 = 5, b w 1 = b w 2 = 0, α 1 = 1, α 2 = 0, Ω = ( 0 , 100 ) 2, and other parameters are shown in Table I.
View largeDownload slide

The shift between pattern structures with increasing proportion of human activity areas under a constant number of human activity areas. The percentage at the top represents the corresponding proportion of human activity areas, increasing gradually from left to right. (a1)–(a4) The spatial distribution of vegetation. (b1)–(b4) The spatial distribution of human activity intensity within the area. Each small box represents a human activity area, and there are 25 in total. (c1)–(c4) The distribution of human activity intensity on the boundary (top, bottom, left, and right). Related parameter values: δ = 0.1⁠, δ t a r g e t = 0.041⁠, κ 1 = 1⁠, κ 2 = 1⁠, κ 3 = 1.0 e − 5⁠, κ 4 = 1.0 e − 10⁠, κ 5 = 1.0 e − 5⁠, u 1 = b n 1 = − 2⁠, u 2 = b n 2 = 5⁠, b w 1 = b w 2 = 0⁠, α 1 = 1⁠, α 2 = 0⁠, Ω = ( 0 , 100 ) 2⁠, and other parameters are shown in Table I.

FIG. 4.
The shift between pattern structures with increasing proportion of human activity areas under a constant number of human activity areas. The percentage at the top represents the corresponding proportion of human activity areas, increasing gradually from left to right. (a1)–(a4) The spatial distribution of vegetation. (b1)–(b4) The spatial distribution of human activity intensity within the area. Each small box represents a human activity area, and there are 25 in total. (c1)–(c4) The distribution of human activity intensity on the boundary (top, bottom, left, and right). Related parameter values: δ = 0.1, δ t a r g e t = 0.041, κ 1 = 1, κ 2 = 1, κ 3 = 1.0 e − 5, κ 4 = 1.0 e − 10, κ 5 = 1.0 e − 5, u 1 = b n 1 = − 2, u 2 = b n 2 = 5, b w 1 = b w 2 = 0, α 1 = 1, α 2 = 0, Ω = ( 0 , 100 ) 2, and other parameters are shown in Table I.
View largeDownload slide

The shift between pattern structures with increasing proportion of human activity areas under a constant number of human activity areas. The percentage at the top represents the corresponding proportion of human activity areas, increasing gradually from left to right. (a1)–(a4) The spatial distribution of vegetation. (b1)–(b4) The spatial distribution of human activity intensity within the area. Each small box represents a human activity area, and there are 25 in total. (c1)–(c4) The distribution of human activity intensity on the boundary (top, bottom, left, and right). Related parameter values: δ = 0.1⁠, δ t a r g e t = 0.041⁠, κ 1 = 1⁠, κ 2 = 1⁠, κ 3 = 1.0 e − 5⁠, κ 4 = 1.0 e − 10⁠, κ 5 = 1.0 e − 5⁠, u 1 = b n 1 = − 2⁠, u 2 = b n 2 = 5⁠, b w 1 = b w 2 = 0⁠, α 1 = 1⁠, α 2 = 0⁠, Ω = ( 0 , 100 ) 2⁠, and other parameters are shown in Table I.

Close modal

The spatial distribution of vegetation is also affected by the distribution of human activity areas, which can be either concentrated or dispersed. To reveal the influence of a human activity area distribution on vegetation patterns, we conducted an experiment where we adjusted the number of human activity areas while keeping the area shares constant. Figure 5 shows the transformation of vegetation toward the target pattern. The results indicate that vegetation patterns induced by human activities differ more from target patterns in concentrated human activity areas compared to relatively dispersed ones. For instance, when there is only one human activity area, the difference is more significant. The perfect transition of vegetation patterns has only been completed in local areas of human activity, while the spatial distribution of vegetation remains the same outside the activity area. In contrast, when human activity areas, such as those with 9 , 16 , 25⁠, or more, are relatively dispersed, the spatial distribution of vegetation is closer to the target pattern. In summary, our findings demonstrate that dispersed human activities are more conducive to the transition between vegetation patterns than concentrated human activities.

FIG. 5.
Conversion of vegetation patterns for different numbers of human activity areas when the proportion of human activity areas is kept constant ( 36 %). The numbers on the green arrow ( 1 , 9 , 16 , 25 ) indicate the number of human activity areas. a1)–(a4) The vegetation patterns. (b1)–(b4) The distribution of human activity intensity within the region. Each small box represents a human activity area, with the number of human activity areas from left to right being 1, 9, 16, and 25, respectively. (c1)–(c4) The distribution of human activity intensity on the boundaries of the region. δ t a r g e t = 0.088. Other parameter values can be found in Fig. 4 and Table I.
View largeDownload slide

Conversion of vegetation patterns for different numbers of human activity areas when the proportion of human activity areas is kept constant (⁠ 36 %⁠). The numbers on the green arrow ( 1 , 9 , 16 , 25 ) indicate the number of human activity areas. a1)–(a4) The vegetation patterns. (b1)–(b4) The distribution of human activity intensity within the region. Each small box represents a human activity area, with the number of human activity areas from left to right being 1, 9, 16, and 25, respectively. (c1)–(c4) The distribution of human activity intensity on the boundaries of the region. δ t a r g e t = 0.088⁠. Other parameter values can be found in Fig. 4 and Table I.

FIG. 5.
Conversion of vegetation patterns for different numbers of human activity areas when the proportion of human activity areas is kept constant ( 36 %). The numbers on the green arrow ( 1 , 9 , 16 , 25 ) indicate the number of human activity areas. a1)–(a4) The vegetation patterns. (b1)–(b4) The distribution of human activity intensity within the region. Each small box represents a human activity area, with the number of human activity areas from left to right being 1, 9, 16, and 25, respectively. (c1)–(c4) The distribution of human activity intensity on the boundaries of the region. δ t a r g e t = 0.088. Other parameter values can be found in Fig. 4 and Table I.
View largeDownload slide

Conversion of vegetation patterns for different numbers of human activity areas when the proportion of human activity areas is kept constant (⁠ 36 %⁠). The numbers on the green arrow ( 1 , 9 , 16 , 25 ) indicate the number of human activity areas. a1)–(a4) The vegetation patterns. (b1)–(b4) The distribution of human activity intensity within the region. Each small box represents a human activity area, with the number of human activity areas from left to right being 1, 9, 16, and 25, respectively. (c1)–(c4) The distribution of human activity intensity on the boundaries of the region. δ t a r g e t = 0.088⁠. Other parameter values can be found in Fig. 4 and Table I.

Close modal

The above-mentioned influence of the distribution of human activity areas (concentrated or dispersed) on the transformation between vegetation patterns is mainly analyzed visually from a structural perspective. To quantitatively characterize its impact on the transformation of vegetation patterns, we calculated the R V D L I corresponding to the transformed vegetation pattern (Fig. 6). The results indicate that when the proportion of human activity areas remains constant ( 36 % )⁠, the R V D L I of the corresponding vegetation patch structure decreases as the number of human activity areas increases. Therefore, dispersal is more favorable than concentrated human activity for the transition of vegetation to pattern structures with low levels of desertification. Moreover, our comparison with the initial pattern structure shows that optimizing human activity can effectively reduce the R V D L I of vegetation, regardless of whether human activity areas are dispersed or not. In summary, if the proportion of human activity areas is constant, then we can consider decentralizing the area, which can not only ensure the reduction of R V D L I of vegetation but also achieve better structural transformation, making it a more preferable choice.

FIG. 6.
R V D L I of vegetation patterns under the influence of different amounts of human activity areas. The numbers 0, 1, 9, 16, and 25 in the first row in the horizontal direction correspond to the number of human activity areas; the 0 % and 36 % in the second row represent the proportion of human activity areas. A bar with square markings indicates a transition to a white-eye pattern ( δ t a r g e t = 0.041 ), while a bar with circular markings indicates a transition to a stripe pattern ( δ t a r g e t = 0.088 ). The orange bar represents the R V D L I of the pattern structure generated without human activity influence. Other parameter values can be found in Fig. 4 and Table I.
View largeDownload slide

R V D L I of vegetation patterns under the influence of different amounts of human activity areas. The numbers 0⁠, 1⁠, 9⁠, 16⁠, and 25 in the first row in the horizontal direction correspond to the number of human activity areas; the 0 % and 36 % in the second row represent the proportion of human activity areas. A bar with square markings indicates a transition to a white-eye pattern ( δ t a r g e t = 0.041 )⁠, while a bar with circular markings indicates a transition to a stripe pattern ( δ t a r g e t = 0.088 )⁠. The orange bar represents the R V D L I of the pattern structure generated without human activity influence. Other parameter values can be found in Fig. 4 and Table I.

FIG. 6.
R V D L I of vegetation patterns under the influence of different amounts of human activity areas. The numbers 0, 1, 9, 16, and 25 in the first row in the horizontal direction correspond to the number of human activity areas; the 0 % and 36 % in the second row represent the proportion of human activity areas. A bar with square markings indicates a transition to a white-eye pattern ( δ t a r g e t = 0.041 ), while a bar with circular markings indicates a transition to a stripe pattern ( δ t a r g e t = 0.088 ). The orange bar represents the R V D L I of the pattern structure generated without human activity influence. Other parameter values can be found in Fig. 4 and Table I.
View largeDownload slide

R V D L I of vegetation patterns under the influence of different amounts of human activity areas. The numbers 0⁠, 1⁠, 9⁠, 16⁠, and 25 in the first row in the horizontal direction correspond to the number of human activity areas; the 0 % and 36 % in the second row represent the proportion of human activity areas. A bar with square markings indicates a transition to a white-eye pattern ( δ t a r g e t = 0.041 )⁠, while a bar with circular markings indicates a transition to a stripe pattern ( δ t a r g e t = 0.088 )⁠. The orange bar represents the R V D L I of the pattern structure generated without human activity influence. Other parameter values can be found in Fig. 4 and Table I.

Close modal
In order to quantitatively characterize the effectiveness of transitions between different patterns under the influence of human activity, the relative error between patterns before and after the transition is calculated as
R e l E r r o r = ‖ n ∗ − n t a r ‖ ‖ n t a r ‖ .
We have conducted a comprehensive analysis of the variation in the relative error concerning different proportions and numbers of human activity areas, illustrating this phenomenon through the transition to gap pattern (⁠ δ t a r g e t = 0.041⁠) as exemplars (see Fig. 7). The outcomes reveal a noteworthy trend: the relative error between the initial pattern and the target pattern diminishes as the proportion of human activity areas increases, while keeping the number of these areas constant. This suggests that augmenting the proportion of human activity areas proves advantageous for facilitating pattern conversion.
FIG. 7.
Contour plots of relative errors under different numbers and proportions of human activity areas. The horizontal axis represents the proportion of human activity areas, the vertical axis represents the number of human activity areas, and the colorbar represents the relative error. Here, δ t a r g e t = 0.041 and other parameter values can be found in Fig. 4 and Table I.
View largeDownload slide

Contour plots of relative errors under different numbers and proportions of human activity areas. The horizontal axis represents the proportion of human activity areas, the vertical axis represents the number of human activity areas, and the colorbar represents the relative error. Here, δ t a r g e t = 0.041 and other parameter values can be found in Fig. 4 and Table I.

FIG. 7.
Contour plots of relative errors under different numbers and proportions of human activity areas. The horizontal axis represents the proportion of human activity areas, the vertical axis represents the number of human activity areas, and the colorbar represents the relative error. Here, δ t a r g e t = 0.041 and other parameter values can be found in Fig. 4 and Table I.
View largeDownload slide

Contour plots of relative errors under different numbers and proportions of human activity areas. The horizontal axis represents the proportion of human activity areas, the vertical axis represents the number of human activity areas, and the colorbar represents the relative error. Here, δ t a r g e t = 0.041 and other parameter values can be found in Fig. 4 and Table I.

Close modal

However, a nuanced observation is crucial: once the proportion of human activity regions reaches a certain threshold, the incremental impact on the relative error becomes marginal. This observation aligns with the conclusions drawn from Fig. 4: beyond a certain point, further increasing the proportion of human activity regions fails to induce a noticeable effect on pattern transformation. Additionally, when the percentage of human activity areas remains constant, we observe an influence of their number on the relative error: an escalation in the number of human activity areas correlates with a reduction in the relative error. Similarly, beyond a certain quantity, further increases in the number of human activity areas result in minimal alterations to the relative error.

In summary, the interplay between the number and proportion of human activity areas significantly influences the relative error. When the number of human activity areas is limited, a substantial proportion of human activity regions is imperative to ensure a small relative error, indicative of a successful conversion. Conversely, with a larger number of human activity areas, a minor proportion suffices to maintain a low relative error. This intricate relationship underscores the necessity for a nuanced understanding of the interdependence between the number and proportion of human activity areas in ensuring effective pattern conversion.

This section delves into the influence of regional boundaries on pattern shifts. While previous research has predominantly focused on regular regions, a question arises: what happens when the shape of the region changes, as in the case of a circular area? To explore this, we adjusted the shape of the region boundaries, employing numerical simulations with a circular area as an illustrative example. Here, we set the radius of the circular area to 50⁠, ensuring its diameter corresponds to the side length of the regular region. The results indicate that altering the boundary shape does not impact the pattern transformation (see Figs. 8 and 9). Additionally, mirroring our observations in regular regions, we examined the impact of both the proportion and the number of human activity areas on pattern structure transformations within the circular region. First, concerning the proportion of human activity areas, maintaining a fixed quantity of areas ( 9 )⁠, a gradual increase in proportion visibly enhances the transformation of vegetation patterns. Second, focusing on the number of human activity areas, with a fixed proportion of 36%, incrementally increasing the quantity ( 1 , 4 , 9 , 16 ) similarly improves the transformation of vegetation patterns. These findings align with conclusions drawn from regular regions (see Figs. 4 and 5), emphasizing the methodological robustness and independence of our results from changes in regional boundaries.

FIG. 8.
Transformation of vegetation patterns for different proportions of human activity areas when the number of human activity areas is constant ( 9 ). The percentage at the top represents the corresponding proportion of human activity areas, increasing gradually from left to right. (a1)–(a4) The spatial distribution of vegetation. (b1)–(b4) The spatial distribution of human activity intensity. Here, δ t a r g e t = 0.041 and other parameter values can be found in Fig. 4 and Table I.
View largeDownload slide

Transformation of vegetation patterns for different proportions of human activity areas when the number of human activity areas is constant ( 9 )⁠. The percentage at the top represents the corresponding proportion of human activity areas, increasing gradually from left to right. (a1)–(a4) The spatial distribution of vegetation. (b1)–(b4) The spatial distribution of human activity intensity. Here, δ t a r g e t = 0.041 and other parameter values can be found in Fig. 4 and Table I.

FIG. 8.
Transformation of vegetation patterns for different proportions of human activity areas when the number of human activity areas is constant ( 9 ). The percentage at the top represents the corresponding proportion of human activity areas, increasing gradually from left to right. (a1)–(a4) The spatial distribution of vegetation. (b1)–(b4) The spatial distribution of human activity intensity. Here, δ t a r g e t = 0.041 and other parameter values can be found in Fig. 4 and Table I.
View largeDownload slide

Transformation of vegetation patterns for different proportions of human activity areas when the number of human activity areas is constant ( 9 )⁠. The percentage at the top represents the corresponding proportion of human activity areas, increasing gradually from left to right. (a1)–(a4) The spatial distribution of vegetation. (b1)–(b4) The spatial distribution of human activity intensity. Here, δ t a r g e t = 0.041 and other parameter values can be found in Fig. 4 and Table I.

Close modal
FIG. 9.
When the proportion of human activity areas remains unchanged ( 36 %), the vegetation patterns of different numbers of human activity areas change. (a1)–(a4) The spatial distribution of vegetation. (b1)–(b4) The spatial distribution of human activity intensity. Each small circle represents a human activity area, with the number of human activity areas from left to right being 1, 4, 9, and 16, respectively. Here, δ t a r g e t = 0.06 and other parameter values can be found in Fig. 4 and Table I.
View largeDownload slide

When the proportion of human activity areas remains unchanged (⁠ 36 %⁠), the vegetation patterns of different numbers of human activity areas change. (a1)–(a4) The spatial distribution of vegetation. (b1)–(b4) The spatial distribution of human activity intensity. Each small circle represents a human activity area, with the number of human activity areas from left to right being 1, 4, 9, and 16, respectively. Here, δ t a r g e t = 0.06 and other parameter values can be found in Fig. 4 and Table I.

FIG. 9.
When the proportion of human activity areas remains unchanged ( 36 %), the vegetation patterns of different numbers of human activity areas change. (a1)–(a4) The spatial distribution of vegetation. (b1)–(b4) The spatial distribution of human activity intensity. Each small circle represents a human activity area, with the number of human activity areas from left to right being 1, 4, 9, and 16, respectively. Here, δ t a r g e t = 0.06 and other parameter values can be found in Fig. 4 and Table I.
View largeDownload slide

When the proportion of human activity areas remains unchanged (⁠ 36 %⁠), the vegetation patterns of different numbers of human activity areas change. (a1)–(a4) The spatial distribution of vegetation. (b1)–(b4) The spatial distribution of human activity intensity. Each small circle represents a human activity area, with the number of human activity areas from left to right being 1, 4, 9, and 16, respectively. Here, δ t a r g e t = 0.06 and other parameter values can be found in Fig. 4 and Table I.

Close modal

In this study, we incorporated optimal control theory into a dryland vegetation water model with cross-diffusion and proposed a new method for studying dryland vegetation restoration. Our investigation unfolds on two fronts: theoretical studies and numerical simulations. First, we have tackled the problem comprehensively by considering local distribution control and boundary control, constructing the corresponding optimal control problem, and conducting a systematic theoretical analysis. This analysis has encompassed outlining the first-order necessity optimality conditions as well as the second-order optimality conditions. To complement our theoretical analyses, we have performed numerical simulations that shed light on the influence of the proportion and number of human activity areas on the transformation of different vegetation patterns.

Our primary focus has been on developing strategies for vegetation restoration in drylands, with the overarching idea of optimizing human activities to shift the vegetation distribution from a structure with high desertification potential to one with low potential, thereby promoting restoration. To this end, we introduced the Regional Vegetation Desertification Level Index ( R V D L I ) and calculated it for various vegetation pattern structures. The results showed that the green-eye pattern had the highest R V D L I⁠, which is consistent with prior research that spot patterns serve as desertification warning signals.8 Conversely, the R V D L I of stripe and white-eye patterns exhibited a decreasing trend, making them suitable restoration targets.

Subsequently, we use the transition from the green-eye pattern to the white-eye pattern as an illustrative case to examine the influence of the proportion of human activity areas on this transition (Fig. 4). Our research results revealed that as the proportion of human activity areas increased, the efficacy of converting between distinct vegetation pattern structures exhibited significant improvement, reaching a notable threshold. Beyond this threshold, further increments in the proportion of human activity areas ceased to yield substantial impacts and, instead, incurred heightened costs. Furthermore, while maintaining a consistent proportion of human activity areas, we scrutinized the repercussions of altering the number of such areas on the transition from a green-eye pattern to a stripe pattern. The results underscored that a greater number of human activity areas, i.e., the more dispersed the human activity areas, notably enhanced the transformation between vegetation pattern structures. This observation was corroborated by the R V D L I data, which consistently revealed diminished R V D L I values for the corresponding vegetation patterns associated with more dispersed human activity areas. Moreover, Fig. 6 reveals that irrespective of the dispersal of human activity areas, the R V D L I of the converted vegetation patterns is substantially reduced compared to their pre-conversion states, underscoring the feasibility of optimizing human activity as a strategy for dryland vegetation restoration.

Due to the challenging climatic conditions and vulnerable ecosystems prevalent in arid regions, numerous areas have witnessed significant vegetation degradation. Consequently, there has been a heightened focus on initiatives aimed at vegetation restoration. Mau et al.,42 for instance, employed periodic modulation of the vegetation landscape through spatial resonance to effectively counteract desertification. However, their findings primarily pertain to spatially forced systems approaching finite-wave-number instabilities under noisy conditions. In contrast, the optimal control methodology employed in this paper does not suffer from the same limitations and finds wide-ranging applications in diverse fields, including population dynamics and epidemiology. For instance, Chang et al.37 utilized sparse control techniques to investigate the control of infectious disease outbreaks. In this paper, our research primarily centers on the principles of local control and boundary control. Local control, in particular, plays a pivotal role in achieving cost efficiency, allowing us to confine the regulation of human activity to specific localized areas. While sparse control can also deliver cost benefits, it is constrained by the stochastic nature of the sparsely controlled region, rendering it incapable of precisely defining the control region. In contrast, local control grants us the ability to pinpoint and delineate the areas subject to regulation.

In conclusion, our research enhances our understanding of how human activities influence changes in vegetation patterns, providing insights for future ecological conservation planning. Desertification management is a long-term endeavor; in addition to the optimal control problem considered in this paper, we can also consider constructing multi-stage optimal control problems, gradually reversing desertification by optimizing human activities in stages. Furthermore, integrating remote sensing and geographic information systems (GIS) technology can strengthen the monitoring and evaluation of vegetation restoration efforts. Remote sensing technologies, such as satellite imagery and drones, can offer detailed and large-scale observations of vegetation changes, while GIS tools can analyze spatial patterns and trends. This integration will provide real-time data for adaptive management strategies, ensuring more accurate and timely interventions. Combined, these technologies will enable more precise tracking of restoration progress and identify areas requiring immediate attention. This idea will be a central focus of our future work, contributing to developing more effective desertification mitigation strategies.

This work was supported by the National Natural Science Foundation of China under Grant Nos. 42275034 and 42075029, the Scientific Activities of Selected Returned Overseas Professionals in Shanxi Province (No. 20230013), and the Graduate Research Innovation Project of Shanxi Province (No. 2024KY017).

The authors have no conflicts to disclose.

Li-Feng Hou: Conceptualization (equal); Data curation (equal); Writing – original draft (equal). Shu-Peng Gao: Conceptualization (equal); Data curation (equal); Writing – original draft (equal). Li-Li Chang: Conceptualization (equal); Data curation (equal). Yong-Ping Wu: Conceptualization (equal); Writing – review & editing (equal). Guo-Lin Feng: Conceptualization (equal); Writing – review & editing (equal). Zhen Wang: Writing – original draft (equal); Writing – review & editing (equal). Gui-Quan Sun: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Writing – original draft (equal); Writing – review & editing (equal).

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

1. The first-order necessity condition

Theorem A.1
Suppose that ( u ∗ , b n ∗ , b w ∗ ) is a local solution of the optimal control problem (2)–(3), ( n ∗ , w ∗ ) is the solution of the corresponding state system (3), and let ( p , q ) be the solution of the following adjoint system:
(A1)
Then, the variational inequality holds:
(A2)
Proof.
The mapping ( u , b n , b w ) ↦ ( n , w ) from K f e to L 2 ( 0 , T ; H 1 ( Ω ) 2 ) is defined as a solution of the state system (3), which has G a ^ t e a u x derivatives in every direction u ^⁠, b n ^⁠, and b w ^ in the tangential cone Tan K fe ( u )⁠, Tan K fe ( b n )⁠, and Tan K fe ( b w ) of K f e at u⁠, b n⁠, and b w⁠.32,43 Therefore, we can get that the directional derivative of the system (3) at ( n ∗ , w ∗ , u ∗ , b n ∗ , b w ∗ ) along u ^⁠, b n ^⁠, and b w ^⁠, respectively, as
(A3)
(A4)
and
(A5)
where ( n 1 , w 1 ) = ( d n d u , d w d u ) ⋅ u ^⁠, ( n 2 , w 2 ) = ( d n d b n , d w d b n ) ⋅ b n ^⁠, and ( n 3 , w 3 ) = ( d n d b w , d w d b w ) ⋅ b w ^⁠.
The first equation of the above three systems (A3)–(A5) is, respectively, an inner product with p⁠, and integration over Ω T yields
(A6)
(A7)
(A8)
Similarly, adding q as an inner product to the second equation of the system (A3)–(A5) and integrating it over Ω T yields
(A9)
(A10)
(A11)
By adding formula (A6) to formula (A9), formula (A7) to formula (A10), and formula (A8) to formula (A11), it can be obtained that
∫ Ω [ p ( x , T ) n 1 ( x , T ) + q ( x , T ) w 1 ( x , T ) ] d x = ∫ 0 T ∫ Ω χ ω u ^ d x d t , ∫ Ω [ p ( x , T ) n 2 ( x , T ) + q ( x , T ) w 2 ( x , T ) ] d x = ∫ 0 T ∫ ∂ Ω b n ^ p d x d t , ∫ Ω [ p ( x , T ) n 3 ( x , T ) + q ( x , T ) w 3 ( x , T ) ] d x = ∫ 0 T ∫ Ω b w ^ q d x d t .
Next, we compute the G a ^ t e a u x derivative of O [ n ∗ , w ∗ , u ∗ , b n ∗ , b w ∗ ] in the direction u ^⁠, b n ^⁠, and b w ^⁠,
0 ≤ O u [ n ∗ , w ∗ , u ∗ , b n ∗ , b w ∗ ] ⋅ u ^ = ∫ Ω κ 1 [ n ∗ ( x , T ) − n t a r ( x ) ] n 1 d x + ∫ Ω κ 2 [ w ∗ ( x , T ) − w t a r ( x ) ] w 1 d x + ∫ 0 T ∫ ω κ 3 u ∗ u ^ d x d t = ∫ 0 T ∫ Ω χ ω ( κ 3 u ∗ + p ) u ^ d x d t ,
0 ≤ O b n [ n ∗ , w ∗ , u ∗ , b n ∗ , b w ∗ ] ⋅ b n ^ = ∫ Ω κ 1 [ n ∗ ( x , T ) − n t a r ( x ) ] n 1 d x + ∫ Ω κ 2 [ w ∗ ( x , T ) − w t a r ( x ) ] w 1 d x + ∫ 0 T ∫ ∂ Ω κ 4 b n ∗ b n ^ d s d t = ∫ 0 T ∫ ∂ Ω ( κ 4 b n ∗ + p ) b n ^ d s d t ,
0 ≤ O b w [ n ∗ , w ∗ , u ∗ , b n ∗ , b w ∗ ] ⋅ b w ^ = ∫ Ω κ 1 [ n ∗ ( x , T ) − n t a r ( x ) ] n 1 d x + ∫ Ω κ 2 [ w ∗ ( x , T ) − w t a r ( x ) ] w 1 d x + ∫ 0 T ∫ ∂ Ω κ 5 b w ∗ b w ^ d s d t = ∫ 0 T ∫ ∂ Ω ( κ 5 b w ∗ + q ) b w ^ d s d t .

We can obtain our variational inequality (A2) by replacing u ^⁠, b n ^⁠, and b w ^ in the above equation with u − u ∗⁠, b n − b n ∗⁠, and b w − b w ∗⁠, respectively.

2. Second-order optimality conditions

In this subsection, we formulate second-order optimality conditions for our control problem (2)–(3) in terms of second-order F r e ´ c h e t derivatives of the associated Lagrange function. The Lagrange function L for the problem is defined as follows:
L [ n , w , u , b n , b w , p , q ] = J [ n , w , u , b n , b w ] + ∫ 0 T ∫ Ω [ − ∂ n ∂ t + Δ n + g 1 ( n , w ) + χ ω u ] p d x d t + ∫ 0 T ∫ Ω [ − ∂ w ∂ t + μ Δ w − μ δ Δ n + g 2 ( n , w ) ] q d x d t + ∫ 0 T ∫ ∂ Ω ( − ∂ n ∂ ν − α 1 n + b n ) p d s d t ∫ 0 T ∫ ∂ Ω ( − μ ∂ w ∂ ν + μ δ ∂ n ∂ ν − α 2 w + b w ) q d s d t ,
where J [ n , w , u , b n , b w ] is the objective functional of the control problem, and p and q are the adjoint variables corresponding to the state equations for n and w⁠, respectively.
The first-order necessity condition (A2) yields the representation
( u ∗ , b n ∗ , b w ∗ ) = { ( u 1 , b n 1 , b w 1 ) if χ ω ( κ 3 u ∗ + p ) > 0 , κ 4 b n ∗ + p > 0 , κ 5 b w ∗ + q > 0 , ( u 2 , b n 2 , b w 2 ) if χ ω ( κ 3 u ∗ + p ) < 0 , κ 4 b n ∗ + p < 0 , κ 5 b w ∗ + q < 0.
Thus, ( u ∗ , b n ∗ , b w ∗ ) is determined by the set of all ( x , t ) such that
{ | χ ω ( κ 3 u ∗ ( x , t ) + p ( x , t ) ) | > 0 , | κ 4 b n ∗ ( x , t ) + p ( x , t ) | > 0 , | κ 5 b w ∗ ( x , t ) + q ( x , t ) | > 0 ,
and the higher-order conditions are only of interest for its complement. We now introduce the sets of strongly active control constraints for ( u ∗ , b n ∗ , b w ∗ ) by
A τ Ω T ( u ∗ ) = { ( x , t ) ∈ Ω T : | χ ω ( κ 3 u ∗ ( x , t ) + p ( x , t ) ) | > τ , f o r ∀ τ ≥ 0 } , A τ Σ T ( b n ∗ , b w ∗ ) = { ( x , t ) ∈ Σ T : | κ 4 b n ∗ + p | > τ a n d | κ 5 b w ∗ + q | > τ f o r ∀ τ ≥ 0 } .
Based on this, we define the τ-critical cone C τ ( u ∗ , b n ∗ , b w ∗ )⁠, which is the set of all ( u , b n , b w ) ∈ L ∞ ( Ω T ) × L ∞ ( Σ T ) × L ∞ ( Σ T ) satisfying
(A12)
Theorem A.2
Suppose that ( u ∗ , b n ∗ , b w ∗ ) is a local solution of the optimal control problem (2)–(3) and ( p , q ) is the corresponding adjoint variable, then
(A13)
where
L ′ ′ ( n ∗ , w ∗ , u ∗ , b n ∗ , b w ∗ , p , q ) ( n , w , u , b n , b w ) 2 = ∫ Ω [ κ 1 n ( x , T ) 2 + κ 2 w ( x , T ) 2 ] d x + ∫ 0 T ∫ Ω [ ( 2 w ∗ n 2 + 4 n ∗ n w ) ( p − q ) + κ 3 χ ω u 2 ] d x d t + ∫ 0 T ∫ ∂ Ω ( κ 4 b n 2 + κ 5 b w 2 ) d s d t .
Proof.
Let ( u , b n , b w ) ∈ C τ ( u ∗ , b n ∗ , b w ∗ ) be arbitrary. Note that we cannot guarantee that ( u ∗ + t u ^ , b n ∗ + t b n ^ , b w ∗ + t b w ^ ) ∈ K f e even for sufficiently small t > 0⁠. We, hence, introduce the set
I m Ω T = { ( x , t ) ∈ Ω T : u 1 ( x , t ) + 1 m ≤ u ∗ ≤ u 2 ( x , t ) − 1 m } , I m Σ T = { ( x , t ) ∈ Σ T : b n 1 ( x , t ) + 1 m ≤ b n ∗ ≤ b n 2 ( x , t ) − 1 m , b w 1 ( x , t ) + 1 m ≤ b w ∗ ≤ b w 2 ( x , t ) − 1 m } ,
where m ∈ N and consider the functions ( u m , b n m , b w m ) := χ m ( u , b n , b w )⁠, where
χ m = { 1 if ( x , t ) ∈ I m Ω T ∪ I m Σ T or { u ∗ ∈ [ u 1 , u 2 ] and u 2 − u 1 ≥ 1 / m ; b n ∗ ∈ [ b n 1 , b n 2 ] and b n 2 − b n 1 ≥ 1 / m ; b w ∗ ∈ [ b w 1 , b w 2 ] and b w 2 − b w 1 ≥ 1 / m } , 0 if u ∗ ∈ ( u 1 , u 1 + 1 / m ) ∪ ( u 2 − 1 / m , u 2 ) ; b n ∗ ∈ ( b n 1 , b n 1 + 1 / m ) ∪ ( b n 2 − 1 / m , b n 2 ) ; b w ∗ ∈ ( b w 1 , b w 1 + 1 / m ) ∪ ( b w 2 − 1 / m , b w 2 ) .
Hence, we have ( u ^ , b n ^ , b w ^ ) = ( u ∗ + t u m , b n ∗ + t b n m , b w ∗ + t b w m ) ∈ K f e for adequately small t > 0⁠. Consequently,
0 ≤ L ( n , w , u ^ , b n ^ , b w ^ ) − L ( n ∗ , w ∗ , u ∗ , b n ∗ , b w ∗ ) = L ′ ( n ∗ , w ∗ , u ∗ , b n ∗ , b w ∗ ) ( n , w , u m , b n m , b w m ) T t + 1 2 ! L ′ ′ ( n ∗ , w ∗ , u ∗ , b n ∗ , b w ∗ ) ( n , w , u m , b n m , b w m ) 2 t 2 + r 2 L ,
where the specific expressions for L ′ ( n ∗ , w ∗ , u ∗ , b n ∗ , b w ∗ ) and L ′ ′ ( n ∗ , w ∗ , u ∗ , b n ∗ , b w ∗ ) can be found in  Appendix A, ( n − n ∗ , w − w ∗ ) = ( t n , t w )⁠, and r 2 L represents the second-order remainder of L⁠. Combining the variational inequality (A2) and the fact that ( u , b n , b w ) ∈ C τ ( u ∗ , b n ∗ , b w ∗ )⁠, it follows that ( u , b n , b w ) and ( u m , b n m , b w m ) vanish at almost all points where
χ ω ( κ 3 u ∗ + p ) ≠ 0 , κ 4 b n ∗ + p ≠ 0 , κ 5 b w ∗ + q ≠ 0.
Therefore,
L ′ ( n ∗ , w ∗ , u ∗ , b n ∗ , b w ∗ ) ( n , w , u m , b n m , b w m ) T = ( χ ω ( κ 3 u ∗ + p ) , u m ) L 2 ( Ω T ) + ( κ 4 b n ∗ + p , b n m ) L 2 ( Σ T ) + ( κ 5 b w ∗ + q , b w m ) L 2 ( Σ T ) = 0
and
0 ≤ L ′ ′ ( n ∗ , w ∗ , u ∗ , b n ∗ , b w ∗ ) ( n , w , u m , b n m , b w m ) 2 + t − 2 r 2 L .
The limit at t ↓ 0 shows that
0 ≤ L ( n , w , u , b n , b w ) ′ ′ ( n , w , u m , b n m , b w m ) 2 .
Notice that ( u m , b n m , b w m ) → ( u , b n , b w ) is almost everywhere as m → ∞ and that ( u m , b n m , b w m ) ≤ χ m ( u , b n , b w ) is everywhere pointwise for all m ∈ N⁠. Thus, we can use Lebègue’s dominated convergence theorem to derive ( u m , b n m , b w m ) → ( u , b n , b w ) in L ∞ ( Ω T ) × L ∞ ( Σ T ) × L ∞ ( Σ T )⁠. Thus, we can infer that
L ′ ′ ( n ∗ , w ∗ , u ∗ , b n ∗ , b w ∗ , p , q ) ( n , w , u , b n , b w ) 2 ≥ 0 f o r ∀ ( u , b n , b w ) ∈ C τ ( u ∗ , b n ∗ , b w ∗ )
as m → ∞⁠.
Theorem A.3
Let ( u ∗ , b n ∗ , b w ∗ ) satisfy the first-order necessary optimal condition of Theorem 3.3 and all the constraints of problem (2) and (3), and there exists a constant λ , τ > 0⁠, such that condition
(A14)
holds. Then, there is a constant ϵ , σ > 0 such that we have the quadratic growth condition
(A15)
for ∀ ( u , b n , b w ) ∈ L ∞ ( Ω T ) × L ∞ ( Σ T ) × L ∞ ( Σ T ) with ‖ u − u ∗ ‖ L ∞ ( Ω T ) 2 + ‖ b n − b n ∗ ‖ L ∞ ( Σ T ) 2 + ‖ b w − b w ∗ ‖ L ∞ ( Σ T ) 2 ≤ ϵ⁠. Specifically, ( u ∗ , b n ∗ , b w ∗ ) is locally optimal.
Proof.
Let ( u ∗ , b n ∗ , b w ∗ ) ∈ L ∞ ( Ω T ) × L ∞ ( Σ T ) × L ∞ ( Σ T ) with
‖ u − u ∗ ‖ L ∞ ( Ω T ) 2 + ‖ b n − b n ∗ ‖ L ∞ ( Σ T ) 2 + ‖ b w − b w ∗ ‖ L ∞ ( Σ T ) 2 ≤ ϵ .
We have
(A16)
where h = ( n − n ∗ , w − w ∗ , u − u ∗ , b n − b n ∗ , b w − b w ∗ )⁠. To use conditions (A14), we perform the decomposition h := h 0 + h 00 with h 0 = ( n 0 − n ∗ , w 0 − w ∗ , u 0 − u ∗ , b n 0 − b n ∗ , b w 0 − b w ∗ ) and h 00 = ( n 00 − n ∗ , w 00 − w ∗ , u 00 − u ∗ , b n 00 − b n ∗ , b w 00 − b w ∗ )⁠, together with
( u 0 , b n 0 , b w 0 ) := { ( u , b n , b w ) if ( x , t ) ∉ A τ Ω T ( u ∗ ) ∪ A τ Σ T ( b n ∗ , b w ∗ ) , 0 if ( x , t ) ∈ A τ Ω T ( u ∗ ) ∪ A τ Σ T ( b n ∗ , b w ∗ ) .
This yields ( u 00 , b n 00 , b w 00 ) = 0 on ( Ω T ∖ A τ Ω T ) ∪ ( Σ T ∖ A τ Σ T )⁠. Obviously, ( u 0 , b n 0 , b w 0 ) satisfies the sign condition of the critical cone; therefore, ( u 0 , b n 0 , b w 0 ) ∈ C τ ( u ∗ , b n ∗ , b w ∗ )⁠. Next, we estimate L ′ ′ ( n ∗ , w ∗ , u ∗ , b n ∗ , b w ∗ ) h 2⁠. Since
(A17)
we estimate the three terms separately.
For L ′ ′ ( n ∗ , w ∗ , u ∗ , b n ∗ , b w ∗ ) h 0 2⁠, invoking (A14) and ( u 0 , b n 0 , b w 0 ) ∈ C τ ( u ∗ , b n ∗ , b w ∗ )⁠, we get
(A18)
For L ′ ′ ( n ∗ , w ∗ , u ∗ , b n ∗ , b w ∗ ) [ h 0 , h 00 ]⁠, we have
(A19)
with c , c 1 > 0⁠. Similarly,
(A20)
To sum up, combining the inequalities (A18)–(A20) yields
(A21)
Now, we choose sufficiently small ϵ > 0 so that ϵ ( c 1 + c 2 ) < τ / 2⁠. Since ( u 00 , b n 00 , b w 00 ) = 0 on ( Ω T ∖ A τ Ω T ) ∪ ( Σ T ∖ A τ Σ T )⁠, we can deduce that
(A22)
Substituting estimates (A21) and (A22) into (A16) then yields
O ( n , w , u , b n , b w ) − O ( n ∗ , w ∗ , u ∗ , b n ∗ , b w ∗ ) ≥ τ ( ∫ A τ Ω T | u − u ∗ | + ∫ A τ Σ T | b n − b n ∗ | + ∫ A τ Σ T | b w − b w ∗ | ) − τ 2 ( ∫ A τ Ω T | u − u ∗ | + ∫ A τ Σ T | b n − b n ∗ | + ∫ A τ Σ T | b w − b w ∗ | ) + λ 4 ( ‖ u 0 − u ∗ ‖ L 2 ( Ω T ) 2 + ‖ b n 0 − b n ∗ ‖ L 2 ( Σ T ) 2 + ‖ b w 0 − b w ∗ ‖ L 2 ( Σ T ) 2 ) + r 2 L ≥ τ 2 ( ∫ A τ Ω T | u − u ∗ | + ∫ A τ Σ T | b n − b n ∗ | + ∫ A τ Σ T | b w − b w ∗ | ) + λ 4 ( ‖ u 0 − u ∗ ‖ L 2 ( Ω T ) 2 + ‖ b n 0 − b n ∗ ‖ L 2 ( Σ T ) 2 + ‖ b w 0 − b w ∗ ‖ L 2 ( Σ T ) 2 ) + r 2 L .
Without loss of generality, let us now set ϵ ≤ 1⁠. Then, by the definition of | u − u ∗ |⁠, | b n − b n ∗ |⁠, and | b w − b w ∗ |⁠, it follows that | u − u ∗ | ≥ ( u − u ∗ ) 2⁠, | b n − b n ∗ | ≥ ( b n − b n ∗ ) 2⁠, and | b w − b w ∗ | ≥ ( b w − b w ∗ ) 2⁠. Since
‖ u 0 − u ∗ ‖ L 2 ( Ω T ) 2 + ‖ b n 0 − b n ∗ ‖ L 2 ( Σ T ) 2 + ‖ b w 0 − b w ∗ ‖ L 2 ( Σ T ) 2 = ∫ Ω T ∖ A τ Ω T ( u − u ∗ ) L 2 ( Ω T ) 2 + ∫ Σ T ∖ A τ Σ T ( b n − b n ∗ ) L 2 ( Σ T ) 2 + ∫ Σ T ∖ A τ Σ T ( b w − b w ∗ ) L 2 ( Σ T ) 2 .
We can further obtain
O ( n , w , u , b n , b w ) − O ( n ∗ , w ∗ , u ∗ , b n ∗ , b w ∗ ) ≥ τ 2 [ ∫ A τ Ω T ( u − u ∗ ) 2 + ∫ A τ Σ T ( b n − b n ∗ ) 2 + ∫ A τ Σ T ( b w − b w ∗ ) 2 ] + λ 4 ( ‖ u − u ∗ ‖ L 2 ( Ω T ∖ A τ Ω T ) 2 + ‖ b n − b n ∗ ‖ L 2 ( S i g m a T ∖ A τ Σ T ) 2 + ‖ b w − b w ∗ ‖ L 2 ( Σ T ∖ A τ Σ T ) 2 ) + r 2 L ≥ min { τ 2 , λ 4 } ( ‖ u − u ∗ ‖ L 2 ( Ω T ) 2 + ‖ b n − b n ∗ ‖ L 2 ( Σ T ) 2 + ‖ b w − b w ∗ ‖ L 2 ( Σ T ) 2 ) + r 2 L .
Also, the remainder r 2 L satisfies
r 2 L ( n ∗ , w ∗ , u ∗ , b n ∗ , b w ∗ , h ) ‖ u − u ∗ ‖ L 2 ( Ω T ) 2 + ‖ b n − b n ∗ ‖ L 2 ( Σ T ) 2 + ‖ b w − b w ∗ ‖ L 2 ( Σ T ) 2 → 0
as
‖ u − u ∗ ‖ L ∞ ( Ω T ) 2 + ‖ b n − b n ∗ ‖ L ∞ ( Σ T ) 2 + ‖ b w − b w ∗ ‖ L ∞ ( Σ T ) 2 → 0.
Thus, for sufficiently small ϵ > 0⁠, we can obtain
O ( n , w , u , b n , b w ) − O ( n ∗ , w ∗ , u ∗ , b n ∗ , b w ∗ ) ≥ min { τ 2 , λ 4 } ( ‖ u − u ∗ ‖ L 2 ( Ω T ) 2 + ‖ b n − b n ∗ ‖ L 2 ( Σ T ) 2 + ‖ b w − b w ∗ ‖ L 2 ( Σ T ) 2 ) = σ ( ‖ u − u ∗ ‖ L 2 ( Ω T ) 2 + ‖ b n − b n ∗ ‖ L 2 ( Σ T ) 2 + ‖ b w − b w ∗ ‖ L 2 ( Σ T ) 2 ) .

In the computation of the optimal control problem, the projected subgradient algorithm with a variable step is mainly used; see the Algorithm for the specific procedure.

Algorithm:

Variable step projection gradient method.

1:  Initialization 
2:  Select initial guess for ( u 0 , b n 0 , b w 0 )⁠, 
  initial step size θ = 1, 
  r e l E r r o r T o l = 1.0 × 10 − 3⁠, h: = 0; 
  Solve the state variables (n0, w0); 
  Solve the adjoint variables (p0, q0); 
  Solve the objective functional O0; 
  Solve the gradient ∇ u O 0⁠, ∇ b n O 0⁠, and ∇ b w O 0⁠; 
  Calculate the relative error n0 and ntar: 
  r e l E r r o r 0 = ‖ n 0 − n t a r ‖ L 2 ( Ω ) ‖ n t a r ‖ L 2 ( Ω )⁠. 
3:  Set m a x c y c l e = 1000⁠. 
4:  end initialization 
5:  main cycle 
6:   while r e l E r r o r h > r e l E r r o r T o l & & h < m a x c y c l e do 
7:             oldu ← uh, o l d b n h ← b n h⁠, o l d b w h ← b w h⁠; 
              o l d ∇ u O ← ∇ u O h⁠, 
              o l d ∇ b n O ← ∇ b n O h⁠, 
              o l d ∇ b w O ← ∇ b w O h⁠; 
8:             h ← h + 1; 
9:  update control: 
           u h = P [ u 1 , u 2 ] ( u h − 1 − θ ∇ u O h − 1 )⁠, 
           b n h = P [ b n 1 , b n 2 ] ( b n h − 1 − θ ∇ b n O h − 1 )⁠, 
           b w h = P [ b w 1 , b w 2 ] ( b w h − 1 − θ ∇ b w O h − 1 )⁠; 
10:  update state variables (nh, wh); 
11:  update adjoint variables (ph, qh); 
12:  update objective functional Oh; 
13:  update gradient ∇ u O h⁠, ∇ b n h⁠, ∇ b w h⁠; 
14:  update r e l E r r o r h = ‖ n h − n t a r ‖ L 2 ( Ω ) ‖ n t a r ‖ L 2 ( Ω ) ; 
15:  if |Oh − Oh−1| < 10−10  or 
              | r e l E r r o r h − r e l E r r o r h − 1 | < 10 − 10  or 
             θ < 10−5 do 
16:    break; 
17:  end if 
18:  if Oh > Oh−1 do 
19:    h ← h − 1; 
20:    θ ← θ/2; 
21:    uh ← oldu, 
    b n h ← o l d b n⁠, 
    b w h ← o l d b w⁠; 
22:    ∇ u O h ← o l d ∇ u O⁠, 
    ∇ b n O h ← o l d ∇ b n h O⁠, 
    ∇ b w O h ← o l d ∇ b w h O⁠; 
23:   end if 
24:     end while 
25:  end main cycle 
1:  Initialization 
2:  Select initial guess for ( u 0 , b n 0 , b w 0 )⁠, 
  initial step size θ = 1, 
  r e l E r r o r T o l = 1.0 × 10 − 3⁠, h: = 0; 
  Solve the state variables (n0, w0); 
  Solve the adjoint variables (p0, q0); 
  Solve the objective functional O0; 
  Solve the gradient ∇ u O 0⁠, ∇ b n O 0⁠, and ∇ b w O 0⁠; 
  Calculate the relative error n0 and ntar: 
  r e l E r r o r 0 = ‖ n 0 − n t a r ‖ L 2 ( Ω ) ‖ n t a r ‖ L 2 ( Ω )⁠. 
3:  Set m a x c y c l e = 1000⁠. 
4:  end initialization 
5:  main cycle 
6:   while r e l E r r o r h > r e l E r r o r T o l & & h < m a x c y c l e do 
7:             oldu ← uh, o l d b n h ← b n h⁠, o l d b w h ← b w h⁠; 
              o l d ∇ u O ← ∇ u O h⁠, 
              o l d ∇ b n O ← ∇ b n O h⁠, 
              o l d ∇ b w O ← ∇ b w O h⁠; 
8:             h ← h + 1; 
9:  update control: 
           u h = P [ u 1 , u 2 ] ( u h − 1 − θ ∇ u O h − 1 )⁠, 
           b n h = P [ b n 1 , b n 2 ] ( b n h − 1 − θ ∇ b n O h − 1 )⁠, 
           b w h = P [ b w 1 , b w 2 ] ( b w h − 1 − θ ∇ b w O h − 1 )⁠; 
10:  update state variables (nh, wh); 
11:  update adjoint variables (ph, qh); 
12:  update objective functional Oh; 
13:  update gradient ∇ u O h⁠, ∇ b n h⁠, ∇ b w h⁠; 
14:  update r e l E r r o r h = ‖ n h − n t a r ‖ L 2 ( Ω ) ‖ n t a r ‖ L 2 ( Ω ) ; 
15:  if |Oh − Oh−1| < 10−10  or 
              | r e l E r r o r h − r e l E r r o r h − 1 | < 10 − 10  or 
             θ < 10−5 do 
16:    break; 
17:  end if 
18:  if Oh > Oh−1 do 
19:    h ← h − 1; 
20:    θ ← θ/2; 
21:    uh ← oldu, 
    b n h ← o l d b n⁠, 
    b w h ← o l d b w⁠; 
22:    ∇ u O h ← o l d ∇ u O⁠, 
    ∇ b n O h ← o l d ∇ b n h O⁠, 
    ∇ b w O h ← o l d ∇ b w h O⁠; 
23:   end if 
24:     end while 
25:  end main cycle 
View Large
1.
H. N.
Le Houérou
,
J. Arid Environ.
34
,
133
(
1996
).
https://doi.org/10.1006/jare.1996.0099
Google Scholar
Crossref
Search ADS
 
2.
J. F.
Reynolds
,
D. M. S.
Smith
,
E. F.
Lambin
,
B.
Turner
,
M.
Mortimore
,
S. P.
Batterbury
,
T. E.
Downing
,
H.
Dowlatabadi
,
R. J.
Fernández
,
J. E.
Herrick
 et al.,
Science
316
,
847
(
2007
).
https://doi.org/10.1126/science.1131634
Google Scholar
Crossref
Search ADS
PubMed
 
4.
J.
Thiery
,
J.-M.
d’Herbès
, and
C.
Valentin
,
J. Ecol.
83
,
497
(
1995
).
https://doi.org/10.2307/2261602
Google Scholar
Crossref
Search ADS
 
5.
E.
Meron
,
Nonlinear Physics of Ecosystems
(
CRC Press
,
2015
).
Google Scholar
Crossref
Search ADS
 
6.
S.
Getzin
,
H.
Yizhaq
,
B.
Bell
,
T. E.
Erickson
,
A. C.
Postle
,
I.
Katra
,
O.
Tzuk
,
Y. R.
Zelnik
,
K.
Wiegand
,
T.
Wiegand
 et al.,
Proc. Natl. Acad. Sci. U.S.A.
113
,
3551
(
2016
).
https://doi.org/10.1073/pnas.1522130113
Google Scholar
Crossref
Search ADS
PubMed
 
8.
J.
von Hardenberg
,
E.
Meron
,
M.
Shachak
, and
Y.
Zarmi
,
Phys. Rev. Lett.
87
,
198101
(
2001
).
https://doi.org/10.1103/PhysRevLett.87.198101
Google Scholar
Crossref
Search ADS
PubMed
 
9.
M.
Rietkerk
,
M. C.
Boerlijst
,
F.
van Langevelde
,
R.
HilleRisLambers
,
J. V.
de Koppel
,
L.
Kumar
,
H. H.
Prins
, and
A. M.
de Roos
,
Am. Nat.
160
,
524
(
2002
).
https://doi.org/10.1086/342078
Google Scholar
Crossref
Search ADS
PubMed
 
10.
E.
Gilad
,
J.
von Hardenberg
,
A.
Provenzale
,
M.
Shachak
, and
E.
Meron
,
J. Theor. Biol.
244
,
680
(
2007
).
https://doi.org/10.1016/j.jtbi.2006.08.006
Google Scholar
Crossref
Search ADS
PubMed
 
11.
S.
Kéfi
,
M.
Rietkerk
,
C. L.
Alados
,
Y.
Pueyo
,
V. P.
Papanastasis
,
A.
ElAich
, and
P. C.
De Ruiter
,
Nature
449
,
213
(
2007
).
https://doi.org/10.1038/nature06111
Google Scholar
Crossref
Search ADS
PubMed
 
12.
G.-Q.
Sun
,
L.-F.
Hou
,
L.
Li
,
Z.
Jin
, and
H.
Wang
,
J. Math. Biol.
85
,
50
(
2022
).
https://doi.org/10.1007/s00285-022-01825-0
Google Scholar
Crossref
Search ADS
PubMed
 
14.
J. A.
Sherratt
,
Nonlinearity
23
,
2657
(
2010
).
https://doi.org/10.1088/0951-7715/23/10/016
Google Scholar
Crossref
Search ADS
 
15.
J. A.
Sherratt
,
SIAM J. Appl. Math.
73
,
1347
(
2013
).
https://doi.org/10.1137/120899510
Google Scholar
Crossref
Search ADS
 
16.
G.-Q.
Sun
,
H.-T.
Zhang
,
Y.-L.
Song
,
L.
Li
, and
Z.
Jin
,
J. Differ. Equ.
329
,
395
(
2022
).
https://doi.org/10.1016/j.jde.2022.05.009
Google Scholar
Crossref
Search ADS
 
18.
K. N.
Suding
,
K. L.
Gross
, and
G. R.
Houseman
,
Trends Ecol. Evol.
19
,
46
(
2004
).
https://doi.org/10.1016/j.tree.2003.10.005
Google Scholar
Crossref
Search ADS
PubMed
 
19.
C. S.
Holling
,
Annu. Rev. Ecol. Syst.
4
,
1
(
1973
).
https://doi.org/10.1146/annurev.es.04.110173.000245
Google Scholar
Crossref
Search ADS
 
20.
B.
Walker
,
C. S.
Holling
,
S. R.
Carpenter
, and
A.
Kinzig
,
Ecol. Soc.
9
,
5
(
2004
).
https://doi.org/10.5751/ES-00650-090205
Google Scholar
Crossref
Search ADS
 
21.
M.
Scheffer
,
S.
Carpenter
,
J. A.
Foley
,
C.
Folke
, and
B.
Walker
,
Nature
413
,
591
(
2001
).
https://doi.org/10.1038/35098000
Google Scholar
Crossref
Search ADS
PubMed
 
22.
M.
Scheffer
,
J.
Bascompte
,
W. A.
Brock
,
V.
Brovkin
,
S. R.
Carpenter
,
V.
Dakos
,
H.
Held
,
E. H.
Van Nes
,
M.
Rietkerk
, and
G.
Sugihara
,
Nature
461
,
53
(
2009
).
https://doi.org/10.1038/nature08227
Google Scholar
Crossref
Search ADS
PubMed
 
23.
M.
Rietkerk
,
S. C.
Dekker
,
P. C.
De Ruiter
, and
J.
van de Koppel
,
Science
305
,
1926
(
2004
).
https://doi.org/10.1126/science.1101867
Google Scholar
Crossref
Search ADS
PubMed
 
24.
K.
Siteur
,
E.
Siero
,
M. B.
Eppinga
,
J. D.
Rademacher
,
A.
Doelman
, and
M.
Rietkerk
,
Ecol. Complex.
20
,
81
(
2014
).
https://doi.org/10.1016/j.ecocom.2014.09.002
Google Scholar
Crossref
Search ADS
 
25.
R.
Bastiaansen
,
O.
Jaïbi
,
V.
Deblauwe
,
M. B.
Eppinga
,
K.
Siteur
,
E.
Siero
,
S.
Mermoz
,
A.
Bouvet
,
A.
Doelman
, and
M.
Rietkerk
,
Proc. Natl. Acad. Sci. U.S.A.
115
,
11256
(
2018
).
https://doi.org/10.1073/pnas.1804771115
Google Scholar
Crossref
Search ADS
PubMed
 
26.
R.
Bastiaansen
,
A.
Doelman
,
M. B.
Eppinga
, and
M.
Rietkerk
,
Ecol. Lett.
23
,
414
(
2020
).
https://doi.org/10.1111/ele.13449
Google Scholar
Crossref
Search ADS
PubMed
 
27.
M.
Rietkerk
,
R.
Bastiaansen
,
S.
Banerjee
,
J.
van de Koppel
,
M.
Baudena
, and
A.
Doelman
,
Science
374
,
eabj0359
(
2021
).
https://doi.org/10.1126/science.abj0359
Google Scholar
Crossref
Search ADS
PubMed
 
28.
E.
Siero
,
K.
Siteur
,
A.
Doelman
,
J. V. D.
Koppel
,
M.
Rietkerk
, and
M. B.
Eppinga
,
Am. Nat.
193
,
472
(
2019
).
https://doi.org/10.1086/701669
Google Scholar
Crossref
Search ADS
PubMed
 
29.
F. T.
Maestre
,
Y.
Le Bagousse-Pinguet
,
M.
Delgado-Baquerizo
,
D. J.
Eldridge
,
H.
Saiz
,
M.
Berdugo
,
B.
Gozalo
,
V.
Ochoa
,
E.
Guirado
,
M.
García-Gómez
 et al.,
Science
378
,
915
(
2022
).
https://doi.org/10.1126/science.abq4062
Google Scholar
Crossref
Search ADS
PubMed
 
30.
R. D.
Bardgett
,
J. M.
Bullock
,
S.
Lavorel
,
P.
Manning
,
U.
Schaffner
,
N.
Ostle
,
M.
Chomel
,
G.
Durigan
,
E. L.
Fry
,
D.
Johnson
 et al.,
Nat. Rev. Earth Environ.
2
,
720
(
2021
).
https://doi.org/10.1038/s43017-021-00207-2
Google Scholar
Crossref
Search ADS
 
31.
F.
Tröltzsch
,
Optimal Control of Partial Differential Equations: Theory, Methods and Applications
(
American Mathematical Society
,
2024
), Vol. 112.
Google Scholar
 
32.
M. R.
Garvie
 and
C.
Trenchea
,
SIAM J. Control Optim.
46
,
775
(
2007
).
https://doi.org/10.1137/050645415
Google Scholar
Crossref
Search ADS
 
33.
E.
Trélat
,
J.
Zhu
, and
E.
Zuazua
,
Math. Models Methods Appl. Sci.
28
,
1665
(
2018
).
https://doi.org/10.1142/S021820251840002X
Google Scholar
Crossref
Search ADS
 
34.
E.
Casas
,
C.
Ryll
, and
F.
Troltzsch
,
SIAM J. Control Optim.
53
,
2168
(
2015
).
https://doi.org/10.1137/140978855
Google Scholar
Crossref
Search ADS
 
35.
C.
Liu
,
S.
Gao
,
M.
Song
,
Y.
Bai
,
L.
Chang
, and
Z.
Wang
,
Chaos
32
,
063115
(
2022
).
https://doi.org/10.1063/5.0087855
Google Scholar
Crossref
Search ADS
PubMed
 
36.
S.
Gao
,
L.
Chang
,
I.
Romić
,
Z.
Wang
,
M.
Jusup
, and
P.
Holme
,
J. R. Soc. Interface
19
,
20210739
(
2022
).
https://doi.org/10.1098/rsif.2021.0739
Google Scholar
Crossref
Search ADS
PubMed
 
37.
L.
Chang
,
W.
Gong
,
Z.
Jin
, and
G.-Q.
Sun
,
SIAM J. Appl. Math.
82
,
1764
(
2022
).
https://doi.org/10.1137/22M1472127
Google Scholar
Crossref
Search ADS
 
38.
J.
Liang
 and
G.-Q.
Sun
,
Nonlinear Dyn.
112
,
8675
(
2024
).
https://doi.org/10.1007/s11071-024-09500-3
Google Scholar
Crossref
Search ADS
 
39.
L.
Chang
,
X.
Wang
,
G.
Sun
,
Z.
Wang
, and
Z.
Jin
,
J. Math. Biol.
88
,
5
(
2024
).
https://doi.org/10.1007/s00285-023-02026-z
Google Scholar
Crossref
Search ADS
 
40.
L.-F.
Hou
,
G.-Q.
Sun
, and
M.
Perc
,
Commun. Nonlinear Sci. Numer. Simul.
126
,
107461
(
2023
).
https://doi.org/10.1016/j.cnsns.2023.107461
Google Scholar
Crossref
Search ADS
 
41.
G.-Q.
Sun
,
C.-H.
Wang
,
L.-L.
Chang
,
Y.-P.
Wu
,
L.
Li
, and
Z.
Jin
,
Appl. Math. Model.
61
,
200
(
2018
).
https://doi.org/10.1016/j.apm.2018.04.010
Google Scholar
Crossref
Search ADS
 
42.
Y.
Mau
,
L.
Haim
, and
E.
Meron
,
Phys. Rev. E
91
,
012903
(
2015
).
https://doi.org/10.1103/PhysRevE.91.012903
Google Scholar
Crossref
Search ADS
 
43.
V.
Barbu
,
Analysis and Control of Nonlinear Infinite Dimensional Systems
, Mathematics in Science and Engineering (
Academic Press
,
1993
), Vol. 190.
Google Scholar
 
© 2024 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).