## 4 Simulation and Analysis  
In this section, some simulations of the system (1)-(3) are performed to illustrate the results stated in Theorems 1 and 2. We consider the following reaction rate expressions:  
\[r_{s1}=\mu_{\text{max}}\frac{S_{1}}{K_{1}+S_{1}},\,\text{hence}\,\,r_{s1}=\mu_ {\text{max}}\frac{x_{1}+S_{1}^{\text{eq}}}{K_{1}+x_{1}+S_{1}^{\text{eq}}}r_{s2 }=\mu_{\text{max}}\frac{S_{2}}{K_{2}+S_{2}},\,\text{hence}\,\,r_{\alpha}=\mu_{ \text{max}}\frac{x_{2}+S_{2}^{\text{eq}}}{K_{2}+x_{2}+S_{2}^{\text{eq}}}.\]The terms \(r_{s1}\) and \(r_{s2}\) satisfy assumption 3. The model (1)-(3) with these reaction rate expressions was calibrated with data of suspended solids concentration from vertical subsurface flow constructed wetlands (CWs) reported in [37]. The CWs are located in Cuenca (Ecuador) and receive domestic wastewater coming out from primary treatment. The influent wastewater comprises a suspended solids concentration of 88.83 mg/L, BOD\({}_{5}\) concentration of 95.75 mg/L, and temperature of 24.6 \({}^{\circ}\)C. The hydraulic loading rate (HLR) is 0.2 md\({}^{-1}\). Model (1)-(3) was fitted by minimizing the sum of the square residuals between experimental and model data of suspended solids concentration. The parameter values obtained were used in the first simulation case and are shown in Table 1.  
To evaluate the obtained convergence results stated in Theorems 1 and 2, model (1)-(4) is simulated with  
\[\overline{S}_{in}=100\;\text{mg}/\text{L};\;\delta_{sin}=A_{1}sin\left(\frac{2 rt}{T_{sin}}\right),\;A_{1}=10;\;T_{sin}=20\;\text{days}. \tag{118}\]  
The bounds \(S_{1}^{l}\), \(S_{1}^{u}\) of \(\Omega_{S1}\) (47); the bounds \(S_{2}^{l}\), \(S_{2}^{u}\) of \(\Omega_{S2}\) (114); and the bounds \(S_{3}^{l}\), \(S_{3}^{u}\) of \(\Omega_{S3}\) (116) are computed and represented by the upper and lower horizontal dotted lines, whereas the equilibrium values \(S_{1}^{eq}\), \(S_{2}^{eq}\) and \(S_{3}^{eq}\) (Equations (5)-(7)) are represented by the middle horizontal dotted lines. The calculation of the aforementioned bounds requires the calculation of bounds \(x_{1}^{l}\), \(x_{1}^{u}\) (24), (25); \(x_{2}^{l}\), \(x_{2}^{u}\) (70), (71); and \(x_{3}^{l}\), \(x_{3}^{u}\) (97). Three simulations are performed, using definitions (118) and parameter values shown in Table 1.  
In the three simulations, convergence to the regions within the calculated bounds is observed: the state \(S_{1}\) convergences to the region within computed bounds \(S_{1}^{l}\), \(S_{1}^{u}\) (see Figure 2a, Figure 3a, Figure 4a,b); the state \(S_{2}\) converges to the region within computed bounds \(S_{2}^{l}\), \(S_{2}^{u}\) (see Figure 2b, Figure 3b, Figure 4c,d); and \(S_{3}\) converges to the region within the computed bounds \(S_{3}^{l}\), \(S_{3}^{u}\) (see Figure 2b, Figure 3b, Figure 4c,d). \begin{table}
\begin{tabular}{c c c c} \hline \hline \multirow{2}{*}{**Parameter**} & **Value for the First** & **Value for the Second** & **Value for the Third** \\  & **Simulation Case** & **Simulation Case** & **Simulation Case** \\ \hline \(Q_{1}/v_{1}\) & 0.702 day\({}^{-1}\) & 0.702 day\({}^{-1}\) & 0.702 day\({}^{-1}\) \\ \(Q_{2}/v_{2}\) & 0.702 day\({}^{-1}\) & 0.702 day\({}^{-1}\) & 0.702 day\({}^{-1}\) \\ \(Q_{D}/v_{2}\) & 0.0001 day\({}^{-1}\) & 0.0001 day\({}^{-1}\) & 0.0001 day\({}^{-1}\) \\ \(Q_{D}/v_{3}\) & 1 day\({}^{-1}\) & 1 day\({}^{-1}\) & 1 day\({}^{-1}\) \\ \(\mu_{max1}\) & 101 mg /(L day) & 101 mg /(L day) & 101 mg /(L day) \\ \(K_{1}\) & 303.1 mg/L & 303.1 mg/L & 303.1 mg/L \\ \(\mu_{max2}\) & 101 mg/(L day) & 202 mg/(L day) & 404 mg/(L day) \\ \(K_{2}\) & 303.1 mg/L & 151.55 mg/L & 30.31 mg/L \\ \hline \hline \end{tabular}
\end{table} 
Table 1: Parameter values for the CSTR network model (1)–(3). Figure 3: Simulation of the CSTR network model (1)–(4), second case. (**a**) Time course of \(S_{1}\); the upper and lower horizontal dotted lines indicate the bounds of \(\Omega_{S1}\).