Model

The Accumulation of MPs in human body

The lifetime accumulation of microplastics for human can be modeled as the rate of the amount of microplastics inside human's bodies respect to the time, in which the balance of intake and loss processes are considered.

We have:

$$\begin{align*} \frac{\mathrm{d} C_{\mathrm{MP}}}{\mathrm{d} t} = (1 - f_{\mathrm{abs}} \cdot f_{\mathrm{a, n}}) \sum_{i = 1}^{n}(\mathrm{Ir}_{i} \cdot C_{i}) + (1 - f_{\mathrm{abs}}) \cdot (f_{\mathrm{dep}} \cdot \mathrm{InR} \cdot C_{\mathrm{a}}) \\ + k_{\mathrm{tis}} \cdot C_{\mathrm{tis}} - k_{\mathrm{loss}} \cdot C_{\mathrm{MP}} \end{align*}$$

and the amount of microplastics inside human's tissue can be expressed as:

$$\frac{\mathrm{d} C_{\mathrm{tis}}}{\mathrm{d} t} = f_{\mathrm{abs}} \cdot f_{\mathrm{a, n} } \sum _{i = 1}^{n}(\mathrm{Ir}_{i} \cdot C_{i} ) + f_{\mathrm{abs} } \cdot (f_{\mathrm{dep} } \cdot \mathrm{InR} \cdot C_{\mathrm{a}} ) - k_{\mathrm{tis} } \cdot C_{\mathrm{tis} }$$

where $i$ refers to the index of the $n$ food types in the study, $C_{i}$ is the microplastics concentration per media in terms of number concentrations, $\mathrm{Ir}_{i}$ represents the related daily ingestion rate, $f_{\mathrm{abs}}$ represents the intestinal absorption, $f_{\mathrm{a, n}}$ refers to the total number of ingested particles ranged from 1 to 10 $\mu$m, $C_{\mathrm{a}}$ is the microplastics concentration in the air, $\mathrm{InR}$ is the inhalation rate, $f_{\mathrm{dep}}$ is assumed as a constant representing the fraction of inhaled particles that deposited in the nasopharyngeal cavities, $C_{\mathrm{tis}}$ refers to the amount of removal of microplastics from the tissue, $k_{\mathrm{tis}}$ is the biliary excretion rate constant, and $k_{\mathrm{loss}}$ represents the loss rate constant of microplastics based on stool frequencies.

Collectively, these equations and the relevant data help us constructed a well simulation and visualization of microplastics accumulation over lifetime inside human's bodies.

Protein Structure Prediction

Protein Docking

To predict the three-location structure of the filtration matrix, we applied Protein-Protein Docking.

The acceptor of Docking is oat protein and the ligand is protein on BCM. For the energy fraction calculation, we used the equilibrated weights.

$$E = 0.40 E_{rep} - 0.40 E_{att}+600 E_{elec}+1.00 E_{DARS}$$

The results of Docking are shown below.

The first ranked protein fraction after binding was accepted for thermal stability simulation by molecular dynamics.

AlphaFold2

In order to obtain the predicted three-dimensional structures of the proteins in our membrane, we used a deep learning-based method -- Alphafold. The Alphafold method generally includes six steps, which are described below.

First, the Alphafold accepts the amino acid sequence of its primary input. Meanwhile, it also derives multiple sequence alignments of the related proteins, which provides the evolutionary information of similar proteins for the algorithm.

Second, the algorithm will build up the neural network architecture, which is based on transformer networks and convolutional layers. The function of these neural networks is to capture the long-range dependencies and relationships between amino acids in the protein sequence, which is useful in predicting the 3-D structure.

Third, the Alphafold will train the data by a large and diverse data set, which includes protein structures from the Protein Data Bank (PDB) and corresponding sequence data, and by observing the data, the model will learn to correlate the sequences with the corresponding protein structures.

Fourth, the Root Mean Square Deviation of $\text{C}\alpha$ Atoms ($\text{C}\alpha$ RMSD), which is used to measure the dissimilarity between the predicted positions and the real positions from the experimental data, will be left during the training of Alphafold. During the optimization, Alphafold will adjust the parameters of the model automatically.

Fifth, after the training process, Alphafold will predict the 3-D coordinate of a protein's atom based on its amino acid sequence and the evolutionary information. Meanwhile, the model uses an attention mechanism to identify the interactions between amino acids and predict the distance between them. Alphafold also has a "refinement" step where the algorithm will improve the initial predictions to enhance the accuracy of the final structure.

