A conceptual model for simulating responses of freshwater macroinvertebrate assemblages to multiple stressors

Simulating macroinvertebrate responses to multiple environmental stressors is an important tool for water quality management, by predicting ecological effects of both stressors and restoration practices. Currently, existing modelling approaches fall short in simulating the responses of macroinvertebrate assemblages to environmental constraints, lacking incorporation of the multiple spatial and temporal scales on which stressors act, including their mutual interactions and uncertainties associated with input data. In answer to these shortcomings, this study aimed to design a conceptual multiscale model for simulating responses of macroinvertebrate assemblages to multiple environmental stressors. To this purpose, we drew up model requirements, selected model building blocks and assembled these into a conceptual model, also documenting the challenges that remain to be solved. This conceptual model offers a direction for simulating responses of macroinvertebrate assemblages to multiple stressors, which in turn can be used to better focus management resources and restoration practices.


Introduction
Simulating macroinvertebrate responses to multiple environmental stressors is an important tool in water quality management, for verifying process understanding, for predicting ecological effects of increased stress as well as restoration practices, and for identifying the regions on which management resources should be focused Pace, 2001). Currently available models for simulating macroinvertebrate responses to multiple stressors are however limited in 1) representing environmental constraints acting on multiple spatial and temporal scales, 2) including the mutual interactions between stressors, and 3) considering the margins of uncertainty associated with the input data. Consequently, these models may not help to gain insight in the reasons for the observed community structure and absence and presence of specific macroinvertebrate species, which in turn does not aid to focus management resources and restoration practices. Below we will discuss these shortcomings in more detail.
1) Environmental constraints acting on multiple scales: The scales on which stressors affect macroinvertebrate assemblages range from the catchment down to local habitats (Allan et al., 1997;Fitzpatrick et al., 2001;Frissell et al., 1986). Yet, to simulate macroinvertebrate responses, current models tend to focus mainly on habitat scale parameters, such as substrate types, nutrient concentrations and water depth (Bennetsen et al., 2016;Gobeyn et al., 2017) and to a lesser extent on the full array of habitat to catchment scale parameters such as discharge regime, land use type and connectivity, as well as oxygen regime and nutrient concentrations as is done in Boets et al., 2015;Kail et al., 2015;Elias et al., 2016 andStoll et al., 2016. In addition, stressors originating from traditionally separated scientific disciplines, like hydrology, chemistry, morphology and system conditions (the large-scale long-term conditional characteristics like climate, geology and geomorphology), are often considered in isolation, since model types maintain a focus on one or a few disciplines, e.g., chemistry is often represented well in processbased models (Poepperl, 2003;Traas et al., 1998), while morphology is better represented in statistical habitat models (Garbe et al., 2016;Guse et al., 2015;Tomsic et al., 2007). Hence, such approaches ignore the combination of large-scale pressures and local scale ecological stressors that jointly affect macroinvertebrate assemblage composition. Therefore, there is a need for models that account for the entire range of spatial scales on which stressors act. 2) Stressor interactions: The multiple stressors affecting macroinvertebrate assemblages are often treated in an additive way. However, in reality, these stressors may interact and a combination of multiple stressors may therefore increase the experienced stress (synergistic interaction) or decrease the experienced stress (antagonistic interaction) (Coors and de Meester, 2008;Jackson et al., 2016). An additive approach obviously does not account for these interactions that may change the response of macroinvertebrate assemblages to multiple stressors. 3) Uncertainty: The uncertainty in the simulated responses associated with the use of expert knowledge as input data is often not communicated well, which can lead to over-interpretation of the model results (Refsgaard and Henriksen, 2004). Yet, expert knowledge may sometimes be the only option, since data on stressor-response relations is often limited. Expert knowledge can also be used in combination with data-driven techniques to quantify relationships that are not covered by the data (Adriaenssens et al., 2004;Boets et al., 2015;Skeffington et al., 2014), although it may be tedious to formalize the knowledge of experts into a workable format, and problems with consistency may arise (Mouton et al., 2009).
To address the abovementioned shortcomings of the current available models for simulating macroinvertebrate responses to multiple stressors, two main types of models can be used: statistical and processbased models (Janssen et al., 2015). In view of the intended application, both model types have advantages and limitations, which are inherent to their design. Statistical models are based on empirical relationships between stressors and responses and may be easier to use, but care should be taken with extrapolation of the derived stressorresponse relationships, which might not be valid outside the domain covered by the data (Guisan and Zimmermann, 2000;Rose et al., 2015). In addition, complex relationships and interactions are not necessarily captured by such models. Process-based models on the other hand offer a way to extrapolate results beyond the range of input data, because of their mechanistic nature. They picture the complex ecological processes in a detailed way, but are therefore more hungry for input data (Buckley et al., 2010;Guisan and Thuiller, 2005). Hence, both model types have characteristics that can be used, and possibly a combination of these model types may offer opportunities for improved macroinvertebrate response modelling.
Some studies indeed combined several model types, addressing relevant stressors on multiple scales for simulating organism responses Morales et al., 2006). Moreover, these studies did report on the uncertainties associated with the model outcomes or proposed extensions of the model to do so. Yet, also in these studies, not all of the abovementioned shortcomings were dealt with, for example, interactions between stressors were still not included. Also, these approaches cannot easily be applied in cases with limited data availability. Therefore, currently, still no ideal model is available for simulating responses of macroinvertebrate assemblages that takes into account interacting stressors acting over multiple scales, that considers the associated uncertainty and that can be applied to a variety of situations with varying data availability. Hence, this study aimed to design a conceptual multiscale model for simulating responses of macroinvertebrate assemblages to multiple environmental stressors. To this purpose, requirements for a multiscale model for simulating macroinvertebrate responses to multiple stressors are proposed, followed by a selection of available building blocks that meet these requirements. Next, a conceptual model is presented, consisting of these selected building blocks, and finally, challenges are documented that remain to be solved before an actual model that answers the requirements posed can be build. This way we may bring the modelling of responses of macroinvertebrate assemblages to multiple stressors further, which in turn will aid the focus of management resources and restoration practices.

