{'fraction_non_alphanumeric': 0.04335714838410152, 'fraction_numerical': 0.023110056485517526, 'mean_word_length': 4.216198035786358, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 21, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 21, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The effects of invasive epibionts on crab-mussel communities: a theoretical approach to understand mussel population decline

Blue mussels (Mytilus edulis) are an important keystone species that have been declining in the Gulf of Maine. This could be attributed to a variety of complex factors such as indirect effects due to invasion by epibionts, which remains unexplored mathematically. Based on classical optimal foraging theory and anti-fouling defense mechanisms of mussels, we derive an ODE model for crab-mussel interactions in the presence of an invasive epibiont, Didemnum vexillum. The dynamical analysis leads to results on stability, global boundedness and bifurcations of the model. Next, via optimal control methods we predict various ecological outcomes. Our results have key implications for preserving mussel populations in the advent of invasion by non-native epibionts. In particular they help us understand the changing dynamics of local predator-prey communities, due to indirect effects that epibionts confer.

1. Introduction 1.1. Background. Blue mussels (Mytilus edulis) are an ecologically and economically important species [1,2,3]. They play several roles in marine ecosystems: as important prey for many species, such as crabs, shorebirds, sea stars, and gastropod molluscs [3,4,5]; as nutrient recyclers and pollution indicators [6]; and as a keystone species, serving as habitat for benthic infaunal organisms [3,7]. However, M. edulis has declined in the Gulf of Maine by over 60% since the 1970s [8]. Mussel post-larval settlement, consistent with this observation, has also declined [11]. The reasons for this decline are unclear, but are almost certainly complex. Thus a clearer understanding of the ecological factors that influence mussel populations is needed.
A primary cause of a species population decline is predation. Invasive predators, like the green crab (Carcinus maenas) and the Asian shore crab (Hemigrapsus sanguineus), readily prey on the blue mussel [12,13,15]. However, mussel size limits crab predation, with crabs consuming mussel prey below 70 mm in shell length [16]. Furthermore, mussels have also adapted to crab predation by thickening their shells in response to novel predator presence, in extremely short time periods [15]. Furthermore, substrate complexity reduces predation on mussels as increasingly complex habitats provide refuge from crab predation [17]. Thus while predation has put considerable pressure on mussel populations, rapidly evolving defense mechanisms, escape from predation via growth, and physical refuges have counteracted predator impacts.
Though mussels do have the aforementioned protections against predation, they are still in decline. Curiously, in the 1970s, an introduced ascidian species Didemnum vexillum, arrived in the Gulf of Maine [18]. D. vexillum is a colonial ascidian that is dominant as a competitor for substratum, prolifically laying down mat-like structures on any hard substrate [22]. Consequently, it acts as an epibiont (fouling organism) on M. edulis [11,12]. Epibionts impact predator-prey communities indirectly by affecting predation rates on basibionts. D. vexillum in particular has chemical anti-predatory defenses. If a crab predator attempts to break off pieces of the D. vexillum colony to reach the mussel, D. vexillum releases secondary metabolites and sulfuric acid that deters the crab [22]. This mechanism by which the epibiont protects the mussel from crab predation is known as associational resistance [24]. While it appears to protect mussels from crab predation, D. vexillum also negatively affects mussel fecundity and fitness, resulting in fewer progeny [11]. Thus D. vexillum has both a positive and negative impact on mussel populations.
Given this complex relationship, we ask, • Could the introduced epibiont D. vexillum change the predator-prey dynamics in an established local crab-mussel community? • Could the net effect of positive and negative impacts from D. vexillum epibiosis, lead to mussel population decline? Our analysis is (to the best of our knowledge) the first mathematical investigation of predator-prey dynamics under pressure of fouling from epibionts in a crab-mussel community. Herein, • we derive a predator-prey model for crab-mussel interactions, given that clearly a certain size of mussel is preferred or "optimal" for the crab, • we consider the effects of an invasive epibiont by meshing OFT and antfouling defense of the mussels, • we model the effects of associational resistance and reduced fecundity, due to the epibiont, and • we next investigate dynamical aspects of the model, and use optimal control theory to predict various outcomes. In adult mussels the protective periostracum, which inhibits epibiont settlement when present [26] tends to wear off due to age, decay and abrasion. Consequently, the periostracum is more prevalent on newer regions of the shell, while absent on older regions [27]. This means juvenile mussels are less likely to be overgrown with epibionts than are adult mussels. When crabs forage for mussels, they typically prefer a medium sized adult. But this preferred size tends to get easily overgrown by epibionts. Epibionts can alter the prey size choice of predators, including crabs, in experiments [22,28,29,30], though this has not yet been tested with D. vexillum. We hypothesize in the current manuscript that D. vexillum can change the feeding preference of crabs away from mid-size adult mussels (that are easily overgrown and therefore less likely to be eaten) towards juvenile mussels (which are less likely to be overgrown and so are easier targets), even though the latter are not the crab's preferred food source. To elucidate our approach we survey some classical results from OFT.
1.2. Optimal Foraging Theory. Optimal foraging theory (OFT) predicts how animals behave when they forage for food. It is well known that predators optimize feeding strategies to maximize energy intake [31]. Essentially, predators try to gain the most energy from their prey by expending the least amount of energy in the hunting process. For a crab foraging for mussels this amounts to maximizing (1) e h = energy gained from mussel intake handling time of mussel .
This translates to a medium-sized adult mussel as the optimally preferred prey by adult crabs. While large mussels have a potentially high source of energy, they take a much longer time to open and consume than smaller mussels. Small mussels, conversely, take a short time to consume, but they offer very little reward. Even so, juvenile mussels are readily consumed by many species, including green crabs and dogwhelk [12]. We refer the reader to the mathematical treatment by Krivan [32,33], who describes, via a three species ODE model, in which a predator hunts (optimally) for two prey species u and v, where u is favored to v. The term u 1 denotes the attack rate with which prey u is hunted, and u 2 denotes the attack rate with which prey v is hunted. The analysis presented in [33] draws from classic results in OFT and shows that, in order to maximize e h , the optimal pair of (u 1 , u 2 ) is given by u 1 = 1, u 2 = 0 if u > u * or u 1 = 1, u 2 = 1 if u < u * , or u 1 = 1, 0 < u 2 < 1 if u = u * , where u * is the critical density for switching.
1.2.1. OFT in the presence of epibionts. Predators are known to switch prey if preferred prey drop below a threshold density [35]. For example, fish species have been shown to switch habitats if foraging in one habitat becomes less fruitful than in another habitat [36]. Theoretical studies also support this [32]. In our context, If e1 h1 > e2 h2 , the adult mussels are preferred to juveniles as the optimal prey for crabs. If e1 h1 < e2 h2 , juveniles mussels are optimal to hunt. Conjecture 1. An increase in epibiont density will cause crabs to switch from adult mussels to juveniles even though the juvenile is less preferred and the adult mussel density remains high.
To this end we first split the mussel class into adults and juveniles (we assume the juveniles have the protective periostracum whereas the adults do not). A crab species preying on two separate classes of adult and juveniles mussels (with adults being the preferred food type), places this in a classic one predator-two prey setting [32]. Our hypotheses for OFT are as follows: (1) In the absence of epibionts, crab-mussel interactions follow classical OFT.
That is, adult mussels will be attacked with rate 1, whilst the less preferred juveniles will not be attacked (u 1 = 1, u 2 = 0). We claim this is the only optimal strategy for the crab, as long as d 1 > e2 h2 , where d 1 is the death rate of prey type 1.
(2) There is a change to the classical case, under pressure of epibiosis from D.
vexillum. (3) If a crab were presented with a preferred adult mussel overgrown with D.
vexillum, it would switch to a prey of a less optimal size, even if the overall adult mussel density was high (assuming their was uniform overgrowing of all adult mussels). (4) The switch would be to juvenile mussels, which we know are almost never overgrown because of their intact protective peristrocum. That is in the (e = K) case, we have (u 1 = 0, u 2 = 1). We claim this is the only optimal strategy for the crab, as long as d 1 > e1 h1 . (5) This will in turn directly affect the feedback loop to the adult mussel population, given that juveniles are transitioning to adults.
The above is rigorously shown in appendix section 8.1. mussels size crab preference best size mussels size crab preference Figure 1. Classic crab preference (without epibiont) shown in left panel. Crabs were presented with clean and overgrown mussels [11,12], and the handling times of the overgrown ones were up to 6 times greater than the clean ones (seen in center panel). These experimental results [11,12] lead us to conjecture that crabs will switch to the smaller uncovered juvenile mussels under pressure of epibiosis, shown in right panel.

Mathematical Formulation
Our goal is to derive a mathematical model that best captures our hypotheses. To this end, we make the following assumptions, (1) An epibiont has invaded into a local crab-mussel community, and is growing logistically. It will eventually reach a critical carrying capacity.
(2) We model the pressure from epibiosis, in terms of the attack rates u 1 , u 2 . That is we assume these are dependent on the epibiont density. As a simple first approach we assume Thus, without any epibiont presence (e = 0), the adult mussel is the only prey eaten and the juvenile is not eaten at all because there is not enough energy gain for the effort involved, so u 1 = 1, u 2 = 0, in line with classical theory [33]. However, this starts to change as the epibiont starts to overgrow the mussels. When the epibiont is at carrying capacity, e = K, we assume the adult mussels are completely overgrown, and thus is not consumed at all. The crab switches completely to juveniles, so that u 1 = 0, u 2 = 1.
(3) We model the decreasing fecundity in mussels due to epibiont cover by considering a growth rate a = a(e). We consider a = a(e) = a Hence as epibionts get to carrying capacity e = K, this growth rate is cut in half and becomes a/2. (4) We assume the intraspecies competition is present only in adult mussels, and not juveniles [49]. (5) We assume the search rates λ 1 , λ 2 to be the same, and normalized to 1, so λ 1 = λ 2 = 1. Based on the above assumptions, we have the following system of differential equations, where Here C, M A , M J are the densities of crabs, adult mussels and juvenile mussels population at a given time t respectively. The population density of D. vexillum is e, while d 1 is the mortality rate of the crab, e 1 , e 2 is the energy gain to the crab from preying on the adult mussel and juvenile mussel respectively, h 1 , h 2 are the handling time of the adult mussel and juvenile mussel respectively, b is the rate at which juveniles leave the juvenile class and become adults, a is the rate at which juveniles are produced, δ 1 measures the intraspecific competition among adult mussels, b 1 is the intrinsic rate of growth of the epibiont population, and K is its carrying capacity.

Dynamical Analysis
3.1. Boundedness. The equation for e is bounded trivially by K. Addition of (129)-(131), and given the fact that e 1 < 1 and e 2 < 1, yields: Thus, comparison with a logistic ODE yields: and thus we can state the following theorem: Theorem 3.1. Consider the crab-mussel system (129)-(132). The solutions (C, M A , M J , e), satisfy the following uniform bounds ). 3.2. Equilibrium and Local Stability with no epibiont. We now consider the existence and stability of the equilibrium for the system when there is no epibiont(e = 0). The system is simplified as Two equilibria, (0, 0, 0) and (0, a δ1 , a 2 δ1b ), on the boundary and one interior equilibrium (C * , M * A , M * J ). It is easy to see (0, 0, 0) is unstable.
The interior equilibrium is given by .
We next state the following theorem Theorem 3.2. Consider the crab-mussel system (129)-(132), without the presence of an epibiont, that is when e = 0. There exists an interior steady state (C * , M * A , M * J ), which is locally asymptotically stable under the following criteria, The proof is relegated to the appendix section 8.2.
Remark 1. Note, the second condition d 1 > e2 h2 is not a result of the Routh-Hurwitz criterion, rather it follows from lemma 8.1 in appendix section 8.1. We enforce it so that the attack rates should be as predicted via classical OFT.
Consider the interior equilibrium, i.e. (C * , M * A , M * J , e * ). It is easy to see e * = K in the equilibrium state.Then we have u 1 (e * ) = 0, u 2 (e * ) = 1, a(e * ) = a 2 . To get (C * , M * A , M * J , e * ) explicitly, it is equivalent to solve the following equations: Thus the interior equilibrium is given by Therefore, the feasibility criteria for this system is We next state the following theorem Theorem 3.3. Consider the crab-mussel system (129)-(132), when the epibiont has reached equilibrium, that is e = K. There exists an interior steady state (C * , M * A , M * J , K), which is locally asymptotically stable under the following criteria, The proof relies on the Routh-Hurwitz criterion [37], and is relegated to the appendix section 8.3.

