Designing the Next Generation of Condition Tracking and Early Warning Systems for Shellfish Aquaculture

Abstract

Early detection of environmental disturbances affecting shellfish stock condition is highly desirable for aquaculture activities. In this article, a new biophysical model-based early warning system (EWS) is described, that assesses bivalve stock condition by diagnosing signs of persistent physiological dysfunctioning. The biophysical model represents valve gape dynamics, controlled by active contractions of the adductor muscle countering the passive action of the hinge ligament; the dynamics combine continuous convergence to a steady-state interspersed with discrete closing events. A null simulation was introduced to describe undisturbed conditions. The diagnostic compares valve gape measurements and simulations. Indicators are inferred from the model parameters, and disturbances are assessed when their estimates deviate from their null distribution. Instead of focusing only on discrete events, our EWS exploits the complete observed dynamics within successive time intervals defined by the variation scales. When applied to a valvometry data series, collected in controlled conditions from scallops (Pecten maximus), the EWS indicated that one among four individuals exhibited signs its physiological condition was degrading. This was detected neither during experiments nor during the initial data analysis, suggesting the utility of an approach that quantifies physiological mechanisms underlying functional responses. Practical implementations of biological-EWS at farming sites are then discussed.

1. Introduction

Bivalves are ectothermic organisms that have colonized most aquatic and marine ecosystems. They are mainly sedentary species, live under a wide range of environmental conditions, and exhibit a variety of behaviors manifested by shell movements [1]. As a consequence, numerous authors have suggested bivalves would be good sentinels of aquatic environment quality (e.g., [2], and references therein) and this, in turn, inspired several long-term contaminant bioaccumulation programs using the shellfish species (e.g., [3]). However, many bivalve species are also extensively fished and cultivated [2,4]. This implies a need for a different type of monitoring that is better suited for evaluating shellfish stock growth and health within the constraints of aquaculture installations.

Since bivalves open and close their shells as part of their normal physiological functions and in response to stress, a large variety of valve movement monitoring approaches for natural [5] or farmed [6] environments have been proposed. Measurements of the valve movement frequency and amplitude (shell opening distance and how ‘wide’ the valves are open) have been conducted using various devices, since the early 20th century [7]. Valve gape dynamics are considered a good means to monitor environmental quality because valvometry measurements are both rapid and sensitive, and they provide a direct estimate of the health status of the organism [8]. Valve gape is also a cost-effective variable to measure. The signal has also been shown to integrate past living conditions that have influenced behavior and physiology [8,9,10]. Many studies have shown that different bivalve species respond to environmental factors. Variations in the valve opening amplitudes and movement frequencies have, thus, been attributed to changes in temperature for Crassostrea virginica [11], pH for several bivalve species [12,13], salinity for Mytilus edulis [14], the presence of toxic algae for Pinctada fucata [10] and Mytilus galloprovincialis [15], cadmium chloride concentrations in Corbicula japonica [16], contaminants of heavy metals (e.g., copper) in Corbicula fluminea [17], alkaline solutions (e.g., biofouling compounds) in Mytilus edulis [18], suspended particulate matter concentrations in Pecten maximus [19], and noises in Magallana gigas [20].