Requirements and building blocks for an ideal macroinvertebrate response model
In answer to the major challenges described above, this section highlights the requirements for the design of an ideal multiscale, multifactorial model for simulating responses of macroinvertebrate assemblages to multiple stressors. In addition, other requirements that were already included in previous models are mentioned again here as essential parts of an ideal macroinvertebrate response model, which include the regional species pool, dispersal and biotic constraints.
Together with the listed requirements, selected building blocks are described that are suitable for constructing a multiscale model for simulating responses of macroinvertebrate assemblages to multiple stressors. These building blocks originate from a set of currently available methods.

Water body type-and region-specific analysis
Macroinvertebrate assemblage composition differs between regions and locally between water body types that differ in e.g. size of the stream, stream order and slope (Bailey, 1983;Johnson and Hering, 2009;Lake et al., 2007). Hence, the potential species pool is a regional characteristic which should be included in the model in order to obtain a realistic simulation of the composition of local macroinvertebrate assemblages. Therefore, a model frame that can be fed by water body type and region-specific data is needed.
Various frameworks for macroinvertebrate response modelling have been described that can be applied to multiple areas by training them with local data. Previous models that describe a general framework which can be adjusted with localized input are the Mussel Dynamics Model (Morales et al., 2006), the Habitat Evaluation Tool  and the eco-hydrological model cascade . These approaches used a combination of model types, like mechanistic hydrological simulations combined with statistical habitat suitability calculations, showing flexibility in the application area of such a model.
Based on the flexibility of these existing methods, we propose to use a general modelling framework, which can be adjusted to be used locally, in different water body types and in different regions. To design a model for a specific water body type and region, target assemblages should be formulated for these specific conditions. These may follow from regional typologies of macroinvertebrate assemblages (Martin and Brunke, 2012;Metcalfe, 1989;Palmer et al., 2005;. If possible, these typologies are preferably based on the distributions of species and not on arbitrarily chosen, environment-or geography-based boundaries, to which species distributions may not adhere (Lorenz et al., 2004;Magnusson and Gering, 2004;de Vries et al., 2020). Following such specific typologies, adjustments of the general framework can be made in the choice for relevant environmental stressors as prediction variables for the responses of the target species assemblage. For example, the response of a lowland stream assemblage is strongly depending on flow velocity, whereas this is of less importance to simulate the response of a lake assemblage.

Dispersal and biotic constraints
Dispersal and biotic constraints that form the local species assemblage from the regional species pool depend on the species capacities for colonisation and competition (Lake et al., 2007). These capacities are determined by all evolutionary adaptations that species developed. Dispersal constraints determine the abundance and distribution of macroinvertebrates on newly available locations. In addition, biotic constraints influence the assemblage locally, for example by competition for space and resources (Leibold et al., 2004).
There are several model types that represent dispersal and biotic constraints, although mostly separately. For a focus on biotic constraints, such as trophic-and non-trophic interactions, food web models can be used. Dispersal constraints are addressed for example in spatially explicit species distribution models (Radinger et al., 2014). However, for modelling macroinvertebrate dispersal constraints, these models may not be easily applied, as knowledge on dispersal abilities and distances is limited. Yet, dispersal abilities may be deduced from the presence of certain traits (De Bie et al., 2012;Poff, 1997;Tonkin et al., 2014).
In the statistical habitat model type, biotic constraints are normally not included (Hirzel & Le Lay 2008). However, a few studies on habitat models have shown that biotic constraints can also be partly included, in a simplified way. Dispersal constraints can be accounted for in habitat models by adding a general geographical dispersal limitation term to model dispersal from source pools (Bennetsen et al., 2016) and by considering species pools and species-specific (re)colonisation potential in the presence of migration barriers . Biotic constraints can be included indirectly by accounting for 'biotic resistance' to biological invasions of sites with a high macroinvertebrate diversity . In simulations for single species that are influenced by the presence of another specific organism, the distribution map of this organism can be included in the models as a predictor (Hirzel & Le Lay 2008).
The existing approaches described above provide building blocks to include dispersal and biotic constraints in a conceptual model for simulating macroinvertebrate responses to multiple stressors.

Environmental constraints acting on multiple scales
The environmental constraints, filtering species out of the regional species pool due to the suitability of the regional environment and the local habitat (Poff, 1997), are often represented by a limited number of environmental variables. Although it remains often unclear what the mechanistic pathways are that make an environmental variable relevant for macroinvertebrate presence and abundance, there is general agreement on a set of variables acting over multiple spatial and temporal scales that together cover most of the relevant environmental drivers for macroinvertebrate assemblages (Feld and Hering, 2007;Friberg, 2010;Verberk et al., 2012;Verdonschot et al., 1998;Villeneuve et al., 2018). In order to simulate macroinvertebrate responses to this set of variables acting on multiple scales, these should be covered in the model simulations, with a meaningful expression format.

Statistical models
Statistical models seem to be an appropriate model type to address environmental constraints acting on multiple scales. Within the statistical model type, habitat models focus on picturing the environmental constraints for a local species assemblage. These models derive biotic responses from the organism's environmental preferences and regional boundary conditions that determine habitat suitability, making it possible to include stressors acting over multiple spatial and temporal scales. This model type assumes that all suitable habitat will eventually be occupied by organisms, and therefore, these models do not directly simulate species abundance and diversity.
Statistical models are not limited to any constraint type, since predictors can be chosen to represent both direct and indirect effects of any environmental stressor (Guisan and Zimmermann, 2000). These statistical models are thus not limited in accounting for either biotic, dispersal or environmental constraints driving the local species assemblage, and therefore offer a suitable means to fulfil our requirements. Below we elaborate further on the application of statistical models in the conceptual model for simulating stress responses of macroinvertebrate assemblages.
Alternatively, environmental constraints can be modelled in more detail by using a combination of process-based models and a habitat model. Independent process-based modules can be used to simulate dynamics in an environmental category (e.g. hydraulics, chemistry), after which the outcomes are fed into a model which determines habitat suitability from the simulated environmental conditions Tomsic et al., 2007).