Remark 2.
Note that the condition d 1 > e1 h1 is not a result of the Routh-Hurwitz criterion, rather it follows from lemma 8.1 in appendix section 8.1. In a sense we enforce it so that the attack rates should be as predicted via classical OFT.
3.4. Global stability. We now derive some results on the global stability of the internal equilibrium.
Theorem 3.4. Consider the model (10)- (12). There exists an > 0, s.t. the internal equilbrium point, (C * , M * A , M * J ), is globally asymptotically stable under the following parametric restriction Proof. Consider M A = M A + , Under this transformation we have the following transformed system It is enough to show the new system (29)-(31) is globally asymptotically stable. The equilibrium of (29)-(31) is given by Note that solutions to this new system are feasible if We define the following Lyapunov function, We Since M A is bounded by a δ1 − from the feasibility conditions, we show Therefore we only need to show and this requires we have 1 2 < e 1 < 1.Then the system (10)- (12) is globally stable if

Hopf Bifurcation
Now we will investigate the Hopf bifurcation for the system in terms of parameter a. In this paper, we will follow the method developed by Liu [38]. Firstly, let us consider the system (129)-(131), without the presence of an epibiont (e = 0), that is when e = 0. The Hopf bifurcation at a = a * can occur if A 2 (a * ), A 0 (a * ), and φ(a * ) = A 2 (a * )A 1 (a * ) − A 0 (a * ) are smooth functions of a in an open interval of a * ∈ R such that: da | a=a * = 0. We check the above, in appendix section 8.4, to state the following theorem