Although several studies [5,6,21] have shown the interest and feasibility of valvometry monitoring systems (see

Table A1. Bivalve alarm systems using valvometry mentioned in scientific literature. Systems marked with an asterisk (*) appear to have been commercialised. The reference given as a technical source contains the best technical description of the system, but is not necessarily the earliest publication.

Name	Developer or Owner	Technical Source
Dreissena-Monitor *	U. Cologne (Germany)	Borcherding 1992 [55]
Micrel	Ifremer (France)	Floch 1998 [56]
MosselMoniter *	Delta Consult B.V. (Netherlands)	Kramer and Foekema 2001 [5]
MolluSCAN eye *	U. Bordeaux (France)	Sow et al., 2011 [23]
MusselTracker	San Jose State U. (USA)	Miller and Dowd 2017 [57]
Sense-T	U. Tasmania and CSIRO (Australia)	Andrewartha et al., 2015 [6]

), Tran et al. [22] remarked that most efforts were focused on developing and improving valvometry technology and not data analysis. After this, a number of articles were published proposing new methods (see

Table A2. Analytical methods employed with, or proposed for, valvometry-based water quality monitoring systems. To the best of our knowledge, none have been implemented on a permanent basis within shellfish farms. The table lists the primary bivalve species studied, gives a brief description of calculational frameworks and the experimental observations used in calculations. All the articles are qualified by the developmental status of the systems: empirical studies report on tests within experimental and/or laboratory systems; data analysis means articles describe new methods of calculation applicable to valvometry data series; operational means the article describes a use-case.

Bivalve	Model or Statistical Framework (Measurements Used)	Status	Source
associated with MolluSCAN
review article	non-parametric model with kernel density estimates (valve gaping with experimental correlation to desired parameter)	operational	Andrade et al., 2016 [21]
Crassostrea gigas	estimated survival function with a semi-Markov model (least visited opening width)	data analysis	Azaïs et al., 2014 [27]
Oysters	estimates kernel density of valve gape state to estimate probability an individual will survive or not, post-observation period (opening amplitudes)	data analysis	Coudret et al., 2015 [24]
C. gigas	discrete signal processing of tail conditional probabilities (valve gaping and closing speed)	data analysis	Durrieu et al., 2015 [25]
C. gigas	same as Durrieu et al., 2015 (valve gaping, closing speed, and enzyme biomarkers)	data analysis	Durrieu et al., 2016 [47]
Corbicula fluminea	statistical comparison of behavioral and biochemical estimates (valve gaping, agitation, tissue contamination, and proteomics)	empirical	Miserazzi et al., 2020 [16]
C. gigas	non-parametric kernel regression, then correlation with parameter of interest (valve gape distances, movement speeds, and valve opening patterns per 24 h interval)	data analysis	Sow et al., 2011 [23]
C. fluminea	non-parametric kernel regression, then correlation with parameter of interest (valve gaping, siphon extension, and contaminant assays)	empirical	Tran et al., 2003 [22]
C. fluminea	logistic regression of valve response time as a function of contaminant dose (valve closure and contaminant detection threshold)	empirical	Tran et al., 2004 [58]
C. fluminea	similar to Tran et al., 2004 added maximum likelihood estimators (valve closure and contaminant detection threshold)	empirical	Tran et al., 2007 [59]
associated with MosselMonitor
Dreissena polymorpha; Mytilus edulis	comparison of descriptive statistics (average amount of time valves closed hourly and percentage of maximum opening distance)	operational	Kramer et al., 1989 [5]
Anodonta woodiana	comparison of descriptive statistics (average amount of time valves closed hourly and percentage of maximum opening distance)	operational	Giari et al., 2017 [33]
Other developments
C. gigas	suggests statistical correlation with environmental conditions (valve gape, heart rate, and body T)	empirical	Andrewartha et al., 2015 [6]
C. gigas, M. edulis	statistical correlation with arbitrary threshold set to 20% (synchronous movements and amplitude variance)	empirical	Bouget and Mazurié 1999 [60]
review article	suggests possible experimental correlation with ocean pCO${}_2$ (valve gaping behavior)	empirical	Clements and Comeau 2019 [12]
M. galloprovincialis	statistically-tested comparisons of periodicity of gaping and environmental parameters (relative valve openness)	empirical	Comeau et al., 2018 [61]
M. galloprovincialis	mixed-model ANOVA of valve gaping and algal toxin concentrations (valve gaping parameters)	empirical	Comeau et al., 2019 [18]
M. galloprovincialis	unnamed time-series analysis (averaged resistances from sensors classified as “not defective” by analog system)	empirical	Gnyubkin 2009 [49]
C. fluminea	3-parameter Hill equation fit of valve closing time series to construct dose-response profiles (fitted daily valve closing activity curves)	empirical, data analysis	Jou and Liao 2006 [26]
C. fluminea	see Jou and Liao 2006, with siphon movements added (fitted daily valve and siphon activity curves, as well as circadian rhythms)	empirical	Jou et al., 2016 [17]
C. fluminea	see Jou et al., 2006 (behavioral toxicity assays and valve daily activity)	empirical	Liao et al., 2009 [62]
M. californianus	statistical inference on estimates (valve gaping, body T, and posture)	empirical	Miller and Dowd 2017 [57]
M. edulis	statistical inference on estimates, linear regression (valve gaping and distance traveled in contaminant exposure tests)	empirical	Redmond et al., 2017 [8]
M. edulis	time series analysis of response curve and statistical inference (valve gape distance and closing)	empirical	Riisgård et al., 2006 [63]

). For example, Sow et al. [23] and Coudret et al. [24] proposed a statistical method based on a Gaussian kernel identification and smooth regression analysis. Durrieu et al. [25] and Jou and Liao [26] developed similar statistical analyses to detect changes in closing events frequency; and Azais et al. [27] described discrete transitions between four states of valve gape dynamics: two steady-states, opened and closed, and two transitory phases (i.e., opening and closing) using a Markov chain model.

However, all these methods are still phenomenological. Very few studies have tackled the problem of the quantification of the processes underlying the biological responses. We argue that without understanding processes that can explain observed behaviors, it is not possible to fully assess problems of bivalve conditions, because this means identifying dysfunctional patterns within highly variable individual dynamics. Among the most advanced reflections was the work performed by Jou et al. [17], who tried to explicit response functions. These authors struggled with false warnings in their controlled experimental conditions, and they also anticipated the problems with false positive or negative detection of dysfunction conditions within natural environments. The solution proposed was to increase the number of monitored factors, but this is unlikely to work without some sort of explicit function of the cause-to-consequence mechanisms [28]. A warning system system requires a precise definition of time and spatial scale (i.e., resolution), as well as an alarm criteria, in order to characterize changes of the observed dynamic regime [29] to work effectively.

The overall purpose of our study is to develop a new type of valvometry-based early warning system. This system aims to assess the physiological state of the organisms, while identifying local or global sources of disturbances. For this, we use a quantitative modeling framework designed to simulate valve gape dynamics as an indicator of their physiological state [30]. In this earlier study, we quantified bivalve gape dynamics based on biophysical processes that govern opening and closing of the shell. Simulations of this model were tested against the experimental results obtained from a species characterized as extremely sensitive (Pecten maximus, L. [31]). Our objective here is to use this model to elaborate a condition tracking and warning system for bivalve stocks in extensive aquaculture environments. We have revised its formulation to include controls by environmental factors, and a generic limitation function has been introduced. Then we combined sequential valvometry records with null simulations. We address the notion of synchronicity between individuals to be able to define survey units, as well as the problem of variation scales in time and space. Finally, we discuss how data can be assimilated in the model and how an early warning system (EWS) would decide if disturbances have begun to affect the stock condition.

2. Constructing a Model-Based Biological-EWS

2.1. Basic Definition of the Early Warning System for Bivalve Aquaculture

Extensive shellfish aquaculture consists of an immeasurable amount of individuals, spread over large areas. In contrast, most monitoring systems collect a very small number of data series, which are analyzed statistically to reveal anomalies that can be interpreted as environmental disturbances. The source of disturbances is often identified by correlating indicators of the anomaly with environmental variables that are collected locally or recovered from remote surveys [32]. EWS must have a specific operating mode, with the ability to analyze data in real time and to account for both the current state and past changes. Although many of the biological early warning systems (biological-EWS) only use thresholds (e.g., [33]), a complete EWS must be based on a forecasting trend [33,34]. Besides, many techniques were used to analyze data in the general domains of process modeling and machine learning, but none of them have attempted to simulate biological processes within a dynamic modeling framework. At present, this is the only means to forecast the physiological state of the animal under different environmental regimes. Ultimately, it is in the best interest of shellfish farmers to quantify the condition of shellfish, since this would increase their marketability [35].