Finally, the model will verify the accuracy of the final predictions by several different metrics, including Global Distance Test (GDT) scores, TM-scores, etc. Such metrics quantify how well the predicted structures match the experimental data.

The above figure shows the predicted structure of the protein. The Alphafold gave us a relatively accurate structural document of our protein, which assists future use.

Simulations

Thermal Stability

In order to test the thermal stability of the BC-APT membrane, a quantum chemistry approach called ReaxFF was used to simulate the stability of the membrane at 400 Kelvin, the results of which are shown below:

Specifically, the ReaxFF method models chemical reactions with atomistic potentials based on a range of reactive force fields. Similar to the nonreactive force field, the reactive force field categorized the energy of the whole system into several different partial energy contributions, as shown below:

$$\begin{align*} E_{\text{system}} = E_{\text{bond}} + E_{\text{over}} + E_{\text{under}} + E_{\text{val}} + E_{\text{pen}} + E_{\text{tors}} +\\ E_{\text{conj}} + E_{\text{vdWaals}} + E_{\text{Coulomb}} \end{align*}$$

Moreover, the basic assumption of ReaxFF is that the bond order $\text{BO}_{ij}'$ between two atoms can be calculated by the interatomic distance $r_{ij}$ directly, as shown in the following equation:

$$\text{BO}_{ij}' = \text{exp}\left [ p_{bo, 1} \cdot (\frac{r_{ij}}{r_{o}} )^{p_{bo, 2}} \right ] + \text{exp}\left [ p_{bo, 3} \cdot (\frac{r_{ij}^{\pi}}{r_{o}} )^{p_{bo, 4}} \right ] + \text{exp}\left [ p_{bo, 5} \cdot (\frac{r_{ij}^{\pi \pi }}{r^{o}} )^{p_{bo, 6}} \right ]$$

Where $p_{bo, 1}$ and $p_{bo, 2}$ represent the exponential terms of sigma term with the unity below $\sim 1.5 \text{\AA}$ and negligible above $\sim 2.5 \text{\AA}$, $p_{bo, 3}$ and $p_{bo, 4}$ mean the exponential terms of the first pi bond with the unity below $\sim 1.2 \text{\AA}$ and negligible above $\sim 1.75 \text{\AA}$, and $p_{bo, 5}$ and $p_{bo, 6}$ represent the exponential terms of the second pi bond with the unity below $\sim 1.0 \text{\AA}$ and negligible above $\sim 1.4 \text{\AA}$.

The bond orders $\text{BO}_{ij}'$ should be corrected for overcoordination, and the final bond orders $\text{BO}_{ij}$ in the molecule are acquired by multiplying $\text{BO}_{ij}'$ to some correction factors. Hence, the following equation is used to calculate the bond energies with the bond order $\text{BO}_{ij}$:

$$E_{\text{bond}} = -D_{\text{e}} \cdot \text{BO}_{ij} \cdot \text{exp} \left [ p_{\text{be, 1}}(1 - \text{BO}_{ij}^{p_{\text{be, 1}}}) \right ]$$

Where $D_{e}$ and $p_{be, 1}$ are bond parameters involved in $\text{C}-\text{C}$, $\text{C}-\text{H}$, and $\text{H}-\text{H}$ bonds.

To control the overcoordination of the bond orders, $E_{\text{over}}$, which will converge to 0 rapidly for undercoordinated system, is used to punish this phenomenon.

$$E_{\text{over}} = p_{\text{over}} \cdot \Delta_{i} \cdot \left ( \frac{1}{1 + \text{exp}(\lambda_{6} \cdot \Delta_{i})} \right )$$

where $p_{\text{over}}$ is a constant parameter and $\Delta_{i}$ is the degree of deviation of the sum of the corrected bond orders around an atomic center from its valency $\text{Val}_{i}$ (normally $\text{Val}_{i} = 4$ for carbon and $\text{Val}_{i} = 1$ for hydrogen), herein, $\Delta_{i} > 0$. Thus, we have:

$$\Delta_{i} = \sum_{j = 1}^{n\text{bond} } \text{BO}_{ij} - \text{Val}_{i}$$

The valence energy $E_{\text{val}}$ is used to calculate the energy contribution for the $\pi$-electron between the centers of undercoordinated atoms.