Optimal Control
In this section our goal is to investigate mechanisms in our crab-mussel-epibiont system, that, if controlled, could lead to optimal levels of crab or mussel densities. We assume that the attack rates u 1 , u 2 are not known a priori and enter the system as time-dependent controls. They no longer depend on the epibiont density. Instead we assume that the handling time depends on the epibiont density e in the following way, where These responses are for the range 0 ≤ e ≤ K. Increase in epibiont density still negatively effects mussel fecundity and the handling time for adult mussels increases with increasing epibiont density. This has a twofold advantage. We can visualise the system from the crab's point of view. That is, the crab can "optimally" control its attack rate, to reach the best possible population density. Also we can visualise the system from the mussel's point of view. That is, the mussel can induce defenses or other mechanisms, that would alter the attack rate of the crab, thus enabling the mussel population density to reach optimum levels. Our model takes the following form, We next derive optimal strategies for three objective functions, where we maximize both crab and mussel populations. To simplify the calculation, we will consider the case when e = K, which is when the epibiont achieves carrying capacity. In this case our system reduces to

5.1.
Maximizing crab denisty w.r.t. attacking rates. To maximize density of the crab, the density of juvenile mussels (thus leading to more adult mussels, its favored food) should be maximized. Crab attack rates should be miminized on the juvenile mussels, as they are less preferred by the crab. Thus we choose the following objective functional, . and we search for the optimal controls in the set U where The goal is to seek an optimal (u * 1 , u * 2 ) s.t., We can state the following existence theorem, Theorem 5.1. Consider the optimal control problem (47)- (49).
Proof. The compactness (closed and bounded in ODE case) of the functional J 1 follows from the global boundedness of the state variables via theorem 3.3, and the boundedness assumption on the controls. Also the functional J 1 is concave in the argument u 2 . This is easily verified via standard application [39]. These in conjunction give the existence of an optimal control via application of classical one predator-two prey theory [40].
In order to derive necessary conditions on the optimal control, we use Pontryagin's maximum principle (PMP). The Hamiltonian for our problem is given by We use the Hamiltonian to find a differential equation of the adjoint λ i , i = 1, 2, 3.
with the transversality condition given as Considering the optimality conditions, the Hamiltonian function is differentiated with respect to control variables u 1 and u 2 resulting in We find a characterization of u * 1 by considering three cases: When the control is at the upper bound,u 11 is strictly greater than 1. When the control is at the lower bound, the solution of u 11 is strictly less than 0. Similarly for u * 2 . Thus a compact way of writing the optimal control is where u 11 and u 21 are given by with We can thus state the following theorem, Theorem 5.2. An optimal control (u * 1 , u * 2 ) ∈ U for the system (47)-(49) that maximises the objective functional J 1 is characterised by (59).