However, several questions have arisen from the development of biological-EWS [32], highlighting three challenges that need to be addressed:

• first is to determine adequate temporal (defined by the individual dynamics) and spatial (defined by the environmental factors) scales of variations;
• second is to assess the inter-individual variability and consider how this source of variability can be filtered out to detect potential sources of disturbances;
• third is to provide precise and operational definitions of the practical components of the EWS. Setup, measurement protocols, maintenance, task flow, computer time, and outputs should be carefully established, in such a way that the system can be tested and be adapted to a large range of conditions.

Overall, except for specific mechanisms, a biological-EWS based on statistical analysis may reveal signs of disturbance, but does not permit the establishment of a diagnosis of the reasons for the physiological dysfunction [36]. To solve this problem, the biological-EWS should contain means to test hypotheses about the underlying mechanisms of the observed dynamics.

In our approach, the model and its implementation explore how to meet these challenges. The computer codes have been written in Scilab 6.1.0, a free and open source software for modeling and analysis (distributed under a GPL licence, and available at http://www.scilab.org (accessed on 29 September 2021)).

2.2. Valvometry to Infer Shell Valve Movements

Valvometry is the measurement of the valve gape distance in bivalve shells [7]. The same term may also refer to a set of different measurement techniques, from early “sooted glass” techniques [7] to automated devices, such as strain gauges [37], Hall-effect sensors (HES) [10], impedance electrodes [22], and fiber optic sensors [38]. Of all the sensors tested, HES have several distinct advantages. They make very high frequency measurements, are relatively inexpensive, and are small, light, and simple to install. A HES consists of a semi-conductor glued on one valve and a magnet glued on the other. Therefore, only one electrical cable links the HES to the data acquisition device. There are, however, some disadvantages specific to this type of sensor [1,30]:

• The technology of the sensor is based on the electromagnetic properties of the system composed of a magnet and the Hall-effect semi-conductor. It cannot be used if any other significant magnetic sources interfere with the system.
• The use of the HES system requires intensive calculations to treat the data series for two main reasons. First, the dynamic geometric property changes during valve movements complicates all of the calculations; second, the resulting inverse function that evaluates the valve gape distance from a series of Hall voltage measurements (see [1] for full details) implies finding the positive real roots of a polynomial function.
• The calibration step, which converts the Hall voltages into valve gape distances is difficult to assess before the data series of the Hall current variations are collected. The admitted procedure is to perform measurements, then to sever the adductor muscle(s) of the bivalve and calibrate the Hall voltage as a function of known distances. This is accomplished by inserting a series of wedges between the two valves [1]. However, in monitoring systems it is necessary to calibrate the HES on living organisms, prior to beginning survey measurements start. It can be performed by chemically relaxing the adductor muscle(s) [39] before inserting the calibrated wedges between the valves.
• As for any other electronic sensors, the sources of electrical noise [40] may affect the resulting measurements [41] and the electromagnetic properties are modified by environmental variability (e.g., temperature [42]). In addition, the precision of any valve gape distance estimates made in terms of sensitivity will decrease with the distance between the magnet and the semi-conductor [1].

Despite these disadvantages, the problems linked with sensor functioning can be addressed during the monitoring program with a site-specific calibration and data acquisition testing phase. For this reason, we have developed a specific protocol [1]. In the present article, we will not discuss these issues further, considering that whatever type of sensors and associated calibration function is used, estimates of the valve gape distances and angles are unbiased and sensitive enough to detect variations required for the monitoring program.

2.3. Modeling the Dynamics of Bivalve Shell Valve Movements

The basis of the model used in this study was described in a previous article ([30]; also see the summary in the Appendix A). It is the first biophysical model to represent processes involved explicitly in variations of the gape angle. Briefly, our model describes how the distance between the two valves, $\delta$, changes as a function of the angle (called the gape angle, $\alpha$, in radians) formed by two rays placed on each valve (Figure 1). The origin is located at the umbo, and the plan formed by the two rays is perpendicular to the axis materialized by the hinge. When the gape angle $\alpha =0$, the shell is closed. The gape angle is used because this parameter is less sensitive to changes in specific protocols, such as where the sensors are placed on individual organisms, and with differences in the sizes of organisms used. Any changes in $\alpha \left(t\right)$ are explained in terms of the passive forces, due to the ligament action, and the active forces (opening force ($F_O$) and closing force ($F_C$)) from the adductor muscle, which act to close the two valves. The adductor muscle(s) is anchored onto the shell interior some distance away from the valve edge and the hinge ligament. When the adductor muscle is active, it uses energy to shorten the muscle, pulling the valves closer. Closing the valves also means the muscle works against the elastic hinge ligament (stretching it). The valves open passively (the muscle relaxes and the hinge ligament returns to its shorter, unstretched state) with no additional energy input.

The primary dynamic model (Appendix A, Equation (A1)) is hybrid and simulates both continuous and discrete dynamics of the valve movements. The only factor applied to the dynamics is an instantaneous transition of the adductor muscle state. For the continuous dynamics, the existence of equilibrium solutions of the model depends on the existence of a balance between the passive and active components (see Appendix A).

For the passive component, due to the action of the ligament, the trivial equilibrium is reached when $\alpha \left(t\right)={\alpha }_{max}$, and this is the only stable steady state. Inversely, the active component, due to the action of the adductor model, has a non-trivial stable equilibrium solution, expressed as:

$$\left\{p_U^{*}=\frac{rs}{rs+as+ar},p_A^{*}=\frac{ar}{rs+as+ar},p_F^{*}=\frac{as}{rs+as+ar}\right\}$$

(1)

Combining these two stable state conditions allowed us to determine the value of the stable equilibrium angle:

$${\alpha }^{*}={cos}^{-1}\left(1-{\zeta }^2{\left(1-\sqrt{Q\frac{ar}{rs+as+ar}}\right)}^2\right)$$

(2)

where $\zeta =\sqrt{1-cos\left({\alpha }_{max}\right)})$.

