If you are stuck on this page, please try refreshing.



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:

dCMPdt=(1fabsfa,n)i=1n(IriCi)+(1fabs)(fdepInRCa)+ktisCtisklossCMP \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:

dCtisdt=fabsfa,ni=1n(IriCi)+fabs(fdepInRCa)ktisCtis \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 ii refers to the index of the nn food types in the study, CiC_{i} is the microplastics concentration per media in terms of number concentrations, Iri\mathrm{Ir}_{i} represents the related daily ingestion rate, fabsf_{\mathrm{abs}} represents the intestinal absorption, fa,nf_{\mathrm{a, n}} refers to the total number of ingested particles ranged from 1 to 10 μ\mum, CaC_{\mathrm{a}} is the microplastics concentration in the air, InR\mathrm{InR} is the inhalation rate, fdepf_{\mathrm{dep}} is assumed as a constant representing the fraction of inhaled particles that deposited in the nasopharyngeal cavities, CtisC_{\mathrm{tis}} refers to the amount of removal of microplastics from the tissue, ktisk_{\mathrm{tis}} is the biliary excretion rate constant, and klossk_{\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.40Erep0.40Eatt+600Eelec+1.00EDARS 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.


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 Cα\text{C}\alpha Atoms (Cα\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.


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:

Figure 1: thermal stability simulation result

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:

Esystem=Ebond+Eover+Eunder+Eval+Epen+Etors+Econj+EvdWaals+ECoulomb \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 BOij\text{BO}_{ij}' between two atoms can be calculated by the interatomic distance rijr_{ij} directly, as shown in the following equation:

BOij=exp[pbo,1(rijro)pbo,2]+exp[pbo,3(rijπro)pbo,4]+exp[pbo,5(rijππro)pbo,6] \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 pbo,1p_{bo, 1} and pbo,2p_{bo, 2} represent the exponential terms of sigma term with the unity below 1.5A˚\sim 1.5 \text{\AA} and negligible above 2.5A˚\sim 2.5 \text{\AA}, pbo,3p_{bo, 3} and pbo,4p_{bo, 4} mean the exponential terms of the first pi bond with the unity below 1.2A˚\sim 1.2 \text{\AA} and negligible above 1.75A˚\sim 1.75 \text{\AA}, and pbo,5p_{bo, 5} and pbo,6p_{bo, 6} represent the exponential terms of the second pi bond with the unity below 1.0A˚\sim 1.0 \text{\AA} and negligible above 1.4A˚\sim 1.4 \text{\AA}.

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

Ebond=DeBOijexp[pbe, 1(1BOijpbe, 1)] 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 DeD_{e} and pbe,1p_{be, 1} are bond parameters involved in CC\text{C}-\text{C}, CH\text{C}-\text{H}, and HH\text{H}-\text{H} bonds.

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

Eover=poverΔi(11+exp(λ6Δi)) E_{\text{over}} = p_{\text{over}} \cdot \Delta_{i} \cdot \left ( \frac{1}{1 + \text{exp}(\lambda_{6} \cdot \Delta_{i})} \right )

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

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

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

Eunder=punder1exp(μ1Δi)1+exp(μ2Δi)f1(BOij,π,Δj) 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}) f1(BOij,π,Δj)=11+μ3exp(μ4j=1neighbors(i)ΔjBOij,π) 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 EpenE_{\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.

Eval=f2(BOij)f2(BOjk)f3(Δj){kakaexp[kb(ΘoΘijk)2]} 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 \} f2(BOij)=1exp(μ5BOijμ6) f_{2}(\text{BO}_{ij}) = 1 - \text{exp}(-\mu_{5} \cdot \text{BO}_{ij}^{\mu_{6}}) f3(Δj)=2+exp(μ7Δj)1+exp(μ7Δj)+exp(pv,1Δj)[μ8(μ81)2+exp(μ9Δj)1+exp(μ9Δj)+exp(pv,2Δj)] \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*} Θo=πΘ0,0{1exp[μ10(2SBO2)]} \Theta_{o} = \pi - \Theta_{0, 0} \cdot \left \{ 1 - \text{exp}\left [ -\mu _{10} \cdot (2 - \text{SBO}2 ) \right ] \right \} SBO=Δj2{1exp[5(12Δj)μ11]}+n=1neighbors(j)BOjn,π \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 } Δj,2=Δj if Δj<0 \Delta_{j,2} = \Delta_{j} \text{ if } \Delta_{j} < 0 Δj,2=0 if Δj0 \Delta_{j,2} = 0 \text{ if } \Delta_{j} \ge 0 SBO2=0 if SBO0 \text{SBO}2 = 0 \text{ if SBO} \le 0 SBO2=SBOμ12 if 0<SBO<1 \text{SBO}2 = \text{SBO}^{\mu_{12}} \text{ if } 0 < \text{SBO} < 1 SBO2=2(2SBO)μ12 if 1<SBO<2 \text{SBO}2 = 2 - (2 - \text{SBO})^{\mu_{12}} \text{ if } 1 < \text{SBO} < 2 SBO2=2 if SBO>2 \text{SBO}2 = 2 \text{ if SBO} > 2