5.2.
Maximizing mussel density w.r.t. attacking rates. To maximinize mussel density, the attack rate on the adult mussels should be minimized. We choose the following objective function, We can state the following existence theorem, Theorem 5.3. Consider the optimal control problem (47)- (49). There exists (u * 1 , u * 2 ) ∈ U s.t.
For the details of the proof of the above necessary conditions and forms of u 12 , u 22 the reader is refered to the appendix section 8.5.

5.3.
Maximizing mussel density w.r.t. intraspecific competition rate. In this approach we view the competition coefficient as a control. To reach certain optimal population densities, the mussels would maximise the densities of both adult and juvenile groups whilst minimising intraspecific competition. To this end our system reduces to We choose the following objective function, . and we search for the optimal controls in the set U 1 . Where, We can state the following existence theorem, Theorem 5.5. Consider the optimal control problem (66)-(68). There exist (u * 1 , u * 2 ) ∈ U s.t. (71) The proof is similar to theorem 5.1. We can next state Theorem 5.6. An optimal control (u * 1 , u * 2 ) ∈ U for the system (66)-(68) that maximises the objective function J 3 is characterised by . For the details of the proof of the above necessary conditions and forms of u 12 , u 22 the reader is refered to the appendix section 8.6. We set h 1 = 2 since h 1 = 1 + e K , however, if we just assume h 1 as a constant and increase h 1 and keep other parameters the same, we found the optimal control u 2 always to be 0, and u 1 decreases and gradually become stable. In fact, when h 1 achieves to some critical value, u 1 begins to increase slightly. Based on e = K, as for system (129)-(132) , the optimal foraging strategies are u 1 = 0, u 2 = 1. However, for the control system (43)-(46) , J 1 (u 1 , u 2 ) will be maximized when u 1 = 0.4343, u 2 = 0; J 2 (u 1 , u 2 ) will be maximized when u 1 = 0, u 2 = 0.9121 and J 3 (δ 1 ) will be maximized when δ 1 = 0 with data set we mentioned.