$$E_{\text{under}} = -p_{\text{under}} \cdot \frac{1-\text{exp}(\mu_{1} \cdot \Delta_{i})}{1 + \text{exp}(-\mu_{2} \cdot \Delta_{i})} \cdot f_{1}(\text{BO}_{ij,\pi}, \Delta_{j})$$ $$f_{1}(\text{BO}_{ij,\pi}, \Delta_{j}) = \frac{1}{1 + \mu_{3} \cdot \text{exp}(\mu_{4} \cdot \sum_{j = 1}^{\text{neighbors(i)} } \Delta_{j} \cdot \text{BO}_{ij,\pi } )}$$

The penalty energy $E_{\text{pen}}$ is used to calculate the energy contribution for the $\pi$-electron between the centers of undercoordinated atoms.

Moreover, the following equations are developed to clarify and calculate the valence angle energy contribution.

$$E_{\text{val}} = f_{2}(\text{BO}_{ij}) \cdot f_{2}(\text{BO}_{jk}) \cdot f_{3}(\Delta_{j}) \cdot \left \{ k_{a} - k_{a} \text{exp} \left [ -k_{b}(\Theta_{o} - \Theta_{ijk})^{2} \right ] \right \}$$ $$f_{2}(\text{BO}_{ij}) = 1 - \text{exp}(-\mu_{5} \cdot \text{BO}_{ij}^{\mu_{6}})$$ $$\begin{align*} f_{3}(\Delta_{j}) = \frac{2 + \text{exp}(-\mu_{7} \cdot \Delta_{j})}{1 + \text{exp}(-\mu_{7} \cdot \Delta_{j}) + \text{exp}(p_{v, 1} \cdot \Delta_{j})} \cdot \\ \left [ \mu_{8} - (\mu_{8} - 1) \cdot \frac{2 + \text{exp}(\mu_{9} \cdot \Delta_{j}) }{1 + \text{exp}(\mu_{9} \cdot \Delta_{j}) + \text{exp}(-p_{v, 2} \cdot \Delta_{j}) } \right ] \end{align*}$$ $$\Theta_{o} = \pi - \Theta_{0, 0} \cdot \left \{ 1 - \text{exp}\left [ -\mu _{10} \cdot (2 - \text{SBO}2 ) \right ] \right \}$$ $$\text{SBO} = \Delta_{j} - 2 \cdot \left \{ 1 - \text{exp}\left [ -5 \cdot (\frac{1}{2} \Delta_{j} )^{\mu_{11}} \right ] \right \} + \sum_{n = 1}^{\text{neighbors}(j) } \text{BO}_{jn, \pi }$$ $$\Delta_{j,2} = \Delta_{j} \text{ if } \Delta_{j} < 0$$ $$\Delta_{j,2} = 0 \text{ if } \Delta_{j} \ge 0$$ $$\text{SBO}2 = 0 \text{ if SBO} \le 0$$ $$\text{SBO}2 = \text{SBO}^{\mu_{12}} \text{ if } 0 < \text{SBO} < 1$$ $$\text{SBO}2 = 2 - (2 - \text{SBO})^{\mu_{12}} \text{ if } 1 < \text{SBO} < 2$$ $$\text{SBO}2 = 2 \text{ if SBO} > 2$$

where $\Theta_{ijk}$ represents the deviations in valence angle from its equilibrium value $\Theta_{o}$, the $f_{2}(\text{BO})$ is used to guarantee the valence angle energy contribution will converge and disappear gradually during bond dissociation, the $f_{3}(\Delta_{j})$ copes with the phenomena of over/undercoordination in the central atom $j$, the values of $\Theta_{o}$ and $\Theta_{ijk}$ rely on the sum of $\pi$-bond orders (SBO) around the central atom $j$ on the valence angle energy, and $\mu_{5}$, $\mu_{6}$, $\mu_{7}$, $\mu_{8}$, $\mu_{9}$, $\mu_{10}$, $\mu_{11}$, $\mu_{12}$ are parameters of valence angle energy, which are 1.49, 1.28, 6.30, 2.72, 33.87, 6.70, 1.06, and 2.04 respectively.

For the stability of the system when dealing with double bonds sharing an atom at a valency angle, the following equations of energy penalty are developed for such systems.

