Abstract

In this paper, a decision-making competition game model concerning governments, agricultural enterprises, and the public, all of which participate in the reduction of nitrogen emissions in the watersheds, is established based on bounded rationality. First, the stability conditions of the equilibrium points in the system are discussed, and the stable region of the Nash equilibrium is determined. Then, the bifurcation diagram, maximal Lyapunov exponent, strange attractor, and sensitive dependence on the initial conditions are shown through numerical simulations. The research shows that the adjustment speed of three players’ decisions may alter the stability of the Nash equilibrium point and lead to chaos in the system. Among these decisions, a government’s decision has the largest effect on the system. In addition, we find that some parameters will affect the stability of the system; when the parameters become beneficial for enterprises to reduce nitrogen emissions, the increase in the parameters can help control the chaotic market. Finally, the delay feedback control method is used to successfully control the chaos in the system and stabilize it at the Nash equilibrium point. The research of this paper is of great significance to the environmental governance decisions and nitrogen reduction management.

1. Introduction

According to the European Nitrogen Assessment, the total economic loss caused by the reactive nitrogen in 27 countries of the European Union amount to 70–320 billion euros per year. The economic cost of the reactive nitrogen pollution is about twice that of Europe’s “willingness” to pay for carbon controls. To integrate the research on global nitrogen emissions and nitrogen pollution, the European Union puts forward the “nitrogen and Europe” research plan, and all countries of the world were invited to participate. Due to the intensification of nonpoint source nitrogen pollution in agriculture enterprises, the problem of nitrogen pollution in river basins is becoming more and more serious. Two-thirds of the coastal rivers and bays in the United States are degraded from nutrient pollution, and nitrogen inputs in these waters continue to increase [1]. Therefore, controlling the input of nitrogen and phosphorus from human activities is essential in reducing eutrophication in watersheds [2]. As the reduction of nitrogen emission starts to attract the attention of the whole world, the nitrogen emission trading market and the nitrogen emission limits of various industries are gradually being formed, and agricultural enterprises will face this major challenge. At the same time, governments should not only consider the environmental benefits but also control the normal operation of the whole market. Therefore, as the main players in the reduction of nitrogen emissions in the river basin, any party in government and agricultural enterprises and the social public decisions will be influenced by the other two parties.

Some literature mainly focuses on watershed pollution based on game theory. Initially, there has been much discussion on the treatment of environmental problems. Ni and Wang [3] used a cooperative game to analyze the allocation of pollution control costs in watershed pollution and explored a reasonable allocation method. Gao et al. [4] analyzed the interaction among upstream governments, downstream governments, and the central government in the Eastern Route of South-to-North Water Transfer Project based on evolutionary game theory. Secondly, there is cross-border water pollution in the river basin. Jorgensen [5] took the upstream and downstream areas in the river basin as the main body of the game and analyzed whether the cooperation between the upstream and downstream could solve the problem of unreasonable pollution through the differential game method. Frisvold and Caswell [6] used the static bargaining game method to study the impact of pollution control policies on the game relationship between two countries in terms of environmental pollution control. Third, there is a conflict of interest between subjects in the river basin. Bárcena-Ruiz [7] used the idea of a differential game to analyze whether two governments should solve the problem of river basin pollution by setting environmental taxes to be the same.

The majority of nitrogen emission reduction in the watershed is based on bounded rationality. At present, the research on dynamic competitive game with bounded rationality comes mainly from the following authors. Puu [8] first found a variety of complex dynamic phenomena in the Cournot duopoly model such as the singular attractor with a fractal dimension. Yali [9] studied a delayed duopoly game considering increasing marginal costs based on bounded rationality and demonstrated that state delay is helpful in enlarging the stability region of the system. Peng et al. [10] and Elsadany [11] discussed the correction of a duopoly game with bounded rationality based on the strategy of maximizing the output expectations of enterprises. Yao and Xu [12] established an advertising market competition model that considered the bounded rationality of participants and analyzed the complex decision-making behaviors of decision-makers in the dynamic game process. Ding et al. [13] proposed a linear dynamic system in a duopoly game involving renewable resource extraction with the strategy of bounded rationality. Yao et al. [14] and Elabbasy et al. [15] both constructed a nonlinear triopoly game model with heterogeneous players, and the three different decision-makers were bounded rational, adaptive, and naive. Research by Zhao [16] investigated a novel Cournot duopoly game model of carbon emission reduction based on the hypothesis of participant’s bounded rationality.