Discussion and Conclusion
Epibiotic invasive species often have anti-predator defenses that are behavioral, chemical, or mechanical [41,42], giving them a survival advantage in a novel habitat because potential predators avoid using them as a food source [43,44]. While this provides a benefit to the basibiont, it impacts other members of the community, including predators of the basibiont as they may show lower preference for basibionts that are overgrown by invasive epibionts [12,45]. However, the effects of epibionts on basibionts are not always positive. Many times the epibiont may attract predators resulting in consumption of the epibiont, which automatically leads to consumption of the basibiont. This is refered to in the literature as "shared doom" [24]. Epibionts can also negatively affect basibiont fecundity and fitness, resulting in fewer offspring [11]. In essence, invasion of predator-prey communities by epibionts is complex, and warrants a thorough mathematical investigation of their impact on predator-prey interactions and populations.
Population cycles are common in predator-prey communities, and although these are possible in our model without epibionts, extensive numerical simulations indicate that at carrying capacity e = K, a Hopf bifurcation is not possible. This points to the epibiont having a stabilizing influence in that it can eliminate population oscillations. A rigorous proof of this is an interesting future direction. Within our study, theorem 3.4 tells us that if the energy gain from the adult mussel is in a certain critical region 0 < e 1 < 1, then one has global stability; even very large perturbations would still allow the system to return to its base state.
Our central question focuses on the effect of the introduced epibiont on the population densities of the local crab-mussel communities. Could high epibiont density lead to lower mussel populations (and so subsequently lower crab populations)? To answer this we compare the equilibrium levels of the juvenile mussel population, "no epibiont" case versus "epibiont reaches carrying capacity" case. If the epibionts do have an adverse effect then we would have M * J | (e=K) < M * J | (e=0) . Comparing these yields, Although we know e 1 > e 2 , under high epibiont density (e = K), we have h 1 >> h 2 , thus even if a < b, (74) could easily hold meaning that there is an adverse effect on the juvenile mussel density via epibiont presence, leading to fewer adults subsequently, and so epibionts could clearly be a factor in mussel population declines as seen via data from the Gulf of Maine [8]. Such decline could eventually lead to crab population decline as well, if the crab species is a specialist on mussels. However, the effects of epibionts on mussel fecundity could also be a cause of predator decline. In order to understand the effects of epibiosis, we investigate the equilibrium density of the crab populations for the "no epibiont" case versus the "epibiont reaches carrying capacity" case. What we note is Clearly, as epibiont cover reduces mussel fecundity from a, to a/2, this directly affects the crab population. In the (e = K) case there is an increase by a factor of only a 2 , as opposed to a factor of a in the (e = 0) case. Thus reduced fecundity in mussels due to epibiont cover, can also reduce crab populations as well.
We assume logistic growth in the epibiont density. Although in the Gulf of Maine epibiont density fluctuates seasonally, our model could be a useful predictive tool in periods where logistic growth is seen. In locations such as Japan and New Zealand, the epibiont D. vexillum grows logistically [11] due to water temperatures staying above the threshold for D. vexillum viability.
We also use optimal control theory to visualise various optimal scenarios to maximize each crab and mussel densities. Herein, we change the problem slightly, and assume the attack rates u 1 , u 2 are not known a priori, but are time-dependent. Our objective is to explore various scenarios that a species of crab or mussel might attempt to optimise, by manipulating the attack rates. Epibionts are assumed to be present, and their effect is modeled via increasing the handling time h 1 of adult mussels, as epibiont density increase.
Simulations suggest (Fig. 4a) that even under high epibiont density (which in this scenario amounts to doubled h 1 ), the crab should not attack juvenile mussels, but attempt to attack adults. Fig. 4b demonstrates, that what is optimal for the mussel is if u 1 = 0, and so the adult mussel must induce defenses to reduce u 1 , even if it realistically cannot drive that rate to zero. This confirms experimental results of rapid shell thickening by mussels, seen via [15]. Fig. 5a looks at the attack rate on adut mussels u 1 , as h 1 changes. Here, we are trying to maximize mussel populations, and u 1 decreases as h 1 increases, as expected. However, u 1 is approximately 0.17; that is, it does not change significantly if handling times become very large. Curiously, it goes up ever so slightly as shown in Fig. 5b. This likely corresponds to the adult mussel thickening its shell just enough to increase handling time by the crab. Fig. 8 looks at the attack rate on juvenile mussels u 2 , as h 1 changes. Here again we are trying to maximize mussel populations. When handling time on adults is low, juvenile mussels are protected from crab predation due to the predators preference for larger mussels. However, when handling time on adults is high (greater than 2.1 in this simulation), it is likely that crabs would switch to the juvenile; therefore, juveniles must be able to disperse or seek refuge in order to bring attack rates on the juvenile mussels to zero (in turn maximizing their population size). Young mussels drift in the water column until they reach a size of approximately 2.5mm, then they settle on a filamentous algal substrate [3]. Some mussel species settle on algal substrate until they are 30 mm in length [48]. This substrate acts as refuge and must be available for juveniles in order to maximize mussel populations. However, with degradation of suitable habitat, the opportunities for escape from crab populations becomes diminished. Major disturbance events, either natural or anthropogenic, in conjunction with invasion by substrate-smothering colonial species and voracious predators, are likely to decrease opportunities for escape. Endeavors to model predator-prey systems incorporating prey refuge may yield surprising results on stability [23], [25]. Thus it would be very interesting to model refuge effects for the juvenile mussels herein.
As a future direction in modelling the crab-mussel-epibiont interaction, we would also like to examine interference effects [19,21,20]. This effect has often seen to be stabilising [34], and thus modeling interference among the crab population, at high epibiont density is also realistic. Note, Theorem 5.6 suggests that eliminating intraspecific competition among mussels is optimal from their point of view, and yields a maximum density, if there was no competition present. Future modelling endeavors may also investigate if high epibiont cover promotes cannibalism in crabs. That is, under high epibiont cover of adult mussels, would a crab prefer to cannibalise its own conspecifics [46,47], rather than switching to juvenile mussels? Another interesting future direction is to look at the foraging of crabs as they move in and out of patches containing mussels, some of which might be protected by mussel farmers, akin to marine protected areas [9,10]. A spatially explicit approach to this end, modeling a changing habitat based on mussel density, would also be interesting [14].
The empirical literature shows that while epibionts alter the prey choice of predators, including crabs and sea stars [28,29,30,22], there are no prey switching experiments using D. vexillum. Our goal is to provide firm modeling grounds for the scope of such experiments in the future. Thus a logical next step for empirical studies is to conduct experiments with D. vexillum to confirm our switching hypothesis, as well as look at switching scenarios under varying levels of overgrowth (with both living and artificial epibionts such as in [28]). An interesting research question therein would be to ask if one sees the inverted parabola shaped curve, typical of OFT scenarios when measuring crab size versus mussel preference. If our switching hypothesis is confirmed, this should not be the case, as smaller size juvenile mussels should be preferred to adults under heavy epibiont cover. All in all we hope our results will help devise suitable strategies and measures that will enable a boost in dwindling mussel populations, particularly as new complexities arise in ecosystems, driven by rapid increase in invasions.