$$E_{\text{pen}} = \mu_{13} \cdot f_{4}(\Delta_{j}) \cdot \text{exp}\left [ -\mu _{14} \cdot (\text{BO}_{ij} -2)^{2}\right ] \cdot \text{exp}\left [ -\mu_{14} \cdot (\text{BO}_{jk} -2)^{2} \right ]$$ $$f_{4}(\Delta_{j}) = \frac{2+\text{exp}(-\mu_{15} \cdot \Delta_{j})}{1 + \text{exp}(-\mu_{15} \cdot \Delta_{j}) + \text{exp}(\mu_{16} \cdot \Delta_{j})}$$

where $f_{4}(\Delta_{j})$ copes with over/undercoordination in central atom $j$ on the penalty energy, and $\mu_{13}$, $\mu_{14}$, $\mu_{15}$, $\mu_{16}$ are the parameters of penalty energy, which are 36.0, 7.98, 0.40, and 4.00 respectively.

To ensure the dependence of the energy of torsion angle $\omega_{ijkl}$ works accurately for bond orders close to 0 and bond orders larger than 1, the following equations were developed.

$$E_{\text{tors}} = f_{5}(\text{BO}_{ij}, \text{BO}_{jk}, \text{BO}_{kl}) \cdot \sin{\Theta_{ijk}} \cdot \sin{\Theta_{jkl}} \cdot$$ $$\left [ \frac{1}{2} V_{2} \cdot \text{exp}\left \{ p_{l} (\text{BO}_{jk} - 3 + f_{6}(\Delta_{j}, \Delta_{k}) )^{2} \right \} \cdot (1 - \cos{2\omega _{ijkl}}) + \frac{1}{2} V_{3} \cdot (1 + \cos{3\omega _{ijkl}}) \right ]$$ $$f_{5}(\text{BO}_{ij}, \text{BO}_{jk}, \text{BO}_{kl}) = \left [ 1 - \text{exp}(-\mu_{17} \cdot \text{BO}_{ij} ) \right ] \cdot$$ $$\left [ 1 - \text{exp}(-\mu _{17} \cdot \text{BO}_{jk} ) \right ] \cdot \left [ 1 - \text{exp}(-\mu _{17} \cdot \text{BO}_{kl} ) \right ]$$ $$f_{6}(\Delta_{j}, \Delta_{k}) = \frac{2 + \text{exp}\left [ -\mu_{18} \cdot (\Delta_{j} + \Delta_{k}) \right ] }{1 + \text{exp}\left [ -\mu _{18} \cdot (\Delta_{j} + \Delta_{k}) \right ] + \text{exp} \left [ \mu_{19} \cdot (\Delta_{j} + \Delta_{k}) \right ] }$$

where $V_{2}$ and $V_{3}$ are two sets of terms that are based on the bond order of the central bond $\text{BO}_{jk}$, the value of $\sin{\Theta_{ijk}} \cdot \sin{\Theta_{jkl}}$ guarantees that the torsion energy contribution will vanish when either $\Theta_{ijk}$ or $\Theta_{jkl}$ reaches $\pi$, the $f_{5}(\text{BO}_{ij}, \text{BO}_{jk}, \text{BO}_{kl})$ allows the gradual disappearance of the torsion energy contribution when one of the torsion angle bonds is dissociated, the $f_{6}(\Delta_{j}, \Delta_{k})$ prohibits the occurrence of excessive torsion energy contributions when atom $j$ and $k$ are overcoordinated, and $\mu_{17}, \mu_{18}, \mu_{19}$ are parameters of torsion energy, which are 3.17, 10.00, and 0.90 respectively.

The following equations depict the contribution of conjugation effects toward the molecular energy.

$$E_{\text{conj}} = f_{7}(\text{BO}_{ij}, \text{BO}_{jk}, \text{BO}_{kl}) \cdot \mu_{20} \cdot \left [ 1 + (\cos^{2}{\omega_{ijkl}} -1) \cdot \sin{\Theta_{ijk}} \cdot \sin{\Theta_{jkl}} \right ]$$ $$f_{7}(\text{BO}_{ij}, \text{BO}_{jk}, \text{BO}_{kl}) = \text{exp}\left [ -\mu_{21} \cdot (\text{BO}_{ij} - \frac{3}{2} )^{2} \right ] \cdot$$ $$\text{exp}\left [ -\mu_{21} \cdot (\text{BO}_{ij} - \frac{3}{2} )^{2} \right ] \cdot \text{exp}\left [ -\mu_{21} \cdot (\text{BO}_{kl} - \frac{3}{2} )^{2} \right ]$$

where $\mu_{20}, \mu_{21}$ are parameters of conjugation energy, which are -1.14, and 2.17 respectively.