${\alpha }_{max}$ (${\alpha }_{max}\in ]0,\pi /2\left[\right)$, the maximum opening angle, was measured, hence $\zeta$ was calculated. The parameters vector $\left\{Q,r,s,a\right\}\in {\mathbb{R}}^{+*4}$ must still be identified. Nevertheless, this calculation has shown that there is a steady state for which the effort performed by the contraction of the adductor muscle compensates for the passive force imposed by the ligament. This effort can be sustained, as long as the parameters r, s, and a (recovery, deactivation, and activation rates, respectively) remain constant. The equilibrium value depends on the fraction of fibers that are activated, $p_A^{*}$. Hence, based on the relative importance of the product $ar$, the activation rate is multiplied by the recovery rate. Conversely, the fraction of fibers that can be mobilized to perform fast closing, $p_U^{*}$, depends on the relative importance of the product, $rs$, multiplying the recovery rate by the deactivation rate. Finally, the equilibrium fraction of fibers that are de-activated, $p_F^{*}$ depends on the relative importance of the produce $as$, multiplying the activation rate by the deactivation rate.

Then, considering that there is a set of N environmental factors that can influence the adductor muscle, we hypothesize that a minimal standardized limiting function $f\left(\xi \right)$ regulates physiological rates within the muscle sub model. This function is expressed as:

$$f\left(\xi \right)=\prod _{i=1}^N{\xi }_i^{b_i}{e}^{b_i(1-{\xi }_i)}$$

(3)

where ${\xi }_i=(e_{max,i}-e_i)/(e_{max,i}-e_{opt,i})$ with $e_{max,i}$ is the upper tolerance limit for the environmental factor $e_i$, $e_{opt,i}$ is the optimal value at which $f\left(e_{opt,i}\right)=1$ and $b_i$ is a parameter that conditions the curvature of each controlling factor i. This means that there is no defined minimum value for the limitation function, even if there is a lower environmental boundary, below which the organism cannot ensure basic physiological requirements. In some cases, the limitation function gets very close to zero, in order to ensure continuity with the lethal condition of the lower environmental boundaries. In other cases, this continuity is not relevant (e.g., when a phase change, from liquid to ice occurs). For any maximum environmental boundary, $e_{max,i}$ (upper tolerance limit), the value of the limitation function $f\left(\xi \right)$ (Equation (1)) was set to zero. The domain of application of this function is, therefore, [−$\infty ,e_{max,i}$]. A multiplicative effect, as used in $f\left(\xi \right)$, is commonly used for eco physiological studies and is justified by the requirement that if one environmental factor places an individual in one of its physiological boundaries, then whatever the values of the other factors are, the physiological function cannot be ensured any longer (i.e., $f\left(\xi \right)\simeq 0$).

The muscle submodel (Appendix A, Equation (A3)) modulates three rates $\left\{a,s,r\right\}$, all in $tim{e}^{-1}$, where r is the muscle recovery rate, s is the fiber deactivation rate, and a is the fiber activation rate. In our biological-EWS, we consider that these rates could be affected by the environmental variability, hence deviate from optimal values for the physiology of the organism. The functions that modulates these rates as a function of environmental changes are then written as:

$$\{r=r_{opt}f\left(\xi \right);s=s_{opt}g\left(\xi \right);a=a_{opt}f\left(\xi \right)\}$$

(4)

where $g\left(\xi \right)=(\nu +f(\xi \left)\right(1-\nu \left)\right)$ with $\nu \ge 1$ is the amplification coefficient for the parameter s. The limitation value for a and r (the fibers’ activation and recovery rates, respectively) at the optimal environment value is equal to 1 and decreases to reach 0 when either the upper or lower environmental tolerance limit of the organisms is reached. For s, the fiber deactivation rate, the limitation function is reversed, with a minimal value equal to 1 at the optimal environment, and a maximum value defined by the amplification coefficient, $\nu$, reached at the environmental tolerance limits. An example of the response effect curves are presented (Figure 2). The functions $f\left(\xi \right)$ and $g\left(\xi \right)$ are asymetrical, leading to asymetrical rates of variations (expressed by the first derivatives), with, in the presented example, sharper effects above the optimal value than below. This led to not including any acclimation mechanism (i.e., optimal rates of the biophysical processes are fixed) and variations of the environment are considered over time frames that are much longer than the scale of the valve dynamics.

2.4. Dynamic Criteria Describing the Physiological Condition of Bivalves

In the quantitative framework described above, two categories of dynamic criteria describe the physiological condition of bivalves: continuous and discrete. The physiological condition is first characterized by the value of the non trivial equilibrium of the continuous dynamics, which quantifies the potential for the individual to maintain a partial opening. It is then characterized by the frequency and amplitude of the discrete events, which quantify the reaction capacity of the individual. It comes second, because the reaction capacity is conditioned by the continuous state, which determines what fraction of the fibers can be mobilized for discrete closing events.