Acknowledgements
JL and RP would like to acknowledge valuable support from the NSF via DMS-1715377 and DMS-1839993.

Appendix
8.1. Optimal Strategy in our setting. Here, we give a rigorous reasoning for our switching hypothesis. If we following standard OFT, we can consider a fitness function We endeavor to maximize R(u 1 , u 2 ), the net rate of energy intake during foraging. The optimal strategy for a crab (according to classical OFT) relies on the density of mussels. That is for each (M A , M J ) , we get a set of optimal controls S(M A , M J ) known as the strategy map.
This is (129)-(132), which is actually a control system with controls (u 1 , u 2 ) relying on the state of the system. Now we look for controls belonging to the strategy map S(M A , M J ). Then we calculate the derivatives of S(M A , M J ) to investigate the maximizing controls u 1 and u 2 .
The sign of ∂R ∂u1 and ∂R ∂u2 depend on the e 1 h 2 − e 2 h 1 . A tricky point here is that attack rates depend critically on the density of adult and juvenile mussels. That is of (u 1 = 1, u 2 = 0), or (u 1 = 0, u 2 = 1) are feasible as attack rates if the mussel densities are above a certain density. However if M A , or M J fall below a certain critical level, theory predicts that the less preferred prey should also be attacked, and one might have a situation of (u 1 = 1, u 2 = 1). What we show next, is that if certain parametric restrictions are met, (u 1 = 1, u 2 = 0), or (u 1 = 0, u 2 = 1) are the only optimal choices for the crab, irrespective of mussel density. Proof. If e1 h1 > e2 h2 , ∂R ∂u1 > 0, the maximum of R(u 1 , u 2 ) is thusly achieved for u 1 = 1. And since the sign of ∂R ∂u2 does not depend on u 2 it follows if ∂R ∂u2 = 0, R(u 1 , u 2 ) will be maximized either with u 2 = 0 or u 2 = 1.Then we get the strategy map e1−d1h1 from the earlier stability calculations. We note, as long as d 1 > e2 h2 , and if this is enforced (u 1 = 1, u 2 = 0) is the only optimal strategy for the crab.
If e1 h1 < e2 h2 in order to maximize R(u 1 , u 2 ), we need u 2 = 1. The strategy map will switch to e2−d1h2 from the earlier stability calculations. We note, h1 , and if this is enforced (u 1 = 0, u 2 = 1) is again, the only optimal strategy for the crab.