Aside from the overlapping, intermolecular forces such as dipole-dipole attractions and dispersion forces should be taken into consideration and such forces are included in the categories of van der Waals and Coulomb forces. To prevent these two forces from disturbing the energy contribution during bond dissociation, the following equations are developed

$$E_{\text{vdWaals}} = D_{ij} \cdot \left \{ \text{exp}\left [ \alpha _{ij} \cdot (1 - \frac{f_{8}(r_{ij})}{r_{\text{vdW}}} ) \right ] -2 \cdot \text{exp} \left [ \frac{1}{2} \cdot \alpha _{ij} \cdot (1 - \frac{f_{8}(r_{ij})}{r_{\text{vdW}}} ) \right ] \right \}$$ $$f_{8}(r_{ij}) = \left [ r_{ij}^{\mu_{22}} + (\frac{1}{\mu_w{}} )^{\mu_{22}} \right ]^{\mu_{22}}$$

where $D_{ij}$ is the bond parameter, the $f_{8}(r_{ij})$ provides a shielded interaction to avoid repulsion and atoms sharing valence angle, and $\mu_{22}$ is the parameter of van der Waal energy, which is 1.69.

As with the van der Waals interactions, the following equation describes the Coulomb interactions between every atom pairs

$$E_{\text{Coulomb}} = C \cdot \frac{q_{i} \cdot q_{j}}{\left [ r_{ij}^{3} + (1\gamma_{ij})^{3}\right ]^{1/3 } }$$

where $\gamma_{ij}$ is the parameter in the Electron Equilibrium Method.

After defining every energy contribution to the force field, the whole system is optimized by an approach from van Duin et al, and related quantum chemistry data is added to the force field training set. During the optimization, the energy of every molecule is minimized continuously, and the quantum chemistry data is kept fixed with appropriate bond length or torsion angle restraints.

In order to control the amount of heat in the system, we used Berendsen thermostats to reach the demanded temperature. This method leads the system to the target temperature by changing the velocities of every atom during the simulation, allowing a rapid convergence in a relatively small and constant time interval.

The Berendsen thermostat uses the following scaling factor inside its algorithm to indicate the proportional scaling of the velocities per time step.

$$\lambda = \left [ 1 + \frac{\Delta t}{\tau _{b}} (\frac{T_{0}}{T} - 1) \right ]^{1/2}$$

where $\Delta t$ is the time step, $\tau_{B}$ is the time constant in the Berendsen thermostat, $T_{0}$ is the desired temperature, and $T$ is the instantaneous temperature. In our case, our time constant is 100 fs (strong damping and rapid convergence).

As shown in Figure 1 above, at temperatures near its boiling point, our protein exhibits a commendable level of stability, demonstrating a resistance to significant denaturation. This remarkable trait enables it to maintain its structural integrity and functionality under the challenging conditions of high temperatures. At the same time, this relatively stable and clear simulation method gives strong evidence for the thermal stability of our membranes.

Prediction of Filtration and Absorption Characteristics

For the simulation of BC-APT membrane permeation and adsorption efficiencies, we used a method based on Newtonian electromagnetic physics -- the classical force field.

Compared to ReaxFF, the Force Field method does not involve the usage of quantum chemical data and the calculation of bond orders, instead, it is a relatively fast and simple approach under Newtonian physics. However, it is a highly efficient and accurate method for predicting the structures and dynamics of molecules.

The Force Field method gives the following relationship between the total potential energy and the valence interactions and nonbonded interactions

$$E = E_{\text{R}} + E_{\theta} + E_{\phi} + E_{\omega} + E_{\text{vdw}} + E_{\text{el}}$$

where $E_{\text{R}}$ is the bond stretching of the valence interaction. Valence interactions also include $E_{\theta}$, the bond angle bending; $E_{\Phi}$, the dihedral angle torsion; and $E_{\omega}$ the inversion terms. The nonbonded interactions contain $E_{\text{vdw}}$, the van der Waal's terms, and $E_{\text{el}}$, the electrostatic terms.

The bond stretch interaction can be depicted as the following Morse function.

$$E_{\text{R}} = D_{\text{IJ}}\left [ e^{-\alpha (r-r_{\text{IJ} }) } - 1 \right ]^{2}$$ $$\alpha = \left [ k_{\text{IJ} } / 2D_{\text{IJ} } \right ]^{1/2}$$