Criteria are calculated from the dynamic model parameters, which are assumed to be attributes of the species. However, they are determined at the organism level, for each studied individual. Monitoring and biological-EWS systems implement replicates most of the time (e.g., [33]), but very few study the synchronization process of the discrete events and their control by external factors. This involves not only studying co-occurrences of the discrete event, but also the distribution of time between two events and the deviation from a null model. In our previous valvometry study with scallops [30], we found that under controlled constant conditions, the occurrence dynamics of single closing events could not be distinguished from a pure, stochastic process; there was an average occurrence rate of approximatively 2.4 h${}^{-1}$. This was confirmed by the fact that the observed probability distribution of time intervals between two events did not deviate significantly from an exponential law. This pattern defined the null model. Paradoxically, in this same study, the amplitude did not follow the same stochastic process dynamics, and instead seemed to have been controlled by the organisms. Even if the resulting distribution could not be characterized by any known analytical distribution law, the observed probability distribution deviated from any corresponding simulated one. However, the unexpected detection of a stochastic closing occurrence pattern suggests that fast closing events in bivalves are performed in the absence of synchronizers. This can be interpreted as if there is no internal synchronizer, then only external environmental variables can play such a role. Furthermore, the detection of external synchronizers in the individual valvometry series implies that several individuals, placed in the same conditions, would exhibit the same behavior. Thus, individual replicates should be implemented in any biological-EWS system.

2.5. Diagnosing Disturbances from Valvometry Measurements

Building on what we have learned in the previous section, a set of criteria can be defined to diagnose disturbances from valvometry data series:

• A sustained increase in the steady-state opening angle,
• A sudden and repeated increase in the closing peak frequency, and
• A significant increase in the intensity of closing peaks, up to extreme events of complete and sustained closure of the shell.

The first criterion indicates a decrease in physiological performance, hence an increasing physiological difficulty to maintain the adductor muscle fibers activated and contracted. The plots in Figure 3 show how these changes can occur under environmental control for a generic bivalve species. The maximum opening angle was standardized to 1 and the equilibrium value was arbitrarily defined (without losses of generality) between 0 and 1. The valve gape dynamics were then simulated, applying (in successive intervals of 24 h) a series of deviations from the optimal conditions for the adductor muscle contraction (fixed in the standardized interval [0, 1], here at 0.7). First the standardized environmental factors decreased from 0.7 to 0.35 (between 24 h and 48 h). The opening standardized angle increased concomitantly from 0.5 (equilibrium at optimal) to ca. 0.76. The amplitude of closing peaks increased, as well, from a maximum of approximately 0.3 to ca. 0.5, on average. Then, the standardized environmental factors decreased again from 0.35 to 0.01 between 48 h and 72 h (close to the lower boundary) and the opening angle at equilibrium increased from 0.96 degrees (equilibrium at optimal) to ca. 0.995 degrees (i.e., close to the maximum valve gape value). In addition, the amplitudes of the peaks decreased to very low values (ca. 0.021 on average).

Next, after a return to the optimal value between 72 h and 96 h, the standardized environmental factors increased from 0.70 to 0.85 between 96 h and 120 h. The standardized opening angle increased consequently from 0.50 to 0.78, and the amplitude of the closing peaks increased from 0.40 degrees to ca. 0.60. Finally, when the standardized environment factors increased from 0.85 to 0.99 between 120 h and 150 h. The standardized opening angle at equilibrium increased consequently from 0.79 degrees to 1.00 degrees and the peaks amplitude decreased down to values very close to 0. All fibers were deactivated at this stage. This simulation illustrates how the steady state opening angle (${\alpha }^{*}$) can be modulated by environmental changes. The minimum value for the steady state opening angle was achieved when the values of the environmental variables are optimal. However, measured data series contain a variability that could not be explained by deterministic control, hence should be considered and quantified when comparing simulations with observations. This is particularly true when estimating parameters because uncertainties are typically calculated by optimization when fitting curves to the data; classically, the sensitivity to detect disturbances is conditioned by the requirement that the amplitude of the estimators changes is greater than these uncertainties.

The simulation shown in Figure 3 demonstrates that continuous and discrete dynamic criteria are not independent, as conceived in the model framework. The maximum amplitude is reached for intermediate equilibrium opening angle between the minimal value (at optimal conditions), and the maximum opening angle, for which the closing amplitude tends to zero. Thus, the baseline condition is described by a set of inter-related parameters (formally $\left\{K,Q,r,s,a\right\}$) that defines both the continuous steady-state value and capacities to close the shell abruptly. Therefore, the three criteria are linked together. Any change in equilibrium modifies the possibilities and intensities of closing peaks, and vice versa. Furthermore, it would be difficult to discriminate changes in one of the criteria regarding the others from one individual data series. An objective to combine the three criteria therefore requires replication. In addition, as mentioned above, the identification of an environmental stress-induced behavior can only be validated if a set of several individuals, submitted to identical environmental conditions present the same valve movement pattern deviating from the null distribution assessed by the baseline condition. The number of individual replicates, however, does not have to be high. Indeed, for a purely stochastic process, the rate of synchronization between two individuals in a 10 s period was only 0.6% and nearly null for three individuals and more.

3. Working with Data

The primary purpose of the warning system is to assess the state of the shellfish stock by monitoring individual behaviors. The EWS proposed in his study is based on a process quantification, in which the formulation of the equations can be interpreted in term of biophysical dynamics. The dynamic model accounts for past states and can be used to predict future states explicitly. So far, our presentation of the development of the warning system has focused on a valve gape model without linking simulations with measurements. In essence, this is a virtual valvometry system that can be adapted to any bivalve species. The model aims to describe and quantify a known situation (i.e., all parameters are identified and refer to a specific condition of a simulated individual). In an ideal situation, under the hypothesis that physiological processes are not subject to acclimation, parameter values can be determined using a set of controlled experiments that cover the full range of environmental variability expected at a particular location. Then these parameters can be applied to the current situation, which accounts for environmental variables without any further parameter optimization. However, this ideal situation is not realistic and parameter identification still needs to be integrated into the monitoring process. The model should then allow for data assimilation (sensu lato), which not only helps to constrain the dynamics to an optimum trajectory (the deviation from the observations is an indicator of the predictability of the dynamic model) but would also provides sequential series of the optimal parameter estimates that can provide indicators of the (ecophysiological) state of the shellfish. By using these indicators, the critical values of the EWS can be defined.