and (94)
It follows from the Routh-Hurwitz stability criteria that all eigenvalues have negative real part if It is obvious that the first two conditions are always satisfied under feasibility condition (16). Furthermore, It is enough to solve bd 1 (ad 1 h 2 1 − ae 1 h 1 + d 1 δ 1 h 1 + δ 1 e 1 ) > 0, that is, e 1 − d 1 h 1 < d1δ1 a + d1δ1e1 ah1 . Therefore, the system (10)- (12) is asymptotically stable if 8.3. Proof of Theorem 3.3. The equilibrium state, of the system (129)-(132), for the epiboint is e = K. At the interior equilibrium state, the parameters u 1 = 0, u 2 = 1 and a(e) = a 2 . Since e will not effect the solution of C, M A and M J once u 1 , u 2 and a(e) are determined, then it is enough to nvestigate the following The Jacobian matrix about (C * , M * A , M * J ) is Since all the parameters are positive, it is obvious that J 22 < 0, J 23 > 0, J 31 < 0, J 32 > 0, and J 33 < 0. Under the feasibility condition (26), J 13 > 0. And the characteristic equation is given by By Routh Hurwitz stability criteria, all eigenvalues have negative real part if It is easy to check B 2 > 0 and B 0 > 0 under the feasility criterion (26). And To make J 22 J 33 −J 23 J 32 > 0, it is enough to show 4 bd1 δ1(e2−d1h2) d 1 δ 1 h 2 −ad 1 h 2 > 0, which gives us e 2 − d 1 h 2 < 16bd1δ1 a 2 . Furthermore, Since J 22 < 0 and J 33 < 0, 8.4. Proof of theorem 4.1. Now let a, the growth rate of juvenile mussels, as the bifurcation parameter. Therefore, if condition (16) holds, A 0 (a * ) are always positive.
and C * , M * A are given by (13) and (14). Furthermore, it is easy to verify that 8.5. Proof of theorem 5.4. The Hamiltonian of the system is given by We use the Hamiltonian to find a differential equation of the adjoint λ i , i = 1, 2, 3.
with the transversality condition gives as (121) λ 1 (T ) = λ 2 (T ) = λ 3 (T ) = 0 8.6. Proof of theorem 5.6. The Hamiltonian of our problem is given by The differential equations for λ 1 (t), λ 2 (t), λ 3 (t), are standard and are derived as in Theorem 5.4. the transversality condition is Considering ∂H ∂δ1 = −M A 2 λ 2 − δ 1 , we derive the optimal control for J 3 (δ 1 ) optimal control change with increasing handling time of adult mussels Figure 8. In this simulation we look at how u 2 changes w.r.t. h 1 .
Here h 2 = 1. We want to see the change in the control u 2 as h 1 increases. The control u 1 = 0 no matter how large h 1 is. What we notice is that u 2 suddenly goes down to 0 from 1, at a critical value h * 1 = 2.1.