Our study is closely related to reduction mechanism, which can be divided into mandatory emission reduction mechanism [17, 18] and incentive emission reduction mechanism [19, 20]. Wang et al. [21] analyzed the relationship between supply chain enterprise operation and government policy. A three-stage Stackelberg game model of decentralized supply chain and a two-stage Stackelberg game model of centralized supply chain were used to study the government’s carbon emission tax policy. De Jonge [22] proposed that the instruments of legislation, subsidies, green taxation, and emission trading can help achieve reduction targets for NOx. Research by Svensson and Elofsson [23] showed that the net nitrogen reductions achieved through environmental policy efforts and the costs of the nitrogen reductions should be considered. According to previous research, most scholars mainly focused on the decision-making game of emission reduction, price, or output of enterprises, whereas seldom discussed the decisions of government policies and the public. However, when the government and the public are involved in nitrogen emission reduction work, under a series of environmental policies and public supervision, studying the complex dynamic behavior of a game involving governments, agricultural enterprises, and the public with bounded rationality will have great practical significance.

The contributions of this study are as follows: first, a reasonable decision-making mechanism for nitrogen emission reduction is proposed considering the influence of government and the public decisions in emission reduction. Second, how to establish a reasonable adjustment strategy of output, supervision intensity, and policy intensity is explored. Third, the integration of nonlinear dynamic theory and nitrogen emission reduction management is fulfilled. Moreover, the effects of the market price of nitrogen trading and the subsidy standard of nitrogen emission reductions on the decisions of three players are simulated, and it is important to adjust the parameters that will become beneficial for enterprises to reduce nitrogen emissions. This study can provide a theoretical guidance for reducing nitrogen pollution in the watersheds.

The organization of the paper is as follows: in Section 2, a competition model concerning governments, agricultural enterprises, and the public is formulated; in Section 3, we analyze the equilibrium point of the game model and give the existence condition and local stability range of the equilibrium point. In Section 4, the complex dynamic behavior under a change in the adjustment speed of three players’ decisions is analyzed. In Section 5, we apply the delay feedback method to control chaos in the system. Finally, some research conclusions are summarized in Section 6.

2. Model

This main purpose of this chapter was to introduce the aforementioned dynamic game model. Considering nitrogen emission trading, nitrogen emission reduction subsidies, and marginal emission reduction costs, this paper analyzes whether governments, agricultural enterprises, and the public make an optimal decision according to their own decision rules in the game. Lanoie et al. [24] discussed the impact of environmental policies on environmental innovation performance, and the intensity of environmental policies was divided into three categories: weak, narrow, and strong. Thus, the optimal goal of the government is to choose the appropriate policy intensity to maximize the total utility, and the optimal goal of the agricultural enterprise is to make the appropriate output decision to maximize its profit when the pollution situation meets the government’s policy intensity. In terms of the public, the best goal is to choose the appropriate supervision intensity to maximize its total utility. The meanings of specific parameters and variables are shown in Table 1.

The following assumptions are made to develop the model:(1)This paper is mainly aimed at agricultural enterprises, and the price of the enterprises in the period t is determined by through the inverse demand function , where a and are positive constants; the production cost of the enterprise is a linear function, namely, .The nitrogen emissions generated in the production process of an enterprise are linearly related to its output, denoted as . The emission reduction of enterprises is related to their own technical level of nitrogen emission reduction, the government policy intensity, and the public’s supervision intensity. Wu et al. [25] proposed that public participation had significant positive effects on the reduction of both binding and nonbinding environmental pollutant emissions. Therefore, the emission reduction of enterprises is [21, 26] as follows:where (i = 1, 2, 3); thus, the final nitrogen emission of enterprises is . The cost of nitrogen emission reduction is . When participating in nitrogen emission trading, the tradable emission permits that enterprises need can be described as ; then, the fee for nitrogen emission trading is .(2)In terms of governments, the revenue function of the government in this paper mainly includes four parts: tax, nitrogen pollution treatment costs, supervision cost, and subsidy expense. According to Wang et al. [21] and Alexeev et al. [26], the nitrogen pollution treatment costs are . The government provides subsidies and incentives for enterprises to reduce nitrogen emissions, so the subsidy expense can be described as .(3)Regarding the public, the research of Carreira et al. [27] showed that the degree of public participation in corporate environmental behavior depended on government’s policy intensity. The supervision cost of the public on enterprises is negatively correlated with the government’s policy intensity [28], so we assume the supervision cost is . Based on Newig et al. [29], the utility function of the consumer can be described as