Expression format of environmental constraints in statistical models
Currently, improvements can be made concerning the choice and expression format of the environmental constraints chosen to predict the composition of local macroinvertebrate assemblages. Frequently, only a limited number of environmental variables is selected, which is problematic, as it is inherently assumed that none of the ignored environmental variables would affect macroinvertebrate assemblage composition. This is further complicated by the fact that environmental stressors are highly interrelated and may act directly as well as indirectly on macroinvertebrate assemblages. Therefore a key challenge is to disentangle the contribution of individual and combined stressors to the overall effect on macroinvertebrate assemblages (Calapez et al., 2017;Elbrecht et al., 2016;Smeti et al., 2019;Tockner et al., 2010;Villeneuve et al., 2018).
To come to the required improvements, for a number of environmental variables, ecologically meaningful expression formats are described below. Environmental variables should be chosen to represent multiple spatial and temporal scales and multiple environmental categories, such as hydrology, chemistry, morphology and system conditions. In addition, ecologically relevant expression formats should represent extreme events and variability in the environment, as opposed to average values only, e.g. it is favoured to explicitly specify stressors such as the number and size of peak flows, rather than a general measure such as discharge. These prerequisites result in a broad set of environmental stressors that could be considered when simulating macroinvertebrate assemblage responses on multiple spatial and temporal scales (Fig. 1). However, to avoid the risk of overfitting due to the high number of stressors, the selection has to be limited to those variables for which an effect on the assemblage can be reasonably expected based on available ecological knowledge (Burnham and Anderson, 1998).
Ecologically relevant expression formats have been described for macroinvertebrates in general Poff 1997), as well as for specific communities (Verberk et al., 2012;. Examples of expression formats are also given by Bovee et al. (1998), Kail et al. (2015) and Weijers et al. (2012), who described models that include expression formats more specifically addressing spatial and temporal variability and frequencies of extreme events.
Extreme events and high variations in environmental variables may have large impacts, as macroinvertebrates prefer a constant environment, especially considering oxygen concentration and flow conditions (Verdonschot and van den Hoorn, 2010). Within the category hydrology, the magnitude, frequency, timing and duration of extreme events were listed as ecologically relevant expression formats Poff, 1997), instead of using mean velocities or discharges only. For oxygen concentration, in current simulation models mostly mean values are used, but these do not represent diurnal oxygen dynamics (van der Lee et al., 2018). Occurrence of anoxic events can be lethal and environmental variables that represent the frequency and duration of such events are therefore more expressive than the mean oxygen concentration (Grieshaber et al., 1994). Although technically demanding, continuous measurements of oxygen concentrations can provide such information.
Drivers of local nutrient concentrations act on a regional scale. For example, there is a strong relationship between nutrient concentrations in surface waters and the surrounding land use (Boyer et al., 2002;De Wit, 1999). Moreover, extremely high nutrient concentrations may not directly pose stress on macroinvertebrates, but indirectly via increased autotrophic biomass production, lower oxygen concentrations and changed food availability, all impacting stream integrity (Allan, 2004;Allan et al., 1997;Delong and Brusven, 1998). Therefore, as opposed to using highly variable local nutrient concentrations, regional characteristics such as land use may be used to integrate nutrient-related stress over space and time.
Substrate patchiness and variation in current velocity at the stream bed are also examples of meaningful, but rarely pictured variables (Townsend, 1989). Hence, collection of data on substrates at reach scales would contribute to a meaningful dataset to be used for prediction. Furthermore, the impact of stressors acting on macroinvertebrates are life-stage specific (Kefford et al., 2007;Stuijfzand et al., 2000;van der Lee et al., 2020). Therefore, the life cycle of aquatic organisms should be considered to include stressors specific to sensitive aquatic and terrestrial life stages.
The relevant expression formats of environmental constraints described above answer to our requirements, and can be included in a conceptual model as suitable building blocks.