Numerical Explorations of Bifurcations of alternate models
In the case e = K, we do not see a Hopf bifurcation numerically. It is worthwhile considering certain alternate models for the epibiont dynamics as future work. We motivate this via considering the following model, The only change here to crab-mussel system (129)-(132) is that we assume the carrying capacity of the epibiont is density dependent, and depends primarily on the adult mussel density that is K = K = k 1 M A + k 2 . Here k 2 represents alternate substrate that the epibiont can grow on. The four dimensional system has 11 parameters with four dependent variables. The following parameters are used in numerical simulations: The system evolve the stable limit cycles for parameter set 134. Time series for all species is shown in the figure 9a, while limit cycles in 2-D phase space are shown in fig 9b, 10a and 10b. To observe the more qualitative behavior of the model, one-parameter bifurcation diagram is drawn with respect to parameter d 1 and parameter a in the figures . A supercritical Hopf bifurcation occurs at d 1 = 0.3567 which emanates stable limit cycles. There is another supercritical hopf bifurcation at d 1 = 0.444. Between these two Hopf bifurcation, model has periodic solutions. After second Hopf bifurcation point model has stable solutions but crab populations are going to extinct. The dynamics is shown in one-parameter bifurcation diagram fig 11a. The qualitative dynamics has been also obtained for range of parameter a drawn in the fig 11b. Initially, for low parameter value a < 1.265 the crab population is too low but as parameter a increases, model exhibits stable coexistence. Further, it undergoes through supercritical Hopf bifurcation at parameter a = 3.386 which emanates stable limit cycles (green filled circle).
The parameter region has been obtained by drawing two-parameter (a, d 1 ) bifurcation diagram in the fig 12a. The parameter region for which one species goes extinction is shown by shaded region (extreme left), region for which stable coexistence is possible shown by red and region for which perioidc solution is possible shown by blue color in the diagram fig 12a. Another two-parameter (d 1 , δ) bifurcation diagram is drawn in the fig 12b. The parameter region for which stable coexistence occurs and region for which perioidc solution is possible is depicted in the diagram fig 12b. As shown by these bifurcation graphs, model has periodic solutions for biologically feasible choice of parameters and one can find the Hopf bifurctaion point for each of the parameters used in the model.
These results show that a Hopf bifurcation is possible, if one considers a density dependent carrying capacity for the epibiont. These results are robust in nature as (a) The biomasses of all species exhibited the periodic coexistence against the Time series is shown for parameter set 134. (b) A stable limit cycle in the Twodimensional phase space E, Mj for parameter set 134 Figure 10. Stable limit cycle in the 2-D pahase space plot different sets of parameters will yield the same qualitative behavior. The periodicity in the system is beneficial for harvesting and coexistence of all the species involved. (b) One-parameter bifurcation diagram with respect to parameter a Figure 11. One-parameter bifurcation diagrams to depict stablilty, Hopf Bifurcation point and periodic solutions with respect to parameter set 134.
(a) Two-parameters (a, d1) bifurcation diagram with respect to parameter set 134 (b) Two-parameters (d1, δ) bifurcation diagram with respect to parameter set 134 Figure 12. Two-parameter bifurcation diagrams to depict parameter region for the stable coexistence and periodic coexistence with respect to parameter set 134.