where Θijk\Theta_{ijk} represents the deviations in valence angle from its equilibrium value Θo\Theta_{o}, the f2(BO)f_{2}(\text{BO}) is used to guarantee the valence angle energy contribution will converge and disappear gradually during bond dissociation, the f3(Δj)f_{3}(\Delta_{j}) copes with the phenomena of over/undercoordination in the central atom jj, the values of Θo\Theta_{o} and Θijk\Theta_{ijk} rely on the sum of π\pi-bond orders (SBO) around the central atom jj on the valence angle energy, and μ5\mu_{5}, μ6\mu_{6}, μ7\mu_{7}, μ8\mu_{8}, μ9\mu_{9}, μ10\mu_{10}, μ11\mu_{11}, μ12\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.

Epen=μ13f4(Δj)exp[μ14(BOij2)2]exp[μ14(BOjk2)2] 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 ] f4(Δj)=2+exp(μ15Δj)1+exp(μ15Δj)+exp(μ16Δj) 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 f4(Δj)f_{4}(\Delta_{j}) copes with over/undercoordination in central atom jj on the penalty energy, and μ13\mu_{13}, μ14\mu_{14}, μ15\mu_{15}, μ16\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 ωijkl\omega_{ijkl} works accurately for bond orders close to 0 and bond orders larger than 1, the following equations were developed.

Etors=f5(BOij,BOjk,BOkl)sinΘijksinΘjkl E_{\text{tors}} = f_{5}(\text{BO}_{ij}, \text{BO}_{jk}, \text{BO}_{kl}) \cdot \sin{\Theta_{ijk}} \cdot \sin{\Theta_{jkl}} \cdot [12V2exp{pl(BOjk3+f6(Δj,Δk))2}(1cos2ωijkl)+12V3(1+cos3ωijkl)] \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 ] f5(BOij,BOjk,BOkl)=[1exp(μ17BOij)] f_{5}(\text{BO}_{ij}, \text{BO}_{jk}, \text{BO}_{kl}) = \left [ 1 - \text{exp}(-\mu_{17} \cdot \text{BO}_{ij} ) \right ] \cdot [1exp(μ17BOjk)][1exp(μ17BOkl)] \left [ 1 - \text{exp}(-\mu _{17} \cdot \text{BO}_{jk} ) \right ] \cdot \left [ 1 - \text{exp}(-\mu _{17} \cdot \text{BO}_{kl} ) \right ] f6(Δj,Δk)=2+exp[μ18(Δj+Δk)]1+exp[μ18(Δj+Δk)]+exp[μ19(Δj+Δk)] 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 V2V_{2} and V3V_{3} are two sets of terms that are based on the bond order of the central bond BOjk\text{BO}_{jk}, the value of sinΘijksinΘjkl\sin{\Theta_{ijk}} \cdot \sin{\Theta_{jkl}} guarantees that the torsion energy contribution will vanish when either Θijk\Theta_{ijk} or Θjkl\Theta_{jkl} reaches π\pi, the f5(BOij,BOjk,BOkl)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 f6(Δj,Δk)f_{6}(\Delta_{j}, \Delta_{k}) prohibits the occurrence of excessive torsion energy contributions when atom jj and kk are overcoordinated, and μ17,μ18,μ19\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.

Econj=f7(BOij,BOjk,BOkl)μ20[1+(cos2ωijkl1)sinΘijksinΘjkl] 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 ] f7(BOij,BOjk,BOkl)=exp[μ21(BOij32)2] 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 exp[μ21(BOij32)2]exp[μ21(BOkl32)2] \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 μ20,μ21\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

EvdWaals=Dij{exp[αij(1f8(rij)rvdW)]2exp[12αij(1f8(rij)rvdW)]} 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 \} f8(rij)=[rijμ22+(1μw)μ22]μ22 f_{8}(r_{ij}) = \left [ r_{ij}^{\mu_{22}} + (\frac{1}{\mu_w{}} )^{\mu_{22}} \right ]^{\mu_{22}}

where DijD_{ij} is the bond parameter, the f8(rij)f_{8}(r_{ij}) provides a shielded interaction to avoid repulsion and atoms sharing valence angle, and μ22\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

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

where γij\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.

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