Stressor interactions
Stressor interactions can be of major importance in aquatic ecosystems , because they can markedly change the expected responses to combinations of individual stressors. Hence, including these interactions in model simulations increases the ecological realism of the outcomes.
Synergistic and antagonistic effects may occur when stressors act simultaneously (Coors and de Meester, 2008;Jackson et al., 2016). Examples include the synergistic effects of low flow and oxygen depletion on macroinvertebrate communities, causing increased drift rates (Calapez et al., 2017), and the increased negative impact of sediment input at reduced flow on benthic invertebrates (Matthaei et al., 2010). These multiple-stressor interactions have not yet been accounted for in macroinvertebrate-simulation models, instead, stressors are treated one by one and their effects are subsequently combined additively. Therefore, including stressor interactions increases ecological realism in simulation outcomes.
Interactions are specific per stressor combination, water body type, assemblage, and region. A number of these interactions has been described (Calapez et al., 2017;Davis et al., 2018;Elbrecht et al., 2016;Matthaei et al., 2010), but understanding of the mechanisms behind these interactions is still limited (Friberg, 2010;Ormerod et al., 2010). Therefore, it is not possible yet to fully represent stressor interactions in models for simulating macroinvertebrate responses to multiple stressors, even though an ideal future model would answer to this requirement.
Yet, there are some very simplified approaches available that can be used to partly meet the requirements. In previous model applications it was attempted to incorporate interactions between multiple stressors by combining univariate response curves (Bovee et al., 1998;Gobeyn et al., 2017;Tomsic et al., 2007), which describe the relationship between an environmental gradient (e.g. oxygen concentration) and the response of an organism (e.g. a measure of species fitness). In these previous models, responses to single stressors were combined into an aggregated response to multiple stressors. One aggregation method is taking the geometric mean instead of the arithmetic mean of response values (Gobeyn et al., 2017;Tomsic et al., 2007) to avoid compensation of non-suitable values (high amount of stress) by suitable values of other stressors. No compensation of non-suitable values is also possible when applying the 'one out, all out' principle or minimum aggregation, taking the stressor value with the poorest condition as an aggregated response value (Bennetsen et al., 2016;Bovee et al., 1998;Langhans et al., 2014;Morales et al., 2006). Even though these methods do not yet account for the complex interactions that may occur in reality, they represent first steps to account for stressor interactions and can therefore be used as provisional building blocks.