where $D_{\text{IJ}}$ is the bond dissociation energy (kcal/mol) and $r_{\text{IJ}}$ is the standard bond length in angstroms. The Morse function is considered to be more accurate in this circumstance since it includes the anharmonic constant $r_{\text{IJ}}$, which nears the equilibrium, and the finite energy $D_{\text{IJ}}$ for breaking bonds.

The general angle bond terms can be described as the following cosine Fourier expansion in $\theta$

$$E_{\theta} = K_{\text{IJK}} \sum_{n = 0}^{m} C_{n} \cos{n\theta }$$

where the coefficient $C_{n}$ aims to satisfy the boundary conditions including that the function has a minimum at the natural bond angle $\theta_{0}$. Herein, since our membrane contains nonlinear structural molecules, we applied a three-term Fourier expansion with three expansion coefficients.

$$E_{\theta} = K_{\text{IJK}} \left [ C_{0} + C_{1}\cos{\theta } + C_{2} \cos{2\theta} \right ]$$ $$C_{2} = 1/(4 \sin^{2}{\theta_{0}}) \quad C_{1} = -4C_{2}\cos{\theta_{0}}$$ $$C_{0} = C_{2}(2\cos^{2}{\theta_{0}} + 1)$$

The angle bend force constants are produced by Badger's rules, which is the second partial derivative of $E$ with respect to $\theta$

$$K_{\text{IJK}} = \frac{\partial^{2} E}{\partial \theta^{2}} = \beta \frac{Z_{I}^{*}Z_{k}^{*}}{r_{\text{IK}}^{5}} r_{\text{IJ}} r_{\text{JK}} \left [ r_{\text{IJ}} r_{\text{JK}} (1 - \cos^{2}{\theta_{0}}) -r_{\text{IK}}^{2} \cos{\theta_{0}} \right ]$$

where $r_{\text{IJ}}$ and $r_{\text{IK}}$ is the natural bond length, $\beta$ is an undetermined coefficient, and $Z_{I}^{*}$ and $Z_{K}^{*}$ are the effective charges of the atoms.

The torsional terms for two bonds $\text{IJ}$ and $\text{KL}$ connected is described with a cosine Fourier expansion in $\phi$

$$E_{\phi} = K_{\text{IJKL}} \sum_{n = 0}^{m} C_{n} \cos{n\phi_{\text{IJKL} }}$$

where $K_{\text{IJKL}}$ and the coefficients $C_{n}$ are determined by the rotational barrier $V_{\phi}$ to explain the periodic trend.

For the inversion energy contribution, a one- or two-term cosine Fourier expansion in $\omega$ is developed, and it's used for atoms $I$, which are bonded to three other atoms $J, K, and L$

$$E_{\omega} = K_{\text{IJKL}}(C_{0} + C_{1}\cos{\omega_{\text{IJKL}}} + C_{2}\cos{2\omega_{\text{IJKL}}})$$

where $K_{IJKL}$ is the force field constant in kcal/mol and $\omega_{\text{IJKL}}$ is the angle between the IL axis and the $IJK$ plane.

The van der Waal energy contribution cannot be neglected, and in order to maintain higher numerical stability, the following exponential-6 form of $E_{\text{vdw}}$ is developed.

$$E_{\text{vdw}} = \left [ D_{\text{IJ}}(\frac{6}{\zeta - 6})e^{\zeta } \right ]e^{-\zeta (x/x_{\text{IJ} })} - \left [ D_{\text{IJ} } (\frac{\zeta }{\zeta - 6} )x_{\text{IJ}}^{6} \right ] / x^{6}$$

where $D_{\text{IJ}}$ is the well depth term, $x_{\text{IJ}}$ is the distance term, and the shape parameter $\zeta$. The the distance term can be expressed as:

$$x_{\text{I}} = \zeta / B_{\text{I}}$$

Herein, the repulsion exponent $B_{\text{I}}$ = 2 $\sqrt{2\text{IP}}_{\text{I}}$ and $\text{IP}_{\text{I}}$ is the ionization energy for each atom $I$. Meanwhile, for the well depth term, we have:

$$D_{\text{I}} = C_{6\text{II}}(\frac{\zeta - 6}{\zeta}) / x_{\text{I}}^{6}$$

where the dispersion terms $C_{6\text{II}}$ is proportional to the upper bound numerical Hartree-Fock values presented by Fraga, Karwowski, and Saxena. Moreover, for the shape parameter $\zeta$ of the second-period atoms, we have:

