Automatic optimization method of force field parameters of multi-system combined constrained bio-membrane mesoscopic model
By employing a multi-system joint constraint method and utilizing parallel computing and particle swarm optimization algorithms, the force field parameters of a biomembrane mesoscopic model are constructed directly using macroscopic experimental observations as the target. This solves the problems of parameter space nonlinearity and reliance on all-atom simulation data in traditional methods, and achieves efficient and robust parameter optimization and accurate simulation of multibody systems.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- ZHEJIANG UNIV
- Filing Date
- 2026-05-25
- Publication Date
- 2026-06-19
AI Technical Summary
Existing technologies lack rigorous physical derivations when constructing mesoscopic force fields, resulting in highly nonlinear parameter spaces. Traditional methods struggle to avoid unstable evolution and rely on expensive all-atom simulation data, limiting the application of mesoscopic simulations in engineering and biomedical fields.
A multi-system joint constraint approach is adopted, which uses a parallel co-simulation system, particle swarm optimization algorithm (PSO) and LAMMPS simulation engine to directly construct the force field parameters of the multibody system based on macroscopic experimental observations. Soft constraint penalty terms are introduced to avoid non-physical states, and MPI parallel computing is used to accelerate the optimization process.
It achieves efficient and automated force field parameter optimization, with parameters that are portable across different physical systems, reducing computational costs, improving the robustness of simulation results and R&D efficiency, and enabling accurate reproduction of key experimental indicators of multibody systems.
Smart Images
Figure CN122242311A_ABST
Abstract
Description
Technical Field
[0001] This invention belongs to the field of computational physics and mesoscopic simulation technology, specifically relating to a method for automatically optimizing and inversely designing force field parameters for a multibody dissipative particle dynamics model by using machine learning algorithms combined with physical constraint mechanisms. Background Technology
[0002] Studying emergent phenomena in soft matter systems such as biomembranes requires traversing spatiotemporal scales that are difficult to achieve with atomic simulations. Mesoscopic simulation techniques, such as dissipative particle dynamics (DPD) and its many-body variant (MDPD), have become important tools bridging the microscopic and macroscopic descriptions due to their inherent inclusion of fluid dynamics and thermal fluctuations. However, constructing mesoscopic force fields faces significant challenges. The interaction parameters in MDPDs often lack rigorous physical derivations, and traditional methods rely heavily on empirical intuition or bottom-up mappings based on atomic trajectories, leading to a disconnect between the mesoscopic force field and experimentally measurable macroscopic properties (such as interfacial tension and film thickness). MDPD introduces density-dependent many-body interactions, resulting in a highly nonlinear and high-dimensional parameter space. Small parameter changes can lead to drastic abrupt changes in phase behavior or interfacial stability, rendering traditional manual trial-and-error or local optimization strategies ineffective. Existing methods, during the parameter space search, inevitably enter parameter regions that lead to unstable evolution of the simulated system, severely impacting the usability of simulation results and limiting the application of mesoscopic simulations in engineering and biomedicine. While machine learning-based methods (such as neural networks) are effective, they heavily rely on large amounts of microscopic atomic simulation data as "true values," resulting in high computational costs and difficulty in handling long relaxation processes that atomic simulations cannot sample. Therefore, there is an urgent need for a method that can automatically, efficiently, and physically consistently back-calculate MDPD force field parameters using macroscopic experimental observations (such as surface tension, bending stiffness, and monolipid occupancy area) as targets. Summary of the Invention
[0003] The purpose of this invention is to overcome the shortcomings of existing technologies and provide an automatic optimization method for force field parameters of a biomembrane mesoscopic model with multi-system joint constraints. This method constructs the mesoscopic force field without relying on all-atom simulations, avoids non-physical solutions in a high-dimensional nonlinear parameter space, and ensures the portability of the obtained parameters across different physical systems.
[0004] The specific technical solution adopted in this invention is as follows:
[0005] In a first aspect, the present invention provides an automatic optimization method for force field parameters of a biomembrane mesoscopic model under multi-system joint constraints, as detailed below:
[0006] S1: Construct a parallel co-simulation system that treats multiple physical systems as a joint evaluation unit;
[0007] S2: Implement parameter decoupling and hierarchical mapping strategies for the parallel co-simulation system described in S1, and discretize the global force field parameter set according to the interaction type between particles and the corresponding physical response characteristics.
[0008] S3: Preset parameter boundaries of each space after discretization in S2 An internal randomized population of N particles is initialized, where each particle represents a set of potential force field parameters. The position vectors of the particles are initialized. and velocity vector ;
[0009] S4: For each set of force field parameters in the population obtained in S3 or S7, the LAMMPS simulation engine is started simultaneously using the MPI parallel mechanism to calculate the parallel co-simulation system described in S1.
[0010] S5: Post-process the simulated trajectory obtained in S4 to extract macroscopic physical properties and calculate physical constraint indicators;
[0011] S6: Based on the deviation between the simulation results obtained in S5 and the preset experimental reference values, construct the loss function. ;
[0012] S7: Update the individual historical best position of each particle based on the loss function value calculated according to the loss function described in S6. The global optimal position of the entire population According to the PSO speed update formula ,according to Update the next generation parameters; among them, Let K be the velocity vector of particle k at the (t+1)th iteration. Let K be the velocity vector of particle k at the t-th iteration. Let be the inertia weight for the t-th iteration. For individual learning factors, As a social learning factor; A random number between [0, 1], used to increase the randomness of individual searches; A random number between [0, 1], used to increase the randomness of the group search; Let be the current position vector of particle k at the t-th iteration, i.e., the current parameter solution;
[0013] S8: Repeat S4~S7 until the maximum number of iterations is reached or the loss function converges to the preset threshold, and finally output the obtained optimal force field parameter set.
[0014] Preferably, in S1, the parallel co-simulation system includes a water-air interface system for calibrating the fluid state equation and the surface tension of the water phase, an oil-water interface system for anchoring the hydrophobic interaction scale and the oil-water interface tension, and a phospholipid bilayer system for precisely controlling the structural and mechanical properties of the biofilm.
[0015] Furthermore, the water-air interface system is a vacuum-water-vacuum sandwich structure, the oil-water interface system is a water-oil-water sandwich structure, and the phospholipid bilayer system is a water-phospholipid bilayer-water sandwich structure.
[0016] Preferably, in S2, the global force field parameter set is discretized into: a solvent interaction subspace characterizing the intrinsic properties of the fluid, an interface interaction subspace characterizing the properties of multiphase contact, and a membrane force field subspace characterizing the self-assembly properties of phospholipid molecules.
[0017] Preferably, S5 includes:
[0018] Calculation of water-air interfacial tension using pressure tensor and oil-water interfacial tension ;
[0019] Calculate film thickness using density distribution ;
[0020] By simulating box size Calculate the area occupied by monolipids ;in, It is the number of phospholipids in the system;
[0021] By fitting the Helfrich-Canham theoretical spectrum to the discrete Fourier transform of the bilayer height field, Calculate the bending stiffness of the membrane ;in, It is Boltzmann's constant. It's temperature. It is tension. It is a two-dimensional wave vector The model;
[0022] The number of water molecules that permeate into the hydrophobic core region of the membrane is counted to determine whether the membrane has ruptured;
[0023] Calculate the diffusion coefficient of the lipid tail chain to determine whether the membrane is in a fluid state;
[0024] The orientation order of the lipid tail chains is calculated to determine whether they conform to the characteristics of a liquid crystal phase.
[0025] Preferably, in S6, the loss function The calculation formula is as follows:
[0026] ;
[0027] in, , , ; This is a soft constraint penalty function used to eliminate abnormal states caused by non-physical parameters; the subscripts sim and exp represent the simulated value and the experimental reference value, respectively. For the weighting coefficients of the loss term in the water-air interface system, For the weighting coefficients of the loss term in the water-oil interface system, For the weighting coefficients of the loss term in the phospholipid bilayer membrane interface system, The weighting coefficient for the area loss term of single lipids. The weighting coefficient for the lipid membrane thickness loss term. is the weighting coefficient for the lipid membrane bending modulus loss term, and j is the index of the soft constraint penalty function.
[0028] Secondly, the present invention provides an automatic optimization system for force field parameters of a biomembrane mesoscopic model with multi-system joint constraints, comprising:
[0029] The parallel co-simulation system construction module is used to build a parallel co-simulation system that treats multiple physical systems as a joint evaluation unit.
[0030] The discretization module is used to implement parameter decoupling and hierarchical mapping strategies for the parallel co-simulation system, and to discretize the global force field parameter set according to the interaction type between particles and the corresponding physical response characteristics.
[0031] The particle swarm optimization algorithm module is initialized to define the preset parameter boundaries of each space after discretization in the discretization processing module. An internal randomized population of N particles is initialized, where each particle represents a set of potential force field parameters. The position vectors of the particles are initialized. and velocity vector ;
[0032] The simulated trajectory calculation module is used to simultaneously start the LAMMPS simulation engine using the MPI parallel mechanism to calculate each set of force field parameters in the population obtained by the initialization particle swarm optimization algorithm module or update iteration module, and to perform calculations on the parallel co-simulation system.
[0033] The post-processing module is used to post-process the simulated trajectory obtained by the simulated trajectory calculation module, extract macroscopic physical properties and calculate physical constraint indicators;
[0034] The loss function construction module is used to construct a loss function based on the deviation between the simulation results obtained by the post-processing module and the preset experimental reference values. ;
[0035] The update iteration module is used to update the individual historical best position of each particle based on the loss function value calculated by the loss function. The global optimal position of the entire population According to the PSO speed update formula ,according to Update the next generation parameters; among them, Let K be the velocity vector of particle k at the (t+1)th iteration. Let K be the velocity vector of particle k at the t-th iteration. Let be the inertia weight for the t-th iteration. For individual learning factors, As a social learning factor; A random number between [0, 1], used to increase the randomness of individual searches; A random number between [0, 1], used to increase the randomness of the group search; Let be the current position vector of particle k at the t-th iteration, i.e., the current parameter solution;
[0036] The results output module is used to repeatedly run the simulation trajectory calculation module, post-processing module, loss function construction module, and update iteration module in sequence until the maximum number of iterations is reached or the loss function converges to a preset threshold, and finally outputs the obtained optimal force field parameter set.
[0037] Thirdly, the present invention provides a computer program product, including a computer program / instruction, which, when executed by a processor, enables the automatic optimization method for force field parameters of a biomembrane mesoscopic model with multi-system joint constraints as described in any of the first aspects.
[0038] Fourthly, the present invention provides a computer-readable storage medium storing a computer program that, when executed by a processor, implements an automatic optimization method for force field parameters of a biomembrane mesoscopic model with multi-system joint constraints as described in any of the first aspects.
[0039] Fifthly, the present invention provides a computer electronic device, including a memory and a processor;
[0040] The memory is used to store computer programs;
[0041] The processor is configured to, when executing the computer program, implement the automatic optimization method for force field parameters of a biomembrane mesoscopic model with multi-system joint constraints as described in any of the first aspects.
[0042] Compared with the prior art, the present invention has the following advantages:
[0043] 1) Traditional manual parameter tuning struggles to simultaneously account for multiple macroscopic properties. This invention directly uses experimentally measured macroscopic physical quantities as optimization targets, enabling precise reproduction of key experimental indicators for multibody systems. Example data shows that the POPC bilayer force field parameters obtained using this method, which simulate the monolipid occupancy area and membrane thickness, are in high agreement with experimental values and outperform some existing coarse-grained force fields.
[0044] 2) Traditional MDPD force field parameterization typically relies on researchers' experience and iterative trial and error, often taking months. This invention utilizes parallel computing and particle swarm optimization algorithms to achieve fully automated parameter optimization, shortening the parameterization cycle from months to days, greatly improving R&D efficiency. Unlike machine learning methods that rely on massive amounts of all-atom molecular dynamics trajectory data, this invention does not require expensive underlying ground truth data; only a small amount of macroscopic experimental values are needed to complete the inversion, significantly reducing the requirements for computing resources and data reserves.
[0045] 3) Existing technologies often tune parameters for a single system, causing these parameters to fail in other environments. This invention, by simultaneously optimizing three systems—the water-air interface, the oil-water interface, and the phospholipid bilayer—ensures that the obtained parameters are physically self-consistent under different hydrophilic and hydrophobic environments. Data shows that this strategy not only reproduces the mechanical properties of the membrane but also precisely anchors the surface tension of water and the interfacial tension of the oil-water interface, ensuring the robustness of the force field.
[0046] 4) The parameters obtained by the method of this invention can be directly used in dissipative particle dynamics simulation software for simulation calculation of biomembrane self-assembly, mechanical response and transport behavior. Attached Figure Description
[0047] Figure 1 A flowchart for parameter optimization simulation based on PSO algorithm and MDPD modeling; Figure 2 Detailed workflow for co-simulation driven by a parallel computing pool.
[0048] Figure 3 This is a schematic diagram of a water-air interface, an oil-water interface, and a POPC bilayer membrane system.
[0049] Figure 4 MDPD model for coarse-grained water, hexadecane, and POPC molecules.
[0050] Figure 5 The convergence of the global loss function (dark green square) and its contributions to the decomposition of the water-air interface, oil-water interface and lipid bilayer system was achieved in 50 iterations.
[0051] Figure 6 The optimized trajectory for water-air surface tension is shown. Simulated values (cyan squares) are compared with experimental data (orange dashed lines) and the standard Martini model (purple dashed lines) as benchmarks.
[0052] Figure 7 The optimized trajectory for surface tension in an oil-water system is shown. Simulated values (cyan squares) are compared with experimental data (orange dashed lines) and the standard Martini model (purple dashed lines) as benchmarks.
[0053] Figure 8 Simultaneous optimization of the properties of the POPC bilayer structure. Curves trace the evolution of APL (left axis) and film thickness (right axis) toward their respective experimental targets for optimal fit.
[0054] Figure 9 This is a time series of lipid bilayer self-assembly in random solution. The system successfully converged to a stable equilibrium bilayer structure, validating the thermodynamic parameters. Detailed Implementation
[0055] To make the above-mentioned objects, features, and advantages of the present invention more apparent and understandable, specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings. Many specific details are set forth in the following description to provide a thorough understanding of the present invention. However, the present invention can be practiced in many other ways different from those described herein, and those skilled in the art can make similar modifications without departing from the spirit of the present invention. Therefore, the present invention is not limited to the specific embodiments disclosed below. Technical features in various embodiments of the present invention can be combined accordingly without mutual conflict.
[0056] This invention provides an automatic optimization method for force field parameters of a biomembrane mesoscopic model under multi-system joint constraints. This method transforms the construction of microscopic force field parameters into a physically constrained global optimization problem. Through multi-system joint simulation and iterative particle swarm optimization (PSO) algorithm, it automatically searches for the optimal parameter set that can simultaneously reproduce macroscopic experimental properties (such as surface tension, membrane thickness, bending stiffness, etc.) and satisfy thermodynamic stability. Figure 1 As shown, the framework maps physical problems involving the determination of unknown multidimensional parameters into an iterative search loop; the PSO algorithm performs a global search to generate key candidate parameters, which are then rigorously evaluated using LAMMPS simulations and actual performance data to obtain the optimal parameter set.
[0057] Specifically, the method of the present invention has the following key improvements:
[0058] 1) Existing automatic parameterization methods typically aim only at minimizing numerical errors in macroscopic properties. In the high-dimensional nonlinear parameter space of many-body dissipative particle dynamics, there exist numerous "pseudo-solutions"—parameter sets that mathematically satisfy the error requirements but physically exhibit membrane rupture, over-crystallization, or unsteady structures. The method of this invention introduces a soft physics constraint penalty term into the optimization objective function. The forced search algorithm avoids thermodynamically unstable parameter regions, ensuring that the generated force field parameters not only fit the experimental data numerically, but are also physically correct and stable, which greatly improves the robustness of the model.
[0059] 2) Existing mainstream coarse-grained methods (such as Boltzmann inversion, force matching, and relative entropy minimization) heavily rely on expensive all-atom molecular dynamics trajectories as training data and struggle to handle long relaxation processes that cannot be sampled in all-atom simulations. The method of this invention uses a particle swarm optimization algorithm to directly connect with macroscopic experimental observations, eliminating the need for all-atom simulation data as prior knowledge. Utilizing MPI parallel acceleration technology, the manual parameter tuning process, which originally required months, is shortened to a few days of automated computation. This significantly reduces the barrier to entry and computational cost for force field development, enabling the handling of force field construction problems in large-scale, long-term self-assembly systems that are difficult to address in atomic simulations.
[0060] 3) Existing methods typically fit parameters to a single, specific biofilm system, leading to overfitting. For example, to fit the membrane area, solvent properties may be incorrectly sacrificed, resulting in a force field that fails to accurately describe the oil-water interface or the behavior of the solvent alone, lacking versatility. The method of this invention establishes a co-simulation framework encompassing three baseline systems: the water-air interface, the oil-water interface, and the phospholipid bilayer. In each parameter evaluation, the properties of these three systems are simultaneously calculated and optimized. This ensures that the obtained force field parameters have good physical consistency and portability, applicable to complex simulations of different hydrophilic and hydrophobic environments.
[0061] The steps of the method of the present invention will be described in detail below.
[0062] S1: To ensure the physical consistency and portability of the force field parameters, construct as follows: Figure 2 The parallel co-simulation system shown treats multiple physical systems as a single joint evaluation unit. For example... Figure 2 As shown, multiple physical systems (systems A, B, and C) are simultaneously simulated using the LAMMPS engine to evaluate candidate systems. This process includes template configuration, trajectory generation, and post-processing to calculate the loss. The aggregated loss metrics are then fed into the main PSO process to update the search strategy and distribute a new set of parameters for the next iteration.
[0063] As a preferred embodiment of the present invention, the parallel co-simulation system mainly includes three reference physical systems for parallel computing: a water-air interface system for calibrating the fluid state equation and the surface tension of the water phase; an oil-water interface system for anchoring the hydrophobic interaction scale and the oil-water interface tension; and a phospholipid bilayer system for precisely controlling the structural properties (such as area and thickness) and mechanical properties (such as bending stiffness) of the biofilm.
[0064] The three systems mentioned above are simulated as a whole during the optimization process to ensure that the obtained parameters are applicable under different thermodynamic environments.
[0065] Among them, the water-air interface system is a vacuum-water-vacuum "sandwich" structure, with water represented by discrete particles with a number density of about 7; the oil-water interface system is a water-oil-water "sandwich" structure, with hexadecane as an example of oil, and each molecule represented by 4 beads connected by harmonic bonds and harmonic bond angles; taking POPC membrane as an example, the phospholipid bilayer system is a water-phospholipid bilayer-water "sandwich" structure, with phospholipid molecules represented by 12 beads connected by harmonic bonds and harmonic bond angles, including hydrophilic head particles and hydrophobic tail particles the same as oil particles.
[0066] In practical applications, this step can be modified to fit a single or multiple physical systems when constructing them, and can be transferred to other coarse-grained molecular dynamics and dissipative force quantum dynamics variant models.
[0067] S2: Implement parameter decoupling and hierarchical mapping strategies for the parallel co-simulation system obtained in S1. Discretize the global force field parameter set according to the interaction type between particles and the corresponding physical response characteristics.
[0068] As a preferred embodiment of the present invention, the global force field parameter set in this step should be discretized into: a solvent interaction subspace characterizing the intrinsic properties of the fluid, an interface interaction subspace characterizing the multiphase contact properties, and a membrane force field subspace characterizing the self-assembly properties of phospholipid molecules.
[0069] S3, Initialize the Particle Swarm Optimization (PSO) algorithm: This involves defining the predefined parameter boundaries for each space after discretization in S2. An internal randomized population of N particles is initialized, where each particle represents a set of potential force field parameters. The position vectors of the particles are initialized. and velocity vector .
[0070] In practical applications, this step can also be implemented using swarm intelligence algorithms, regression models, or other search algorithms.
[0071] S4: For each set of force field parameters in the population obtained from S3 or S7, the LAMMPS simulation engine is started simultaneously using the MPI parallel mechanism to calculate the parallel co-simulation system obtained from S1.
[0072] S5: Post-process the simulated trajectory obtained in S4 to extract macroscopic physical properties and calculate physical constraint indicators.
[0073] In a preferred embodiment of the present invention, the steps are as follows:
[0074] Calculation of water-air interfacial tension using pressure tensor and oil-water interfacial tension ;
[0075] Calculate film thickness using density distribution ;
[0076] By simulating box size Calculate the area occupied by monolipids ;in, It is the number of phospholipids in the system;
[0077] By fitting the Helfrich-Canham theoretical spectrum to the discrete Fourier transform of the bilayer height field, Calculate the bending stiffness of the membrane ;in, It is Boltzmann's constant. It's temperature. It is tension. It is a two-dimensional wave vector The model;
[0078] The number of water molecules that permeate into the hydrophobic core region of the membrane is counted to determine whether the membrane has ruptured;
[0079] Calculate the diffusion coefficient of the lipid tail chain to determine whether the membrane is in a fluid state;
[0080] The orientation order of the lipid tail chains is calculated to determine whether they conform to the characteristics of a liquid crystal phase.
[0081] S6: Based on the deviation between the simulation results obtained in S5 and the preset experimental reference values, construct the loss function. .
[0082] As a preferred embodiment of the present invention, the loss function It consists of a weighted error term based on macroscopic attributes and a physical constraint penalty term, and its calculation formula is as follows:
[0083] ;
[0084] in, , , ; This is a soft constraint penalty function used to eliminate abnormal states caused by non-physical parameters; the subscripts sim and exp represent the simulated value and the experimental reference value, respectively. For the weighting coefficients of the loss term in the water-air interface system, For the weighting coefficients of the loss term in the water-oil interface system, For the weighting coefficients of the loss term in the phospholipid bilayer membrane interface system, The weighting coefficient for the area loss term of single lipids. The weighting coefficient for the lipid membrane thickness loss term. is the weighting coefficient for the lipid membrane bending modulus loss term, and j is the index of the soft constraint penalty function.
[0085] In practical applications, the objective function constructed in this step can be modified to other observable target physical quantities as needed, such as diffusion coefficient, membrane compressibility coefficient, etc.
[0086] S7: Update the individual historical best position of each particle based on the loss function value calculated from the loss function obtained in S6. The global optimal position of the entire population According to the PSO speed update formula ,according to Update the next generation parameters; among them, Let K be the velocity vector of particle k at the (t+1)th iteration. Let K be the velocity vector of particle k at the t-th iteration. Let be the inertia weight for the t-th iteration. For individual learning factors, As a social learning factor; A random number between [0, 1], used to increase the randomness of individual searches; A random number between [0, 1], used to increase the randomness of the group search; Let be the current position vector of particle k in the t-th iteration, i.e., the current parameter solution.
[0087] In practical applications, this step can also be implemented using swarm intelligence algorithms, regression models, or other search algorithms.
[0088] S8: Repeat S4~S7 until the maximum number of iterations is reached or the loss function converges to the preset threshold, and finally output the obtained optimal force field parameter set.
[0089] The methods and effects of the present invention will be specifically illustrated below through examples.
[0090] Example
[0091] This embodiment provides an automatic optimization method for force field parameters of a biomembrane mesoscopic model with multi-system joint constraints. The method is as follows:
[0092] S1: Construct a parallel co-simulation system that treats multiple physical systems as a single joint evaluation unit.
[0093] In this embodiment, the parallel co-simulation system includes three reference physical systems for parallel computing: a water-air interface system for calibrating the fluid state equation and the surface tension of the water phase; an oil-water interface system for anchoring the hydrophobic interaction scale and the oil-water interface tension; and a phospholipid bilayer system for precisely controlling the structural properties (such as area and thickness) and mechanical properties (such as bending stiffness) of the biofilm.
[0094] The three systems described above are simulated simultaneously as a whole during the optimization process to ensure that the obtained parameters are applicable under different thermodynamic conditions. The three systems are as follows: Figure 3 and 4 As shown.
[0095] Among them, the water-air interface system is a vacuum-water-vacuum "sandwich" structure, with water represented by discrete particles with a number density of about 7; the oil-water interface system is a water-oil-water "sandwich" structure, with hexadecane as an example of oil, and each molecule represented by 4 beads connected by harmonic bonds and harmonic bond angles; taking POPC membrane as an example, the phospholipid bilayer system is a water-phospholipid bilayer-water "sandwich" structure, with phospholipid molecules represented by 12 beads connected by harmonic bonds and harmonic bond angles, including hydrophilic head particles and hydrophobic tail particles the same as oil particles.
[0096] S2: Implement parameter decoupling and hierarchical mapping strategies for the parallel co-simulation system obtained in S1. Discretize the global force field parameter set according to the interaction type between particles and the corresponding physical response characteristics.
[0097] In this embodiment, the global force field parameter set should be discretized into: a solvent interaction subspace characterizing the intrinsic properties of the fluid, an interface interaction subspace characterizing the properties of multiphase contact, and an intramembrane force field subspace characterizing the self-assembly properties of phospholipid molecules.
[0098] S3, Initialize the Particle Swarm Optimization (PSO) algorithm: This involves defining the predefined parameter boundaries for each space after discretization in S2. An internal randomized population of N particles is initialized, where each particle represents a set of potential force field parameters. The position vectors of the particles are initialized. and velocity vector The preset parameters and parameter boundaries are shown in Table 1.
[0099] Table 1. Search space boundaries for optimizing interactions and intramolecular parameters.
[0100]
[0101] S4: For each set of force field parameters obtained in the population from S3 or S7, the LAMMPS simulation engine is simultaneously launched using the MPI parallel mechanism to calculate the parallel co-simulation system obtained in S1. For example... Figure 3 For water-gas and oil-water interface systems, each candidate parameter set first undergoes a 10,000-step equilibration phase, followed by a 40,000-step production run. For the complex POPC phospholipid bilayer system, a 30,000-step equilibration phase is performed, followed by an 80,000-step production run.
[0102] S5: Post-process the simulated trajectory obtained in S4 to extract macroscopic physical properties and calculate physical constraint indicators. In this embodiment, this step is as follows:
[0103] Calculation of water-air interfacial tension using pressure tensor and oil-water interfacial tension ;
[0104] Calculate film thickness using density distribution ;
[0105] By simulating box size Calculate the area occupied by monolipids ;in, It is the number of phospholipids in the system;
[0106] By fitting the Helfrich-Canham theoretical spectrum to the discrete Fourier transform of the bilayer height field, Calculate the bending stiffness of the membrane ;in, It is Boltzmann's constant. It's temperature. It is tension. It is a two-dimensional wave vector. yes The model;
[0107] The number of water molecules that permeate into the hydrophobic core region of the membrane is counted to determine whether the membrane has ruptured;
[0108] Calculate the diffusion coefficient of the lipid tail chain to determine whether the membrane is in a fluid state;
[0109] The orientation order of the lipid tail chains is calculated to determine whether they conform to the characteristics of a liquid crystal phase.
[0110] S6: Based on the deviation between the simulation results obtained in S5 and the preset experimental reference values, construct the loss function. In this embodiment, the loss function It consists of a weighted error term based on macroscopic attributes and a physical constraint penalty term, and its calculation formula is as follows:
[0111] ;
[0112] in, , , ; This is a soft constraint penalty function used to eliminate abnormal states caused by non-physical parameters; the subscripts sim and exp represent the simulated value and the experimental reference value, respectively. For the weighting coefficients of the loss term in the water-air interface system, For the weighting coefficients of the loss term in the water-oil interface system, For the weighting coefficients of the loss term in the phospholipid bilayer membrane interface system, The weighting coefficient for the area loss term of single lipids. The weighting coefficient for the lipid membrane thickness loss term. is the weighting coefficient for the lipid membrane bending modulus loss term, and j is the index of the soft constraint penalty function. , , Within the phospholipid system, the weights are respectively , , The number of water molecules penetrating the hydrophobic core region of the membrane. It must be less than the tolerance threshold of 5. The mean square displacement of the lipid tail chain must be greater than... Orientational order parameter of lipid tail chains The value is between 0.35 and 0.75. If any of the above constraints are violated, the value of the soft constraint penalty function is set to 100.
[0113] S7: Based on the loss function value calculated from the loss function obtained in S6, update the individual historical best position 𝐩best, 𝑘 of each particle, and the global best position 𝐠best of the entire population, according to the PSO velocity update formula. ,according to Update the next generation parameters; among them, Let K be the velocity vector of particle k at the (t+1)th iteration. Let K be the velocity vector of particle k at the t-th iteration. Let be the inertia weight for the t-th iteration. For individual learning factors, As a social learning factor, A random number between [0, 1] used to increase the randomness of individual searches. Let be a random number between [0, 1] used to increase the randomness of the group search. Let be the current position vector of particle k at the t-th iteration (i.e., the current parameter solution). Individual learning factor. Social learning factor Inertial weight A linear decay strategy is adopted, and during the optimization process, from linearly reduced to To ensure stability, the speed limit is set at 20% of the parameter range.
[0114] S8: Repeat S4~S7 until the maximum number of iterations is reached or the loss function converges to a preset threshold (e.g., ...). Figure 5 As shown in Table 2, the optimal force field parameter set is finally output. The particle swarm optimization algorithm has a population size of N = 96 and performs 50 iterations.
[0115] Table 2. Optimized interactions and intramolecular parameters.
[0116]
[0117] Figures 3-8 The parameterization process of MDPD biofilm based on the PSO algorithm is shown (steps S1-S8). As can be seen from the figure, the loss function approaches 0 with increasing iterations. The optimized interfacial tension of the water-air interface, the interfacial tension of the water-oil interface, the APL of the biofilm, and its thickness are approximately 65.6 mN / m, 54.3 mN / m, and 0.586 nm, respectively. 2 3.663 nm, which is close to the experimental value (e.g. Figures 5-8 (As shown).
[0118] Figure 9 The self-assembly simulation of the phospholipid system further demonstrated that the biomembrane parameters obtained in step S8 could accurately capture the bilayer structure of the lipid-water system using a random mixing system. As shown in the figure, the initial configuration consisted of a random mixture of water particles (cyan) and POPC lipids (red / pink heads, yellow tails). Driven by hydrophobic interactions, the lipids rapidly separated from the solvent, forming a transient perforated bilayer structure that minimized the exposure of the hydrophobic tails to water. The system eventually relaxed successfully into a stable, defect-free planar bilayer. The spontaneous formation of a layered phase from the chaotic initial state confirmed that the optimized parameters correctly captured the global free energy minimum of the lipid-water system.
[0119] The embodiments described above are merely preferred embodiments of the present invention and are not intended to limit the invention. Those skilled in the art can make various changes and modifications without departing from the spirit and scope of the invention. Therefore, all technical solutions obtained through equivalent substitution or transformation fall within the protection scope of the present invention.
Claims
1. An automatic optimization method for force field parameters of a biomembrane mesoscopic model with multi-system joint constraints, characterized in that, Specifically as follows: S1: Construct a parallel co-simulation system that treats multiple physical systems as a joint evaluation unit; S2: Implement parameter decoupling and hierarchical mapping strategies for the parallel co-simulation system described in S1, and discretize the global force field parameter set according to the interaction type between particles and the corresponding physical response characteristics. S3: Preset parameter boundaries of each space after discretization in S2 An internal randomized population of N particles is initialized, where each particle represents a set of potential force field parameters. The position vectors of the particles are initialized. and velocity vector ; S4: For each set of force field parameters in the population obtained in S3 or S7, the LAMMPS simulation engine is started simultaneously using the MPI parallel mechanism to calculate the parallel co-simulation system described in S1. S5: Post-process the simulated trajectory obtained in S4 to extract macroscopic physical properties and calculate physical constraint indicators; S6: Based on the deviation between the simulation results obtained in S5 and the preset experimental reference values, construct the loss function. ; S7: Update the individual historical best position of each particle based on the loss function value calculated according to the loss function described in S6. The global optimal position of the entire population According to the PSO speed update formula ,according to Update the next generation parameters; among them, Let K be the velocity vector of particle k at the (t+1)th iteration. Let K be the velocity vector of particle k at the t-th iteration. Let be the inertia weight for the t-th iteration. For individual learning factors, As a social learning factor; A random number between [0, 1], used to increase the randomness of individual searches; A random number between [0, 1], used to increase the randomness of the group search; Let be the current position vector of particle k at the t-th iteration, i.e., the current parameter solution; S8: Repeat S4~S7 until the maximum number of iterations is reached or the loss function converges to the preset threshold, and finally output the obtained optimal force field parameter set.
2. The automatic optimization method for force field parameters of a biomembrane mesoscopic model with multi-system joint constraints as described in claim 1, characterized in that, In S1, the parallel co-simulation system includes a water-air interface system for calibrating the fluid state equation and the surface tension of the water phase, an oil-water interface system for anchoring the hydrophobic interaction scale and the oil-water interface tension, and a phospholipid bilayer system for precisely controlling the structural and mechanical properties of the biofilm.
3. The automatic optimization method for force field parameters of a biomembrane mesoscopic model with multi-system joint constraints as described in claim 2, characterized in that, The water-air interface system is a vacuum-water-vacuum sandwich structure, the oil-water interface system is a water-oil-water sandwich structure, and the phospholipid bilayer system is a water-phospholipid bilayer-water sandwich structure.
4. The automatic optimization method for force field parameters of a biomembrane mesoscopic model with multi-system joint constraints as described in claim 1, characterized in that, In S2, the global force field parameter set is discretized into: a solvent interaction subspace characterizing the intrinsic properties of the fluid, an interface interaction subspace characterizing the properties of multiphase contact, and an intramembrane force field subspace characterizing the self-assembly properties of phospholipid molecules.
5. The automatic optimization method for force field parameters of a biomembrane mesoscopic model with multi-system joint constraints as described in claim 1, characterized in that, S5 includes: Calculation of water-air interfacial tension using pressure tensor and oil-water interfacial tension ; Calculate film thickness using density distribution ; By simulating box size Calculate the area occupied by monolipids ;in, It is the number of phospholipids in the system; By fitting the Helfrich-Canham theoretical spectrum to the discrete Fourier transform of the bilayer height field, Calculate the bending stiffness of the membrane ;in, It is Boltzmann's constant. It's temperature. It is tension. It is a two-dimensional wave vector The model; The number of water molecules that permeate into the hydrophobic core region of the membrane is counted to determine whether the membrane has ruptured; Calculate the diffusion coefficient of the lipid tail chain to determine whether the membrane is in a fluid state; The orientation order of the lipid tail chains is calculated to determine whether they conform to the characteristics of a liquid crystal phase.
6. The automatic optimization method for force field parameters of a biomembrane mesoscopic model with multi-system joint constraints as described in claim 1, characterized in that, In S6, the loss function The calculation formula is as follows: ； in, , , ; This is a soft constraint penalty function used to eliminate abnormal states caused by non-physical parameters; the subscripts sim and exp represent the simulated value and the experimental reference value, respectively. For the weighting coefficients of the loss term in the water-air interface system, For the weighting coefficients of the loss term in the water-oil interface system, For the weighting coefficients of the loss term in the phospholipid bilayer membrane interface system, The weighting coefficient for the area loss term of single lipids. The weighting coefficient for the lipid membrane thickness loss term. is the weighting coefficient for the lipid membrane bending modulus loss term, and j is the index of the soft constraint penalty function.
7. An automatic optimization system for force field parameters of a biomembrane mesoscopic model with multi-system joint constraints, characterized in that, include: The parallel co-simulation system construction module is used to build a parallel co-simulation system that treats multiple physical systems as a joint evaluation unit. The discretization module is used to implement parameter decoupling and hierarchical mapping strategies for the parallel co-simulation system, and to discretize the global force field parameter set according to the interaction type between particles and the corresponding physical response characteristics. The particle swarm optimization algorithm module is initialized to define the preset parameter boundaries of each space after discretization in the discretization processing module. An internal randomized population of N particles is initialized, where each particle represents a set of potential force field parameters. The position vectors of the particles are initialized. and velocity vector ; The simulated trajectory calculation module is used to simultaneously start the LAMMPS simulation engine using the MPI parallel mechanism to calculate each set of force field parameters in the population obtained by the initialization particle swarm optimization algorithm module or update iteration module, and to perform calculations on the parallel co-simulation system. The post-processing module is used to post-process the simulated trajectory obtained by the simulated trajectory calculation module, extract macroscopic physical properties and calculate physical constraint indicators; The loss function construction module is used to construct a loss function based on the deviation between the simulation results obtained by the post-processing module and the preset experimental reference values. ; The update iteration module is used to update the individual historical best position of each particle based on the loss function value calculated by the loss function. The global optimal position of the entire population According to the PSO speed update formula ,according to Update the next generation parameters; among them, Let K be the velocity vector of particle k at the (t+1)th iteration. Let K be the velocity vector of particle k at the t-th iteration. Let be the inertia weight for the t-th iteration. For individual learning factors, As a social learning factor; A random number between [0, 1], used to increase the randomness of individual searches; A random number between [0, 1], used to increase the randomness of the group search; Let be the current position vector of particle k at the t-th iteration, i.e., the current parameter solution; The results output module is used to repeatedly run the simulation trajectory calculation module, post-processing module, loss function construction module, and update iteration module in sequence until the maximum number of iterations is reached or the loss function converges to a preset threshold, and finally outputs the obtained optimal force field parameter set.
8. A computer program product comprising a computer program / instructions, characterized in that, When the computer program / instruction is executed by the processor, it can realize the automatic optimization method for force field parameters of biomembrane mesoscopic model with multi-system joint constraints as described in any one of claims 1 to 6.
9. A computer-readable storage medium, characterized in that, The storage medium stores a computer program, which, when executed by a processor, implements the automatic optimization method for force field parameters of a biomembrane mesoscopic model with multi-system joint constraints as described in any one of claims 1 to 6.
10. A computer electronic device, characterized in that, Including memory and processor; The memory is used to store computer programs; The processor is configured to, when executing the computer program, implement the automatic optimization method for force field parameters of a biomembrane mesoscopic model with multi-system joint constraints as described in any one of claims 1 to 6.