Therefore, the profit function of a government, agricultural enterprise, and the public is

Then, the marginal profit of a government, enterprise, and the public in period t is

Due to incomplete market information and limitations of their own conditions in reality, when governments, agricultural enterprises, and the public make decisions with bounded rationality, they cannot fully predict the future market demand. Therefore, it is assumed that they can only determine their decisions based on the local estimation of marginal profit. If the marginal profit is positive in period t, they will increase their decision quantity in period t+1. Thus, a three-dimensional discrete dynamic game model in the t+1 period is set up as follows [30]:where (i = 1, 2, and 3), respectively, represents the adjustment speed of each bounded rational player. For the convenience of calculation, we assume , , and , so the dynamic adjustment mechanism of the government, enterprises, and the public with bounded rationality is simplified as

3. Equilibrium Points and Local Stability

In order to study the dynamic behavior of the game model, the nonnegative equilibrium point will be discussed in this chapter. In system (6), equilibrium points are obtained by setting x(t + 1) = x(t), y(t + 1) = y(t), and z(t + 1) = z(t), so we can obtain six equilibrium points:

Obviously, , and are bounded equilibrium points. When , is a Nash equilibrium point. To discuss the local stability of the above equilibrium points, we must consider the Jacobian matrix of system (6):

Theorem 1. If the Nash equilibrium point is strictly nonnegative, the boundary equilibrium points , and of system (6) are unstable equilibrium points.

Proof. In order to prove this result, we find the eigenvalues of the Jacobian matrix at each boundary equilibria , and . The Jacobian matrix at iswhose eigenvalues are , and . Since , is the government’s marginal supervisory cost, and represents the increase in environmental benefits with an increase in , so we can get . Otherwise, the government’s policy improvement will not make any sense, namely, . Since , is the highest price of the product in the market, while can be regarded as the total variable costs of the enterprise. In the actual market, the highest price of the product must be higher than its total variable costs, namely, . With conditions where and , can be obtained. Therefore, is an unstable equilibrium point.
The Jacobian matrix at iswhose eigenvalues are , . Since , it is clear that when the condition , is obtained. Then, is an unstable equilibrium point.
The Jacobian matrix at iswhose eigenvalues are , , and . It is clear that when the condition , . Then, is an unstable equilibrium point. Similarly, we can prove that and are also unstable.
The Jacobian matrix at isBy calculating the eigenvalue of the Jacobian matrix , we can find . Thus, is an unstable equilibrium point.
The Jacobian matrix at isBy calculating the eigenvalue of the Jacobian matrix , . It is clear that when the condition , . Thus, is an unstable equilibrium point.

Theorem 2. If the system parameters satisfy , and when the following Jury conditions are performed, the Nash equilibrium point is locally asymptotically stable.

Proof. In order to investigate the local stability of the Nash equilibrium point , the Jacobian matrix at isThe characteristic equation of the matrix iswhereThe local stability conditions of the Nash equilibrium are given by Jury’s conditions, which are the sufficient and necessary conditions :Obviously, the Nash equilibrium point is a stable node in the stability region defined by (17). However, if , and go beyond the stability region, more complex phenomena in terms of the evolution of outputs will occur such as bifurcation and chaos. Moreover, we found that the local stability of the system in the Nash equilibrium point can be decided by every parameter in (17). Based on inequalities (17), the three-dimensional stability domains of the system (6) are simulated when and take different values (as shown in Sections 4.2 and 4.3).

4. Numerical Simulations

In this section, we analyzed the dynamic behaviors of the bounded rational players through various numerical simulations. They could observe the influence of the adjustment speed of , and , the market price of nitrogen trading , and the subsidy standard of nitrogen emission reduction on the model. In order to study the local stability properties of the equilibrium point, it is convenient to take the parameter values as follows: , , , , , , .

4.1. The Impact of the Adjustment Speed on the Stability of the System