where Δt\Delta t is the time step, τB\tau_{B} is the time constant in the Berendsen thermostat, T0T_{0} is the desired temperature, and TT 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=ER+Eθ+Eϕ+Eω+Evdw+Eel E = E_{\text{R}} + E_{\theta} + E_{\phi} + E_{\omega} + E_{\text{vdw}} + E_{\text{el}}

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

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

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

where DIJD_{\text{IJ}} is the bond dissociation energy (kcal/mol) and rIJr_{\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 rIJr_{\text{IJ}}, which nears the equilibrium, and the finite energy DIJD_{\text{IJ}} for breaking bonds.

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

Eθ=KIJKn=0mCncosnθ E_{\theta} = K_{\text{IJK}} \sum_{n = 0}^{m} C_{n} \cos{n\theta }

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

Eθ=KIJK[C0+C1cosθ+C2cos2θ] E_{\theta} = K_{\text{IJK}} \left [ C_{0} + C_{1}\cos{\theta } + C_{2} \cos{2\theta} \right ] C2=1/(4sin2θ0)C1=4C2cosθ0 C_{2} = 1/(4 \sin^{2}{\theta_{0}}) \quad C_{1} = -4C_{2}\cos{\theta_{0}} C0=C2(2cos2θ0+1) 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 EE with respect to θ\theta

KIJK=2Eθ2=βZIZkrIK5rIJrJK[rIJrJK(1cos2θ0)rIK2cosθ0] 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 rIJr_{\text{IJ}} and rIKr_{\text{IK}} is the natural bond length, β\beta is an undetermined coefficient, and ZIZ_{I}^{*} and ZKZ_{K}^{*} are the effective charges of the atoms.

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

Eϕ=KIJKLn=0mCncosnϕIJKL E_{\phi} = K_{\text{IJKL}} \sum_{n = 0}^{m} C_{n} \cos{n\phi_{\text{IJKL} }}

where KIJKLK_{\text{IJKL}} and the coefficients CnC_{n} are determined by the rotational barrier Vϕ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 II, which are bonded to three other atoms J,K,andLJ, K, and L

Eω=KIJKL(C0+C1cosωIJKL+C2cos2ωIJKL) E_{\omega} = K_{\text{IJKL}}(C_{0} + C_{1}\cos{\omega_{\text{IJKL}}} + C_{2}\cos{2\omega_{\text{IJKL}}})

where KIJKLK_{IJKL} is the force field constant in kcal/mol and ωIJKL\omega_{\text{IJKL}} is the angle between the IL axis and the IJKIJK plane.

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

Evdw=[DIJ(6ζ6)eζ]eζ(x/xIJ)[DIJ(ζζ6)xIJ6]/x6 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 DIJD_{\text{IJ}} is the well depth term, xIJx_{\text{IJ}} is the distance term, and the shape parameter ζ\zeta. The the distance term can be expressed as:

xI=ζ/BI x_{\text{I}} = \zeta / B_{\text{I}}

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

DI=C6II(ζ6ζ)/xI6 D_{\text{I}} = C_{6\text{II}}(\frac{\zeta - 6}{\zeta}) / x_{\text{I}}^{6}

where the dispersion terms C6IIC_{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:

ζ=10.02+0.6775n \zeta = 10.02 + 0.6775n

for the third-period atoms, we have:

ζ=8.587+0.897n \zeta = 8.587+0.897n

and for the remaining main group elements, we have:

ζ=8+n \zeta = 8 + n

where nn 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:

Eel=332.0637(QiQj/ϵRij) E_{el} = 332.0637(Q_{i}Q_{j}/\epsilon R_{ij})

where QiQ_{i} and QjQ_{j} are charges in coulomb, RijR_{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×10101 \times 10^{-10} kcal/mol/A˚\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:

Simulation of the absorption of the PET by the nano-fibres
Simulation of the absorption of the PET by BC-APT membrane

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

dqtdt=k1(qeqt) \frac{\mathrm{d} q_{t}}{\mathrm{d} t} =k_{1}(q_{e}-q_{t})

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

After some integration, the previous equation can be transformed to

ln(qeqt)=ln(qe)k1t \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

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

In a similar way, by applying integration, we have

tqt=1k2qe2+1qe \frac{t}{q_{t}}=\frac{1}{k_{2}q_{e}^{2}}+\frac{1}{q_{e}}
Pseudo-first-order model Pseudo-second-order model

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 R2R^{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:


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

Langmuir adsorption isotherm equation is:

qe=kqmaxce1+kce q_{e}=\frac{kq_{max}c_{e}}{1+kc_{e}}

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

Freundlich model Langmuir model

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 R2R^{2} of the Freundlich model is 0.98883, which is smaller than the R2R^{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.


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

On this page:



Copyright @ 2023 BAID-CHINA