$$\zeta = 10.02 + 0.6775n$$

for the third-period atoms, we have:

$$\zeta = 8.587+0.897n$$

and for the remaining main group elements, we have:

$$\zeta = 8 + n$$

where $n$ is the atomic number, and the $\zeta$ for the noble gases is assigned to 15.

The electrostatic interactions can be calculated by the following equation:

$$E_{el} = 332.0637(Q_{i}Q_{j}/\epsilon R_{ij})$$

where $Q_{i}$ and $Q_{j}$ are charges in coulomb, $R_{ij}$ is the distance in angstroms, and $\epsilon$ is the dielectric constant.

The procedure of the simulations was carried out in a Newton-Raphson minimization scheme with the norm gradient convergence criteria of $1 \times 10^{-10}$ kcal/mol/$\text{\AA}$, and a force constant matrix with absent negative eigenvalues was used to verify the minima of the simulations.

The model accepts the spatial structure of the filtration site and some contextual information as inputs. We kept track of each MP particle and statistically analyzed the performance of the filtration site, and the resulting video is shown below:

The results of the simulation show the difference in the adsorptivity between our BC-APT membrane and the nanocellulose membrane. Apparently, the adsorption of our BC-APT membranes is much stronger than that of nanocellulose membranes. When a PET molecule comes close to our membrane, it is directly adsorbed and immobilized very thoroughly without any detachment. With nanocellulose membranes, the PET molecule is not adsorbed well and the membrane does not bind the PET molecule securely. As a result, our BC-APT membranes show high adsorptivity and efficiency.

Similarly, this result is validated by the Force Field method, which accurately handles a large number of different intermolecular attractions and molecular collisions and adsorptions in a very short period of time, providing strong simulation evidence for the final results and comparisons of adsorption efficiencies.

Kinetic Model

To determine the adsorption kinetics of the BAM with CBM membrane, we conducted a series of experiments to measure the adsorption rate of the BAM with CBM.

From the experimental data and previous research, we speculate that the reaction path is either pseudo-first-order or pseudo-second-order.

The pseudo-first-order model generally obeys the following linear relationship

$$\frac{\mathrm{d} q_{t}}{\mathrm{d} t} =k_{1}(q_{e}-q_{t})$$

where $q_{e}$ and $q_{t}$ (mg/g) are the amounts of the amount adsorbed at equilibrium and at time $t$ respectively, and $k_{1}$ is the equilibrium rate constant (L/min).

After some integration, the previous equation can be transformed to

$$\ln{(q_{e}-q_{t})} = \ln{(q_{e})}-k_{1}t$$

this equation also implies that whether the reaction is pseudo-first-order can be determined by its half-life. When the reaction has a half-life, then it should be pseudo-first-order.

The pseudo-second-order model can be expressed as

$$\frac{\mathrm{d} q_{t}}{\mathrm{d} t} = k_{2}(q_{e}-q_{t})^{2}$$

In a similar way, by applying integration, we have

$$\frac{t}{q_{t}}=\frac{1}{k_{2}q_{e}^{2}}+\frac{1}{q_{e}}$$

We fit our data set in both pseudo-first and pseudo-second-order models, which are given in the two diagrams above. The adsorption capacity increases rapidly in the first 1.5 minutes and reaches the adsorption equilibrium in about 4 minutes. In addition, the correlation coefficient $R^{2}$ of the pseudo-first-order model is 0.99615, which is higher than the correlation coefficient of the pseudo-second-order model -- 0.98645. Thus, the pseudo-first-order model is more consistent with the adsorptive kinetics of microplastics on our membrane.

Adsorption isotherm

Freundlich adsorption isotherm equation is:

$$q_{e}=k_{F}{c_{e}^{\frac{1}{n}}}$$

where $q_{e}$ is the adsorption amount of the adsorbent in the equilibrium state (mg/g), $c_{e}$ is the equilibrium concentration of the solute in the liquid phase (mg/L), $n$ is the inhomogeneity factor of the adsorption strength of .... ions, which varies with the inhomogeneity of the adsorbent, and $k_{F}$ is the constant related to adsorption capacity.

Langmuir adsorption isotherm equation is:

$$q_{e}=\frac{kq_{max}c_{e}}{1+kc_{e}}$$

where $q_{e}$ is the adsorption amount of the adsorbent in the equilibrium state (mg/g). $c_{e}$ is the equilibrium concentration of the solute in the liquid phase (mg/L), $q_{max}$ is the maximum adsorption capacity of adsorbent on adsorb rate (mg/g), and $k$ is the langmuir constant related to adsorption energy.