Figure 1 shows the bifurcation diagram with respect to the adjustment speed of a government’s policy intensity while  = 0.21 and  = 0.5. The corresponding largest Lyapunov exponents with respect to are drawn in Figure 2. In the range , the Lyapunov exponents are negative, which means that the Nash equilibrium point is stable. When  = 1.1057, the first bifurcation point in Figure 1 corresponds to the first peak (1.1057, −0.008) in Figure 2, leading to the system gradually entering a period-doubling bifurcation. Finally, when , the maximal Lyapunov exponents are almost greater than zero, indicating that chaotic behavior is occurring and the Nash equilibrium point is becoming very unstable.

Figure 3 shows the bifurcation diagram with respect to the adjustment speed of the enterprises’ output while  = 0.7 and  = 0.5. The corresponding largest Lyapunov exponents with respect to are drawn in Figure 4. As can be seen from Figures 3 and 4, in the range , the Lyapunov exponents are negative, which means that the Nash equilibrium point is stable. When  = 0.2279, the first bifurcation point in Figure 3 corresponds to the first peak (0.2279, −0.0282) in Figure 4. With increasing to 0.2792, the second bifurcation point in Figure 3 corresponds to the second peak (0.2792, −0.0473) in Figure 4, and the system then gradually enters a period-doubling bifurcation. Finally, when , the maximal Lyapunov exponents are almost greater than zero, indicating that chaotic behavior is occurring and the Nash equilibrium point is becoming very unstable.

Figure 5 shows the bifurcation diagram with respect to the adjustment speed of the public’s supervision intensity while  = 0.7 and  = 0.21. The corresponding largest Lyapunov exponents with respect to are drawn in Figure 6. From Figures 5 and 6, when , the Lyapunov exponents are negative, which means that the Nash equilibrium point is locally stable for small values of . However, in the range , the system starts to enter into the chaotic state and complex dynamic behavior occurs.

The strange attractors corresponding to Figure 1 are shown in Figures 710, which shows the changing situation for strange attractors at different values of while  = 0.21 and  = 0.5. When  = 0.008, the decision-making behavior of the government, enterprises, and the public forms a spiral trajectory map and finally forms a gradual stability point. However, with the finiteness of market information and the bounded rationality of the game players, when , the stable point gradually appears as a branch state, as shown in Figure 8. When , it was found that the point was no longer stable and chaos began to appear until a chaos phenomenon in Figure 9 appeared.

Figure 10 shows the strange attractor in a chaotic state while  = 0.7,  = 0.31, and  = 0.5. At this time, the decision-making behavior of players appears to be a complex chaos phenomenon. The strange attractors corresponding to Figure 5 are shown in Figures 11 and 12, which shows the change situation for strange attractors at different values of while  = 0.7 and  = 0.21. When increases to 8.5, chaos has occurred and a vortex shaped attractor appears, as shown in Figure 11. When , it was found that the vortex had evolved into a annular phase diagram, as shown in Figure 12.

In order to further explore the chaotic phenomenon caused by a change in the decision-making adjustment speed, we investigated the sensitivity at the initial value of the system (6). These numerical simulations are performed by setting  = 1.28,  = 0.21, and  = 0.5 (the system is in a chaotic state at this time). It can be seen from Figure 13 that two orbits of x(t), y(t), or z(t) are indistinguishable at the beginning, but after several iterations, the separation between them builds up rapidly; that is, subtle changes in the initial conditions will greatly affect the results.

Through the above numerical simulation analysis, it can be concluded that the adjustment speeds , , and of the bounded rational players may greatly affect the stability of system (6) and lead to complex chaos phenomena in the system. Once trapped in a chaotic market, slight changes in various initial conditions of the government, agricultural enterprises, and the public will greatly affect the final results. In addition, the players cannot effectively predict various changes in reality, which will result in their decisions not being effectively implemented.

4.2. The Impact of the Market Price of Nitrogen Trading on the Stability of the System

When enterprises decide whether to trade emission permits based on their own nitrogen emissions, it is necessary to compare the market price of nitrogen trading with the cost of nitrogen emission reductions. The emission level and emission reduction technology level of enterprises are determined by their production equipment and technology, which cannot be changed quickly. The variable cost of nitrogen emissions can affect the decision-making behavior of enterprises; therefore, it is the market price of nitrogen trading that affects the stability of the system.