3.1. Optimization in the Framework of Hybrid Dynamics

The first task of the optimization process was to develop an accurate method to detect peaks (occurrence, frequency, and amplitude). We applied a method to detect peaks, based on the variation of speed and turning points [30]. To verify the accuracy of this method, we simulated time series using the model, then applied the method and compared the results of detection with the discrete events time imposed in the simulations. The results are presented Figure 4. We allowed for random fluctuations, generated by a Wiener process, of the environmental control around the optimal value (upper panel of Figure 4a).

The lower left panel of Figure 4 shows that the empirical distribution of the time duration between two events (Figure 4d) follows an exponential law, close to the theoretical distribution, regardless of the fluctuation in the environment (Figure 4a). The distribution of peak amplitudes (Figure 4e) is dependent of the series of peak occurrence, hence on the duration between two peaks, and does not follow a priori any known theoretical distribution. Therefore, the reference distribution can only be seen as an empirical distribution characterizing the simulation. The estimated peaks show a good concordance with the simulations, and correlations between estimated and calculated occurrence times (Figure 4b) and between estimated and calculated amplitudes (Figure 4c) are both close to the linear equality. Except for few deviations, the occurrence of the peak (starting value) was well estimated. The estimated amplitudes (differences between the higher and lower value of the peaks) were slightly underestimated, the calculated values been systematically larger (Figure 4c). However, the differences never exceeded 5%.

The main challenge of the optimization process is that $\left\{Q,r,s,a\right\}$ cannot be determined from valvometry measurements alone [30]. Ancillary experimental data of direct measurements of the muscle activity and related forces that are generated should be measured for the species of interest. Even if forces can be estimated using dynamometer measurements [43], the method is very intrusive for individuals hence, it can only be installed for a very limited time and is not practical to deploy on many individuals for a monitoring program.

Therefore, it is necessary to consider simplifying the optimisation process performed from valvometry measurements only. The relative importance of the recovery rate, r, conditions how a hierarchy can be built. Let set a and s being of equal value. If r is much greater than a and s, then $p_A^{*}$ and $p_U^{*}$ will converge toward 0.5, while $p_F^{*}$ will converge to zero. In the best configuration, the activation rate is significant and the deactivation rate is negligible, leading $p_A^{*}$ to converge to 1, while $p_U^{*}$ and $p_F^{*}$ will converge to 0. The following hierarchy $\left\{r>a>s\right\}$ is suitable for the current problem of optimization. From this hierarchy, some transformations must be made to identify the vector $\left\{Q,r,s,a\right\}$. As Q (dimensionless) is a scaling factor between ${ϵ}^2$ and $p_A$, it can be set to 1 with no loss of generality. Then, $\left\{1,r,s,a\right\}$ can be standardized by stating that $\left\{a=k_{1}r,s=k_{2}r\right\}$. It follows that $k_2$ is linked to $k_1$ by:

$$k_2=\left(\frac{1}{{ϵ}^2}-1\right)\left(\frac{k_1}{k_1+1}\right)$$

(5)

hence the identification is restricted to the two parameter vector, $\left\{K,r\right\}$. The identification of $\left\{K,r\right\}$, and consequently $\left\{s,a\right\}$, was accomplished by optimizing the time series segments between two identified discrete peaks. A direct search algorithm [44] was used to minimize the sum of the square root differences between the simulated and observed valve gape angle values.

3.2. Using Valvometry Data and Series of Parameter Estimates to Detect Disturbances

Next, the warning system is constructed with a protocol defined from the identified criteria and the optimization procedure described above. Because of the hybrid nature of the dynamics, the monitored time series from a set of individuals must be segmented into a series of discrete successive sequences. Each sequence was then used to analyze and diagnose the state of the system by comparing observations with simulations. The definition of the duration of the discrete sequences was conditioned by the time scale of the dynamics (i.e., a time slot that encompasses the dynamics including both the continuous and discrete regimes). As the continuous regime was characterized by a convergence to a steady state between two discrete closing events, the minimum time scale was defined as about the inverse of half of the frequency of the peak occurrence (i.e., the Nyquist frequency).

To illustrate our approach, we used a data series presented in our two previous studies [1,30]. They corresponded to the Hall-Effect valvometry measurements that were collected in controlled conditions of temperature and salinity. A 12 h light-dark discrete cycle was imposed on the system. Four individuals (numbers 1 to 4) were observed for a total of 193 h (ca. eight days). Valvometry measurements were recorded at 10 Hz.