Uncertainty
Uncertainty is arising from variable sources of data, being empirical, correlational, or originating from expert knowledge. In addition, stochasticity in natural processes increases variation in the data. Accounting for and communicating these uncertainties, when they are unavoidable, is vital for a correct interpretation of the model outcomes.
To deal with this uncertainty, many good modelling practices have been described (Refsgaard & Henriksen 2004;Rose et al., 2015). It may seem obvious that such practices are followed in the design and application of models, but in practice, this is not always the case (Refsgaard et al., 2005). Therefore, they are mentioned again here, as they provide necessary means to deal with uncertainties in data. Good modelling practices involve a proper documentation of the uncertainties related with the input data, transparency about the assumptions underlying the model design and a sensitivity analysis for evaluating the robustness of the modelling results and for identifying how changes in input parameters influence the behaviour of a model. Sensitivity and uncertainty analyses can be very useful for refining the model and correctly interpreting model outcomes (Rose et al., 2015). These good modelling practices are therefore important steps in the development and application of a model.
In addition, an appropriate model design should be chosen, suitable to deal with these uncertainties. For instance, Bayesian Belief Networks (BBNs) are graphical network structures which can be used to relate environmental variables statistically to predict organism responses using probabilistic relationships. An advantage of this model design is that it explicitly deals with the uncertainty of input variables and small or incomplete datasets by calculating the probability distribution of the outcome (Uusitalo 2007). These models seem also promising because they can include expert-based and other knowledge, which would otherwise be lost, for defining the model structure and for formulating the probabilistic relationships which define the dependencies within the model (McCann et al., 2006). Such model structures have been previously applied to macroinvertebrates (Adriaenssens et al., 2004;Boets et al., 2015;Forio et al., 2015;McLaughlin and Reckhow, 2017;Skeffington et al., 2014) and several other biotic endpoints representing ecological water quality (Johns et al., 2017;Marcot et al., 2001). Thus, this model design is a useful building block for a conceptual model for simulating macroinvertebrate responses to multiple stressors.
An ideal multiscale, multifactorial model for simulating the responses of regional and water body type-specific freshwater macroinvertebrate assemblages to multiple stressors should take all these requirements into account. Yet, there is currently no model available that meets these requirements. Hence, the next step is to select the appropriate building blocks from available methods to come to such an ideal model.