Similarly, we fit the data set to both the Freundlich model and Langmuir model, which are shown in the diagrams above, and the parameters are summarized. Apparently, the correlation coefficient $R^{2}$ of the Freundlich model is 0.98883, which is smaller than the $R^{2}$ of the Langmuir model (0.9957). Therefore, the Langmuir model is more suitable than the Freundlich model for depicting the isotherms of microplastics on our membrane.

References

Mohamed Nor, N. H., Kooi, M., Diepens, N. J., & Koelmans, A. A. (2021). Lifetime accumulation of microplastic in children and adults. Environmental Science & Technology, 55(8), 5084–5096.

Kumar, L., Sehrawat, R., & Kong, Y. (2021). Oat proteins: A perspective on functional properties. LWT, 152, 112307. doi:10.1016/j.lwt.2021.112307

Van Duin, A. C. T., Dasgupta, S., Lorant, F., & Goddard, W. A. (2001). ReaxFF: a reactive force field for hydrocarbons. The Journal of Physical Chemistry A, 105(41), 9396–9409.

Chenoweth, K., Van Duin, A. C. T., & Goddard, W. A. (2008). ReaxFF reactive force field for molecular dynamics simulations of hydrocarbon oxidation. The Journal of Physical Chemistry A, 112(5), 1040–1053.

Smirnov, K. S., & van de Graaf, B. (1996). Consistent implementation of the electronegativity equalization method in molecular mechanics and molecular dynamics. J. Chem. Soc. , Faraday Trans., 92, 2469–2474. doi:10.1039/FT9969202469

Desta, I. T., Porter, K. A., Xia, B., Kozakov, D., & Vajda, S. (2020). Performance and its limits in rigid body protein-protein docking. Structure, 28(9), 1071–1081.

Vajda, S., Yueh, C., Beglov, D., Bohnuud, T., Mottarella, S. E., Xia, B., … Kozakov, D. (2017). New additions to the C lus P ro server motivated by CAPRI. Proteins: Structure, Function, and Bioinformatics, 85(3), 435–444.

Kozakov, D., Hall, D. R., Xia, B., Porter, K. A., Padhorny, D., Yueh, C., … Vajda, S. (2017). The ClusPro web server for protein--protein docking. Nature Protocols, 12(2), 255–278.

Kozakov, D., Beglov, D., Bohnuud, T., Mottarella, S. E., Xia, B., Hall, D. R., & Vajda, S. (2013). How good is automated protein docking? Proteins: Structure, Function, and Bioinformatics, 81(12), 2159–2166.

Fraga, S., Karwowski, J., & Saxena, K. M. S. (1 1976). Handbook of atomic data. place = United States, year = 1976, month = 1(volume = , place = United States, year = 1976, month = 1). Retrieved from https://www.osti.gov/biblio/6501446

Jumper, J., Evans, R., Pritzel, A., Green, T., Figurnov, M., Ronneberger, O., … Others. (2021). Highly accurate protein structure prediction with AlphaFold. Nature, 596(7873), 583–589.

Chenoweth, K., van Duin, A. C. T., & Goddard, W. A. (2008). ReaxFF Reactive Force Field for Molecular Dynamics Simulations of Hydrocarbon Oxidation. The Journal of Physical Chemistry A, 112(5), 1040–1053. doi:10.1021/jp709896w

Rappe, A. K., Casewit, C. J., Colwell, K. S., Goddard, W. A., Iii, & Skiff, W. M. (1992). UFF, a full periodic table force field for molecular mechanics and molecular dynamics simulations. Journal of the American Chemical Society, 114(25), 10024–10035. doi:10.1021/ja00051a040

Mayo, S. L., Olafson, B. D., & Goddard, W. A. (1990). DREIDING: a generic force field for molecular simulations. The Journal of Physical Chemistry, 94(26), 8897–8909. doi:10.1021/j100389a010

Casewit, C. J., Colwell, K. S., & Rappe, A. K. (1992a). Application of a universal force field to main group compounds. Journal of the American Chemical Society, 114(25), 10046–10053. doi:10.1021/ja00051a042

Casewit, C. J., Colwell, K. S., & Rappe, A. K. (1992b). Application of a universal force field to organic molecules. Journal of the American Chemical Society, 114(25), 10035–10046. doi:10.1021/ja00051a041