The average occurrence rate for the four individuals and for the recording period was estimated to be 2.6 ± 0.28 (SE) h${}^{-1}$. Thus, the duration of sequences were fixed at 1 h. In each time interval of one hour, K and r were estimated using optimization. The procedure evaluates each calculation value of a and s with a fixed Q (set equal to 1. $k_1$ and $k_2$ were re-estimated (see Equation (5)) at each iteration. Then, the steady-state value of the angle, ${\alpha }^{*}$ was calculated using Appendix A Equation (A6).

Two issues arose during the optimization procedure:

• When there was no peak, there was no means to detect a convergence to an equilibrium, hence it was impossible to estimate K, and the last value estimated was used for the simulation. This occurred 16, 37, 83, and 47 times over 193 time periods for individuals 1 to 4, respectively.
• When the equilibrium value (continuous baseline) drifted continuously (i.e., increased or, rarely, decreased monotonously), then the estimated values of $\{r,a,s\}$ tended to drift in such a way that the equilibrium estimate (Appendix A Equation (A6)) got closer to the extreme value, ${\alpha }_{m}ax$, or zero, respectively. Then, the equilibrium value was forced to be the average of the series and the parameter optimization was ignored. This occurred 22, 29, 17, and 24 times over 193 periods for individuals 1 to 4, respectively.

The results of four indicators are presented Figure 5. The upper panels (a) and (b), shows the direct continuous and discrete dynamic indicators. On the panel (a), the relative equilibrium angle exhibited stable values near 0.4 for individuals 2, 3, and 4, while it drifted significantly for individual 1, from 0.48 to 0.62. On the panel (b), changes in the number of peaks show that the highest record was reached for individual 2, with 13 peaks in one hour. The maximum number reached for the random occurrence with the observed frequencies, equal to 2.6 ± 0.28 (SE) h${}^{-1}$, and for 800 simulations (ca. four series of 193 sequences) was 12 peaks per hour. This value was reached only for individuals 2 and 3. The lower two panels (c) and (d) corresponds to the parameter estimates. Changes in K, panel (c), exhibit series of peaks for all individuals. Finally, changes in the ratio r/s (lower right, panel d) show that individual 1 was the only specimen to have decreased consistently below the threshold value of one. Hence, for individual 1, the the recovery rate (r) was consistently lower that the deactivation rate (s) after ca. period 30. For individual 3, the ratio r/s decreased initially to the threshold value, between ca. periods 30 and 90, then increased again above the threshold value and reached a level close to the values for individuals 2 and 4. By examining these 4 indicators, we remarked that, despite fluctuations, the condition of individual 1 differed from the conditions of individuals 2, 3, and 4 and there was no global synchronization among them.

4. Discussion

4.1. A System Based on Quantification of the Valve Gape Dynamics

Several monitoring and early warning systems, based on valvometry measurements, have been described in the literature (e.g., [8,23,45,46,47,48,49], see also Table A2); additionally, a few systems have been commercialized (Table A1). The system that we proposed in this study is, to the best of our knowledge, the first to be built on a biophysical model of the adductor muscle dynamics. This type of model offers an important advantage, compared to a system that only seeks statistical indicators of dysfunction. Our proposed approach constructs an explanatory biophysical framework for the patterns observed that can be linked to physiological observations. Furthermore, indicators are derived from quantification of the dynamics. The only other similar approach is that described by Gnyubkin [50], who proposed a biophysical model to simulate the movements of the valves, invoking the antagonist actions of the ligament and the adductor muscle. His model was designed on a single, second-order, differential equation that quantified the dynamics of the valve gape distance. Gnyubkin’s model includes three strictly positive parameters and a function, $P\left(t\right)$, interpreted as the outer force that sets the pendulum in motion (the contractile force of the adductors, which close/open valves and/or keep them in a fixed position). The main problem of this function is that if $P\left(t\right)$ is set to 0 (no contractile force), the stable equilibrium opening is equal to 0. To cope with this, Gnyubkin bounded $P\left(t\right)$ by a constant (“Const”) that fixes the distance of partial opening. In contrast, our model significantly improves on the biophysical framework by allowing the equilibrium opening as an emergent property of the muscle dynamics. This change is why optimization of a series of parameters of the muscle dynamics is possible and from which we can derive an indicator of dysfunction (e.g., the ratio between recovery and deactivation rates). In addition, the dynamic pattern is complex because it is submitted to multiple dynamic regimes (primarily a discrete regime that controls fast clapping movements and a continuous regime that sustains a stable partial opening state) that cannot be grasped with only one (set of) equation(s), nor one unique statistical model (e.g., [27]).

These dynamic regimes led to the definition of disturbance indicators:

• For a continuous regime, the equilibrium opening value is a balance between the passive opening force produced by the ligament and the active force contraction capacity of the adductor muscle. The wider the valve gape angle, the lower the contraction strength. Therefore, this is a relevant indicator of the physiological condition of the organism. As with all ectotherms, the metabolism is sensitive to temperature changes, but bivalves have developed physiological abilities to cope with unfavorable changes: primarily, the diversity of fibers allows different behaviors with minimizing the amount of energy [51]. An example of this is the catch behavior of the smooth muscle, allowing valve closing to avoid desiccation for bivalves in intertidal conditions, even if large temperature changes occur at the same time.
• For a discrete regime, which is mainly characterized by fast closings to expulse water and wastes and any irritating exogeneous particles or organisms. Bivalves also close their shell quickly when perceiving a threat. A behavior of repeated or sustained closing events is an indicator of stress. The deviation of the closing event probability from random process distributions is a good indicator and can lead to the definition of thresholds, but the identification of external disturbances relies on the detection of synchronized behaviors, hence requires replication. There is no standardization regarding replications and about the number of individuals that should present a synchronized behavior to detect an external disturbance. Surprisingly, there is no statistical estimate of such a number, even if the the analysis of the individual disturbed behavior heavily rest on statistics (e.g., [25]).
• A particular regime applicable only for mobile species (i.e., not attached to a substrate and able to swim over short distances) is characterized by a reflex-type series of contractions that create strong valve movements. For scallops, used as a model organism in this study, such activity is used to bury (to orient the opening toward the main flow direction), to jump, or to swim [1]. It causes sequences of slow increases of the valve gape angle followed by a sharp closing and then a slow reopening [1,30]. This requires a substantial amount of energy, hence these episodes are limited in time and tend not to occur if the organism is physiologically deficient. Such behaviors occur in case of punctual stress (e.g., threat by predators) but not exclusively and they may be difficult to exploit because of their scarcity.

The adductor muscle is the active component that counters the passive opening of shell valves by the ligament located at the umbo (Figure 1). The biophysical formulation is the most appropriate because it is straightforward and expresses forces that generate movements. It is, however, a terminal function and should be enriched by two fundamental principles of control. The first control is provided by the concentration of the calcium-troponin complex and the second control is stimulation by the nervous system of the organism. The model proposed in [52], developed for vertebrates, links these two controls, modulating the concentration of the calcium-troponin by the occurrence of successive stimulations. This model indicates that there is also an equilibrium force that can be achieved, a function of the concentration of the calcium-troponin, which is dynamically controlled by stimulations. The second control, crucial for discrete events occurrence, is performed by the nervous system of the organism. In the absence of ancillary data (i.e., the electrical activity of the nervous system), it is very difficult to introduce this control of the adductor muscle, hence at this time, the model is still restricted to a phenomenological description of discrete event occurrences.

Once the dynamic regimes were assessed, the sequences on which estimates could be performed were determined, in such a way that a complete set of dynamic criteria can be studied. Once again, the null model (no disturbances and a random occurrence of closing peaks) was used to define the sequence duration. The continuous dynamics were assumed to be close to a local equilibrium (defined by the environmental condition values), and the occurrence of the closing peaks was considered to follow a distribution law that corresponded to a pure stochastic process (Figure 4). Indicators can be calculated for each sequence as deviation from the null model. Thresholds are defined by the extreme values of the null distribution. Hence the threshold indicates when disturbances affect one or several component of the dynamics. The variability of the environment can be taken into account using the limiting functions described in Equations (1) and (2). Nothing guarantees that the optimal conditions are fulfilled, and the parameters of these functions can be estimated by successive iterations. Particularly, the minimum angle is achieved by and defines optimal parameters. However, in the description of the model, we mentioned that the two limiting functions $f\left(\xi \right)$ and $g\left(\xi \right)$ can be asymmetrical, inducing asymmetrical rates of variations, rising questions about the potential for acclimation, that were not taken into account in the present study, because of a lack of quantitative knowledge about the evolutionary ecology of bivalves.

In the experimental context, environmental conditions were controlled and not constraining for individuals (i.e., were not set close to tolerance limits ). In this context, our EWS diagnosed that individual 1 showed a decline in physiological condition when the other individuals did not, or at least, recovered after a short period of slight decline (e.g., individual 3). The deterioration of individual 1 was assessed mainly because the $r/s$ ratio decreased below the threshold of 1, indicating that the recovery does not compensate de-activation any longer. This result was not identified when analysing the outcome of the experimental work [1] or while constructing and simulating the biophysical dynamic model [30]. It demonstrates the relevance of our approach and the benefits an Early Warning System derived from combining valvometry measurements and a dynamic modeling of underlying physiological mechanisms.

4.2. Scoring the Dynamic Regimes

To facilitate the interpretation of the assessment methods, scoring system were introduced as formal approaches to support decision theory [53]. To interpret indicators, a scoring system was proposed, providing an image of how each individual performs in their environment, over a series of consecutive intervals. The simplest scoring system is to start from 0 and to count +1 if the indicator is deviating from favorable conditions, −1 if it converges to favorable conditions, and 0 if it does not change from one sequence to another. Therefore, such a scoring system is a qualitative expression of the discrete dynamics of parameter changes at the scale of the sequences. A simple process to shift the baseline is to reset negative values to zeros (which then converge to the optimal condition). The main disadvantage is that any random fluctuation modifies the score and then contributes to a possible drift of the system behavior. Figure 6 illustrates this by presenting the results of scoring the successive changes in equilibrium angle. In the upper panel (a) any increase in the equilibrium angle (overall muscle relaxation) is scored +1, any decrease, −1, and 0 otherwise. The score increased fast for individuals 1 and 3 and slower for individuals 2 and 4.

To cope with the variability, there are two options. The first is to integrate changes over a given past sequences window, but the size of the window can only be determined arbitrarily. The second option is to compare changes regarding the estimated variability of the pure random corresponding predictor. The score is modified only if the variation of the indicator if larger than the standard error of the corresponding estimates. The lower panel (b) of Figure 6 shows results of this procedure. Scoring that account for variability makes the drift less pronounced and the variations of individual 1 did not differ from the 3 others. A scoring system is a promising communication tool for our biophysical modeling framework and is complimentary to estimates and indicator calculations. However, it requires a separate, specific development to improve usability.

4.3. Practical Consideration about Data Acquisition and Alerts

Some biological-EWS are based on the detection of behavioral changes by organisms in situ (e.g., [23]) and others use semi-controlled experimental conditions (e.g., [33]). The advantage of the second setup is that it allows for the regulation of the variability by tempering the effect of environmental fluctuations. Our initial development [1,30] used controlled experiments to evaluate the sensor technology and to test simulations of the model we designed to represent the observed valve gape dynamics. We believe that this remains a preferred path for an alert system, since it can simplify the comparison between simulations and observations because it decreases the number of parameters describing the response to environmental variability.

The third challenge mentioned in (Section 2.1), stressing the need to provide precise and operational definitions of the practical components of the EWS, was not entirely addressed in this study. However, we can state that the instrumented group of individuals used in the alarm system should be placed in flow-through tanks, at different locations in the exploited area, or basin. This will respect the conditions of the definitions used in the model. Each unit must, at least, record basic environmental variables temperature, salinity, turbidity, and valve-gape distance variations on the set of selected individuals. The system should be conceived as an evolving tool, allowing it to be enriched with new measures as information about the dynamics at a particular site improves. Finally, another set of individuals constituting a background, or controlled, survey unit should be set-up in controlled environmental conditions, separated from the alarm system individuals.

5. Conclusions

This study presented a new approach for biological-EWS based on valvometry measurements and simulations of a biophysical model. We have shown how our system can distinguish environmental disturbances from physiological weaknesses using the synchronization of responses among a group of individuals kept in identical conditions. To characterize undisturbed situations, we introduced a null model where discrete events are random and physiological parameters are controlled by variations of environmental factors. From this modeling, we inferred four disturbance indicators: (1) the equilibrium angle, (2) number of peaks, (3) amplitude of closing peaks, and (4) the muscle fiber recovery-to- deactivation rates ratio. Disturbances are assessed when indicator values deviate from their null distributions. The biological-EWS described here is advanced enough to be tested in a farming site configuration. At this stage, we anticipate that our framework would evolve in collaborative platforms with different stakeholders involved in improving the conditions of shellfish culture.