Inequalities (17) define the stable range of the Nash equilibrium point of the system under the adjustment speeds , , and . When the initial values of each parameter are fixed, the region of stability for the Nash equilibrium point under different values of is shown as in Figure 14. When the trading price of nitrogen emissions  = 1 increases to  = 10, the stability of the system will decrease. In addition, the Nash equilibrium point will evolve from (1.0545, 1.4680, 1.2720) to (1.0545, 1.3180, 1.3470), indicating that, with an increase in , the government’s policy intensity will remain unchanged, while the enterprises’ output will be reduced and the public’s supervision intensity will increase.

Figure 15(a) shows the bifurcation diagram with respect to the market price of nitrogen trading while  = 0.7,  = 0.21, and  = 0.5 (the system is stable). The Nash equilibrium point then becomes . From this Figure 15(a), it can be observed that the equilibrium point is locally stable for the small values of the parameter . When increases, the Nash equilibrium point becomes unstable, and even complex dynamics phenomena such as period-doubling bifurcation and chaos appear. The main reason behind this is that an increase in increases the variable cost for enterprises. When the technical level remains unchanged, enterprises have to reduce their output.

Figure 15(b) shows the bifurcation diagram with respect to while  = 1.25,  = 0.21, and  = 0.5 (the system being in chaos). The research shows that when the adjustment speed of policy intensity is too large, the public’s supervision intensity will decrease to zero with an increase in , and the system still remains in a chaotic state. Figure 15(c) shows the bifurcation diagram with respect to while  = 0.7,  = 0.315, and  = 0.5 (the system being in chaos). It can be seen that when the adjustment speed of enterprises’ output is too large, as increases, the system gradually evolves from chaos to period-doubling bifurcation until reaching a state of equilibrium. However, when continues to increase, a complex evolution similar to that shown in Figure 15(a) will appear and eventually enter chaos. Figure 15(d) shows the bifurcation diagram with respect to while  = 0.7,  = 0.21, and  = 10 (the system being in chaos). When the adjustment speed of the public’s supervision intensity is too large, each player is still in chaotic state with an increase in . The main reason for this is that has an influence on the decision of enterprises’ output y(t) by affecting the marginal profit of enterprises. Therefore, when is too large, increasing can control chaos.

4.3. The Impact of Subsidizing Nitrogen Emission on the Stability of the System

The subsidy of nitrogen emission reductions is proportional to the amount of emission reductions. The subsidy of emission reductions is not only a source of income for an enterprise but also the government’s fiscal expenditure to encourage enterprises to reduce emissions. Therefore, the subsidy standard of the reduction of nitrogen emissions will affect the system.

When the subsidy standard of the reduction of nitrogen emissions  = 0.3, the stable region of the Nash equilibrium point is shown as in Figure 14(a). If other parameters are fixed, the nitrogen emission reduction subsidy standard varies to  = 0.32 from  = 0.3, and we can see that the area of the stable region increases in the direction of and , as shown in Figure 16. Therefore, the stability of the system will increase with an increase in ; in addition, the Nash equilibrium point will evolve from (1.0545, 1.468, 1.272) to (1.0545, 1.4698, 0.7503), which means that, as increases, the government’s policy intensity will remain unchanged, while the enterprises’ output will increase and the public’s supervision intensity will decrease.

Figure 17(a) shows the bifurcation diagram with regard to the subsidy standard of the nitrogen emission reduction while  = 0.7,  = 0.21, and  = 0.5 (the system being stable) because the Nash equilibrium point at this time becomes  = (1.0545, 1.4417 + , ( − 7.0481 − )). From Figure 18, it can be observed that the equilibrium point is locally stable for small values of the parameter . When increases, the Nash equilibrium point becomes unstable, and even complex dynamic phenomena appear such as period-doubling bifurcation and chaos. The main reason is that the increase in increases the government’s variable cost, and the government has to reduce its policy intensity. When other conditions remain unchanged, enterprises can obtain more subsidies by increasing their output. At this time, the public will reduce their supervision intensity due to the increase in enterprises’ emission reduction. When becomes too large, the market cannot be balanced, and the decision-making of the government, enterprises, and the public cannot reach an equilibrium point any more until chaos appears.