A conceptual model for simulating responses of freshwater macroinvertebrate assemblages to multiple stressors
From the selected building blocks listed in the previous section, we assembled a conceptual model. This conceptual model is based on the filterc oncept, which is a structural approach with hierarchical levels of processes driving the species assemblages at regional, catchment, reach and habitat scales (Lake et al., 2007;Poff, 1997). The filter concept states that the realized local assemblage is formed from the regional species pool by species that overcome dispersal constraints (by capacity for dispersal and spatial connectivity), environmental constraints (filters at hierarchical landscape scales which fit the species preferences) and biotic constraints (such as predation, competition and facilitation). In turn, this theory is derived from the habitat templet concept explaining how characteristic species traits are related to to heterogeneous environments (Southwood, 1977). Consequently, in the conceptual model, stressors affecting macroinvertebrate assemblages ranging from catchment down to local habitat can be considered.
The proposed conceptual method consists of several structural elements, referred to with numbers in Fig. 2. The water body type and regional species pool determine the design of the network and the relationships between nodes (1). The model core consists of three domains, representing the environmental, dispersal and biotic constraints, representing filters shaping the local macroinvertebrate assemblage (Lake et al., 2007). The main domain is an extended habitat-model within a BBN-framework (2). Relationships between the nodes of the network are informed by using existing mechanistic and statistical models. The environmental constraints domain is run for multiple species of interest within the target species pool. Element (3) represents dispersal constraints. Together with the environmental constraints, these elements determine which species from the regional species pool end up in the local species pool. Element (4) represents biotic constraints. The output of this domain feeds into the last step (5) where an impact assessment for the macroinvertebrate assemblage is made. These structural elements are described in more detail below.
The boundary conditions (1) consist of a specific water body type and geographical region. For these preconditions, a target macroinvertebrate assemblage is determined. This target assemblage can be taken from existing regional typologies used for assessing the monitored state of a water body. Furthermore, the water body type and the regional species pool are defining the design of the network, the choice of environmental constraints and the relationships between the nodes.
The main domain describes the environmental constraints (2) shaping the macroinvertebrate assemblage. This part of the model is designed for a single species, as single-species assessment allows using specific data on habitat requirements and other species-specific responses (Ferrier and Guisan, 2006). First, predictors are selected based on a set of environmental stressors acting directly on the species (Feld and Hering, 2007;Frissell et al., 1986;Verdonschot, 2000), taking into account the spatial and temporal scales and variability on which these stressors are acting, as described in detail above. Then, environmental drivers placed hierarchically above these stressors are included, allowing for testing and comparing effects of restoration measures on different spatial scales . The indirect and direct predictors of macroinvertebrate responses make up a network of nodes, which can be connected using probabilistic relationships, forming a Bayesian Belief Network (Landuyt et al., 2013;McCann et al., 2006). Effects of interactions between direct stressors on macroinvertebrate responses can be accounted for in the probabilistic relationships within the design of Conditional Probability Tables, which describe how the values of multiple nodes combine into an output value of the underlying node. This network results in a habitat suitability index for specific species as a measure to express biological responses to multiple environmental constraints.
To move from predictions for single species to predictions for a species assemblage, this domain can be run for multiple species of interest within the target assemblage. The habitat suitability scores for multiple species can be combined in several ways, for example by taking the minimum habitat suitability score, which would act as a bottleneck on the target assemblage (Milhous and Waddle, 2001). The output of the domain on environmental constraints is thus a combined measure of the habitat suitability for a target assemblage, which provides a first input to determine the local species pool.
The third domain describes the dispersal constraints (3) acting on a regional species pool to form the local species assemblage. The input of this domain is the regional species pool, which can be deduced from distribution maps, monitoring data and regional typological data, modified by environmental constraints that influence the spatial connectivity. Whether the species within the target species pool are actually present at the local site of interest, depends on the spatial connectivity in the landscape and on the colonisation capacity of the organisms (Sarremejane et al., 2017). As described above, there is a number of simplified methods to account for the effects of dispersal constraints that can be included in this conceptual model. Colonisation capacity is expressed by certain species traits, allowing the organism to travel a certain distance in the water, actively or passively by drift, or via air (Townsend and Hildrew, 1994). This can be simplified by including a general geographical dispersal limitation term (Bennetsen et al., 2016). Spatial connectivity would ideally be modelled in a spatially explicit way to account for the stream network, including the locations of source populations and migration barriers Kuemmerlen et al., 2019), but can be summarized here by taking the distance which can be travelled via water without obstructions like dams or weirs (van Puijenbroek et al., 2019). Using these simplified measures for dispersal constraints, the locally available species pool can be deduced from the regional species pool. Together with the output from the environmental constraints, this serves as a second input to determine the local species pool.
The previous steps deduce the local species pool from the regional species pool. Next, the effects of biotic constraints (4) can be considered. Constructing an entire model to represent this filter is not aimed for because of the added complexity and thus decreased usability. Instead, trophic and non-trophic interactions can be represented by using simple measures. In case of known trophic interactions with specific other species, the presence of this other species can be used as another predictor representing competition for food sources or predation (Hirzel and Le Lay, 2008). Non-trophic interactions can be included by estimating biotic resistance to invasions . These terms yield a simplified measure of the stress originating from biotic constraints, to be used in the final impact assessment. The three domains yield different output formats. The output of the first domain, quantifying environmental constraints, takes the form of a combined habitat suitability score for a target macroinvertebrate assemblage. The domain on dispersal constraints indicates if the species of interest can be present in the local species pool, being a measure for the colonization potential. The environmental and dispersal constraints together determine the local species pool. From the third domain follows a measure for stress from biotic interactions. This can be combined with the output from the previous domains in a probabilistic network to make a final assessment of the impact of all joint stressors (5), originating from environmental, dispersal and biotic constraints on the macroinvertebrate species assemblage. To combine the multiple measures of stress we can learn from aggregation methods of assessment scores, focusing on the minimum, mean or a mixture of the single scores (Langhans et al., 2014). Because of the BBN-method adopted, uncertainty of the output can be shown in the associated probability curve. The proposed conceptual model is a flexible framework that is meant to be specified per region and water body type. The next steps may include showing its applicability with test cases. Yet, some technical challenges and knowledge gaps remain which should be addressed before an actual model that answers all requirements posed by the conceptual model can be built. One of these remaining challenges is the lack of species-specific stressor data. Especially knowledge on the interactions between stressors for specific species in specific surroundings is limited. Still, data availability is assumed to increase in view of regulations that request widespread monitoring for ecological water quality assessment . By tackling the described challenges, future studies will be part of the development towards the proposed model structure for simulating macroinvertebrate responses to multiple stressors. The application of the conceptual model in real world case studies may help water management by providing an increased understanding of the system, aiding to optimise the use of resources and the selection of restoration practices. Furthermore, actual models based on the conceptual model can be used as a diagnostic tool to detect the main stressors affecting macroinvertebrate assemblages and to analyse different scenarios to predict the response of organisms to stressors and intended restoration measures.

Conclusions
Currently, no ideal model is available for simulating responses of freshwater macroinvertebrate assemblages to multiple stressors. Hence, essential missing steps should be included in the development of future models. Here, we demonstrated that these essential building blocks are available and we proposed a conceptual model structure for simulating macroinvertebrate responses to multiple stressors. Yet, considerable efforts should still be devoted to better include interactions between stressors. This will bring simulating responses of freshwater macroinvertebrate assemblages to multiple stressors further, which in turn will contribute to the focus of management resources and restoration practices.

Author contributions
JV wrote the main manuscript text. All authors were involved in the study design, contributed to the writing of the manuscript and gave final approval for publication.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.