Figure 17(b) shows the bifurcation diagram with regard to while  = 1.25,  = 0.21, and  = 0.5 (the system being in chaos). The research shows that when the adjustment speed of policy intensity is too large, as increases, the public’s supervision intensity will decrease to zero, and the system gradually evolves from chaos to period-doubling bifurcation until reaching an equilibrium state. However, when continues to increase, a complex evolution similar to that shown in Figure 17(a) will appear and eventually enter chaos. Figure 17(c) shows the bifurcation diagram with respect to while  = 0.7,  = 0.315, and  = 0.5 (the system being in chaos). It can be seen that when the adjustment speed of enterprises’ output is too large, as increases, the government and enterprises are still in a chaotic state. Figures 17(d) and 17(e) show the bifurcation diagram with respect to while  = 0.7,  = 0.21, and  = 10 (the system being in chaos). When the adjustment speed of the public’s supervision intensity is too large, the system gradually evolves from chaos to an equilibrium state with an increase in . Therefore, when or is too large, increasing can control chaos.

The above numerical simulation shows that the market price of nitrogen trading and the subsidy standard for the reduction of nitrogen emission are important factors in the dynamic game among governments, enterprises, and the public participating in the reduction of nitrogen emissions in the basin. They not only influence the Nash equilibrium point of the system but also affect the stable region of the system.

5. Chaos Control

Through model analysis and numerical simulation, it is found that when , , or exceeds the critical value, the system (6) will lose stability. At this time, the chaotic system will have a sensitive dependence on the initial conditions, which means that the government, enterprises, and the public would not be able to predict the market development and any small adjustment of the initial conditions. Therefore, it is very important to perform chaotic control on the system (6) to ensure that it is in a stable equilibrium state.

There are many chaos control methods. This section uses the delayed feedback control method proposed by Pyragas [31] to control the chaos of the system (6). It is expressed as , , where is the controlling factor and is the length of the time delay. Substituting  = 1 into the second equation of the system (6), the controlled system can be modeled asand the Jacobian matrix of (17) at the Nash equilibrium point is

Figure 19 shows the bifurcation diagram with regard to the control factor k, while the initial values of the other parameters are fixed, and  = 0.7,  = 0.315, and  = 0.5. Figure 18 shows the largest Lyapunov exponents with regard to the control factor k. From Figure 19, with the increase in k, the decision variables x(t), y(t), and z(t) can evolve from chaos to periodic bifurcation and finally stabilize at the Nash equilibrium levels. With a gradual increase in k, in the range k > 0.384, the controlled system (18) becomes stable without chaotic behaviors. The effects of the control factor k on the controlled system before and after chaos are shown in Figure 20 when k = 0.45, and this figure depicts the change process of the controlled system from chaos to a stable state when the initial values of the bounded rational players are  = (0.7, 0.4, 0.3).

6. Conclusions

In this paper, bounded rationality, nitrogen emission trading, and the subsidy of reductions in nitrogen emissions are considered in terms of a dynamic game involving the government, enterprises, and the public, and a decision-making game model is established based on bounded rationality. At the same time, we analyzed the dynamic behavior of players with bounded rationality, the equilibrium points of the model are discussed, and a three-dimensional stability region of the Nash equilibrium point is presented. Through the discussion, it can be concluded that many parameters such as the market price of nitrogen trading and the subsidy standard of nitrogen emission reductions would affect the stability of the system; when the parameters become beneficial for enterprises to reduce nitrogen emissions, the chaotic market will restore, and the regional stability of the system will decrease with the increase in the parameters. Furthermore, the numerical simulation shows the dynamic evolution process of the decisions of the participants. The results show that when the adjustment speed values of the bounded rational player , , and are small, the system is stable. If one of , , and increases beyond the stability region of the Nash equilibrium point, bifurcation, chaos, and other dynamic behaviors will occur. Finally, it is proven that the delayed feedback control method can effectively control the system in a chaos state to restore the stable equilibrium market.

Data Availability

Some or all data, models, or code generated or used during the study are available from the corresponding author by request.

Disclosure

Jixiang Zhang and Xuan Xi should be regarded as the co-first authors.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Authors’ Contributions

Jixiang Zhang and Xuan Xi contributed to the work equally.

Acknowledgments

This work was supported by the Fundamental Research Funds for the Central Universities (Grant no. NS2019045).