A method and system for simulating nanoparticle diffusion behavior in tissue interstitial spaces

Abstract

The application provides a simulation method and system for diffusion behavior of nanoparticles in tissue gaps, and belongs to the field of molecular dynamics simulation, and comprises the following steps: constructing a multi-component coarse-grained molecular model; performing energy minimization on the multi-component coarse-grained model by using a conjugate gradient method, and obtaining an initial configuration with stable energy by setting a convergence tolerance; maintaining a constant temperature condition by using a Langevin heat bath under a canonical ensemble, performing molecular dynamics equilibrium simulation on the initial configuration to obtain the diffusion behavior of nanoparticles in tissue gaps, and recording motion trajectory data of the nanoparticles; obtaining mean square displacement (MSD) of the nanoparticles according to the motion trajectory data, and plotting MSD-t curves by corresponding MSD and simulation time t; performing linear fitting on the MSD-t curves and solving, to obtain a diffusion coefficient of the nanoparticles in the tissue gaps. The application analyzes the interaction between the nanoparticles and the solid-liquid two-phase medium from the molecular level, and provides a powerful microscopic research tool for elucidating the diffusion mechanism of the nanoparticles.

Description

Technical Field

[0001] This invention belongs to the field of molecular dynamics simulation technology, specifically relating to a method and system for simulating the diffusion behavior of nanoparticles in tissue interstitial spaces. Background Technology

[0002] Drug-carrying nanoparticles (NPs), as key carriers for targeted therapy, still face the serious challenge of low delivery efficiency in clinical applications. Taking tumor treatment as an example, only about 0.7% of nanomedicines reach the tumor site after intravenous injection, and the amount of drug that ultimately enters the target cells to exert its effect is reduced to 0.0014% of the total injected amount. The main reason for this phenomenon is that the transport of nanoparticles in the interstitial space is severely restricted. In diseased tissues, nanomedicines can only migrate in the interstitial space by diffusion. However, the microenvironment of diseased tissues such as tumors undergoes significant changes: the extracellular matrix (ECM) density increases, tissue stiffness increases, and the viscosity of interstitial fluid increases, forming a complex physical barrier. Therefore, simulating the diffusion behavior of nanoparticles in the interstitial space is particularly important.

[0003] Current research on the diffusion behavior of nanoparticles in interstitial spaces remains at the theoretical and experimental stages. Most studies (such as those in PNAS) focus solely on the physical barrier effect of the solid ECM matrix, neglecting the microscopic interactions between the macromolecular components (concentration and length) in the interstitial fluid and the nanoparticles. This lack of comprehensive consideration of the synergistic influence of the solid-liquid two-phase media on diffusion behavior leads to an incomplete understanding of the diffusion mechanism. The classic Stokes-Einstein diffusion equation only applies to homogeneous, ideal, and simple media, failing to explain the anomalous diffusion behavior of nanoparticles in non-homogeneous and complex media such as interstitial spaces (e.g., the nonlinear change of the diffusion coefficient with environmental viscosity), thus limiting its theoretical guidance value. Existing experimental techniques lack sufficient spatial and temporal resolution, making it difficult to accurately capture the dynamic diffusion trajectory of nanoparticles at the molecular scale and quantify the microscopic interactions between nanoparticles and ECM / interstitial fluid macromolecules. Consequently, the dynamic mechanism of the diffusion process is difficult to directly observe and verify.

[0004] In summary, existing methods generally neglect key hydrodynamic factors such as interstitial fluid viscosity, making it difficult to accurately quantify the dynamic process of nanoparticle diffusion behavior in the interstitial space, which directly restricts the design and development of efficient nano-drug delivery systems. Summary of the Invention

[0005] To address the challenges of unclear diffusion mechanisms of nanoparticles in complex interstitial media and the limitations of traditional theories and experimental methods, this invention provides a method for simulating the diffusion behavior of nanoparticles in interstitial spaces.

[0006] To achieve the above objectives, the present invention provides the following technical solution: A method for simulating the diffusion behavior of nanoparticles in tissue interstitial spaces includes the following steps: A multi-component coarse-grained molecular model was constructed, comprising a polymer network simulating the extracellular matrix (ECM), nanoparticles serving as drug carriers, and polymer chains embedded in the pores of the polymer network simulating macromolecular components in interstitial fluid. The energy of the multi-component coarse-grained model is minimized using the conjugate gradient method, and an energy-stable initial configuration is obtained by setting a convergence tolerance. Under a canonical ensemble, a Langevin thermostat was used to maintain constant temperature conditions. Molecular dynamics equilibrium simulations of the diffusion behavior of nanoparticles in the interstitial space were performed on the initial configuration, and the motion trajectory data of the nanoparticles were recorded. The mean square displacement (MSD) of the nanoparticles is obtained based on the motion trajectory data. The MSD of the nanoparticles is plotted against the simulation time t to obtain the MSD-t curve. The MSD-t curve is linearly fitted to solve for the slope of the MSD-t curve. The diffusion coefficient of the nanoparticles in the interstitial space is determined based on the slope of the MSD-t curve.

[0007] Preferably, all components in the multi-component coarse-grained molecular model are composed of coarse-grained CG beads; The polymer network is an expandable topology structure formed by the cross-linking of multiple polymer chains, and its network size and porosity characteristics can be controlled by adjusting the lattice parameters; By adjusting the ratio of the total number of polymer chain CG beads to the total number of network CG beads within the range of 0 to 2, and adjusting the ratio of single chain length to network pore size within the range of 0 to 2, the interstitial environment of different viscosities from physiological to pathological states can be simulated. The nanoparticles are rigid structures composed of several CG beads linked by strong harmonic bonds, and their size and shape are configured according to the simulation requirements; the number of nanoparticles inserted in the multi-component coarse-grained molecular model is adjusted by adjusting the ratio of nanoparticles to the total volume of the system.

[0008] Preferably, after constructing the multi-component coarse-grained molecular model, the method further includes setting a corresponding molecular dynamics interaction field; the molecular dynamics interaction field includes bonding interactions and non-bonding interactions; the bonding interactions are described using harmonic potentials, including: The bond stretching interaction has the following potential energy function: In the formula, The potential energy is the bond length stretching energy. Let b be the bond stretching constant acting between particle i and particle j, and let b represent the bond stretching. The instantaneous distance between the particle pairs is [missing information]. To balance its bond length; The potential energy function of the bond-angle bending interaction is: In the formula, The potential energy is the bond angle bending energy. Let $\frac{i}{j}$ be the bending force constant acting on the bond angle formed by the sequentially connected particles $i$, $j$, and $k$. Indicates that the bond angle is bent. This is the instantaneous value of the bond angle. Its equilibrium value; The nonbonded interactions are described using the Lennard-Jones potential, whose potential function is: In the formula, For Lennard-Jones potential energy, and These are the interaction strength parameter between particle i and particle j and the characteristic distance when the potential energy is zero, respectively. The interaction force field sets different cutoff radii for different particle pairs, specifically including two types: Type I cutoff radius: its value is greater than This is used to preserve the full range of attractive forces in the Lennard-Jones potential in order to simulate long-range tunable interactions between particles; The second type of cutoff radius: its value is precisely set at the minimum value of the Lennard-Jones potential energy, and the potential energy is shifted to make the potential value at that point zero, thereby forming a modified potential field that retains only short-range repulsion.

[0009] Preferably, when using the conjugate gradient method to minimize the energy of the multi-component coarse-grained model and obtaining an energy-stable initial configuration by setting a convergence tolerance, both the energy convergence tolerance and the force convergence tolerance for energy minimization are set to 1e. -25 The maximum number of iterations required to minimize energy is 1,000,000. During the energy minimization process, the potential energy of the system, the size of the simulation box, the pressure tensor, and the total atomic potential energy are output every 10,000 steps.

[0010] Preferably, the process of maintaining constant temperature using a Langevin thermoelectric bath under a canonical ensemble to simulate the diffusion behavior of nanoparticles in the interstitial tissue under molecular dynamics equilibrium conditions, and recording the trajectory data of the nanoparticles, specifically involves: The multi-component coarse-grained model after energy minimization was subjected to equilibrium simulation under a canonical ensemble; a Langevin heat bath was used to maintain an isothermal environment at a reduction temperature of 1.0ε; the integration step size was set to 0.002τ, and the total simulation time was 10,000,000 steps; the thermodynamic parameters and motion trajectory data of the system were monitored and output.

[0011] Preferably, the thermodynamic parameters of the system include energy parameters, temperature parameters, pressure parameters, and system size parameters; the motion trajectory data includes the time series coordinates of the particles, the bonding and topological relationships of the particles, and the specific dynamic information of the nanoparticles.

[0012] Preferably, the diffusion coefficient of nanoparticles in the interstitial space is calculated according to Einstein's diffusion equation. D : ;in t represents the slope of the MSD-t curve.

[0013] The present invention also provides a simulation system for the diffusion behavior of nanoparticles in the interstitial space, comprising: The model building module is used to build a multi-component coarse-grained molecular model, which includes a polymer network simulating the extracellular matrix (ECM), nanoparticles as drug carriers, and polymer chains embedded in the pores of the polymer network to simulate macromolecular components in interstitial fluid. The model optimization module is used to minimize the energy of a multi-component coarse-grained model using the conjugate gradient method, and obtains an energy-stable initial configuration by setting the convergence tolerance. The simulation module is used to perform molecular dynamics equilibrium simulation of the diffusion behavior of nanoparticles in the interstitial space under constant temperature conditions using a Langevin heat bath under a canonical ensemble, and to record the motion trajectory data of the nanoparticles. The diffusion coefficient acquisition module is used to obtain the mean square displacement (MSD) of the nanoparticles based on the motion trajectory data, plot the MSD of the nanoparticles against the simulation time t to obtain the MSD-t curve, perform linear fitting on the MSD-t curve, solve for the slope of the MSD-t curve, and determine the diffusion coefficient of the nanoparticles in the interstitial space based on the slope of the MSD-t curve.

[0014] The present invention also provides a computer device, including a memory, a processor, and a computer program stored in the memory, wherein the processor executes the computer program to implement the steps of any one of the methods for simulating the diffusion behavior of the nanoparticles in the interstitial space.

[0015] The present invention also provides a computer-readable storage medium storing a computer program that, when loaded by a processor, is capable of executing any of the steps in the simulation method for the diffusion behavior of the nanoparticles in the interstitial space.

[0016] The method for simulating the diffusion behavior of nanoparticles in the interstitial space provided by this invention has the following beneficial effects: This invention first constructs a multi-component coarse-grained model comprising a polymer network (solid phase) containing an ECM, nanoparticles, and polymer chains simulating interstitial fluid macromolecules (liquid phase). By adjusting the concentration and length of the polymer chains to simulate different viscosity environments, it compensates for the neglect of interstitial fluid viscosity factors. A stable initial configuration is then obtained through energy minimization using the conjugate gradient method, laying a reliable foundation for subsequent simulations. Subsequently, isothermal molecular dynamics simulations are conducted under a canonical ensemble using a Langevin thermoelectric bath, accurately recording motion trajectories and thermodynamic parameters, overcoming the limitations of experimental resolution. Finally, the diffusion coefficient (MSD) is calculated from the motion trajectory data and plotted as a curve. Through linear fitting combined with the Einstein equation, the diffusion coefficient is solved, achieving quantitative analysis of diffusion behavior. Each step is closely integrated, comprehensively considering the synergistic effect of the solid-liquid phases while precisely quantifying the dynamic diffusion process of nanoparticles, clearly revealing its diffusion mechanism. This provides a reliable microscopic research tool for the design and development of efficient nanoparticle drug delivery systems, effectively solving the core technical bottlenecks of existing methods. Attached Figure Description

[0017] To more clearly illustrate the embodiments and design schemes of the present invention, the accompanying drawings required for this embodiment will be briefly described below. The drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0018] Figure 1 A flowchart illustrating the simulation method for the diffusion behavior of nanoparticles in the interstitial space provided by this invention; Figure 2 This is a schematic diagram of the multi-component coarse-grained molecular models constructed in Examples 1 and 2; wherein, Figure 2 (a) is a schematic diagram of a polymer network. Figure 2 (b) is a schematic diagram of nanoparticles. Figure 2 (c) is a schematic diagram of the polymer chain; Figure 3 This is a simulation diagram of the system under equilibrium state, simulating the diffusion process of NPs in the interstitial space in Example 1. Figure 4 This is a simulation diagram of the system under equilibrium state, simulating the diffusion process of NPs in the interstitial space in Example 2. Figure 5 The mean square displacement-time curves of NPs under different polymer chain concentrations in Example 1 of the present invention are shown, wherein the ratio of the total number of CGs in polymer chains to the total number of CGs in the network are [0, 0.25, 0.5, 1, 2], respectively. Figure 6The mean square displacement-time curves of NPs under different polymer chain lengths in Example 2 of the present invention are shown, wherein the ratios of polymer chain length to network pore size are [0, 0.5, 1, 1.5, 2]. Figure 7 This is a diffusion coefficient diagram of NPs under different polymer chain concentrations in Example 1 of the present invention; Figure 8 This is a diffusion coefficient diagram of NPs under different polymer chain lengths in Example 2 of the present invention. Detailed Implementation

[0019] To enable those skilled in the art to better understand and implement the technical solutions of the present invention, the present invention will be described in detail below with reference to the accompanying drawings and specific embodiments. The following embodiments are only used to more clearly illustrate the technical solutions of the present invention and should not be construed as limiting the scope of protection of the present invention.

[0020] The interstitial space consists of the solid extracellular matrix (ECM) and interstitial fluid. This invention uses the ECM as the solid-phase framework and explicitly introduces macromolecules that determine viscosity from the interstitial fluid (while employing an implicit water model to save computational effort) to construct a solid-liquid composite system, thereby studying the diffusion behavior of nanoparticles in this environment. Molecular dynamics simulations can simultaneously analyze the interaction between nanoparticles and both solid and liquid phases at the molecular scale, providing a powerful microscopic research tool for elucidating their diffusion mechanisms. Therefore, this invention, based on the molecular dynamics software LAMMPS, constructs an ECM / NPs / macromolecule system model, selects a suitable bio-nanomaterial composite force field, sets appropriate simulation parameters and a systematic methodology, obtains the mean square displacement (MSD) of nanoparticles in the composite system, and calculates the diffusion coefficient of nanoparticles in the interstitial space using the Einstein equation.

[0021] By simulating the regulation of interstitial fluid viscosity by changing the concentration of polymer chains (adjusting the ratio of the total number of CG beads in polymer chains to the total number of CG beads in the network) and the length (adjusting the ratio of single chain length to network pore size), the influence of interstitial fluid viscosity changes on the diffusion behavior of nanoparticles in solid-liquid composite systems is systematically studied. The diffusion process can be analyzed by combining mean square displacement, diffusion coefficient, and motion trajectory in the simulation results, thereby exploring the diffusion law of nanoparticles (NPs) in the interstitial space.

[0022] Based on this, the present invention provides a method for simulating the diffusion behavior of nanoparticles in tissue interstitial spaces, specifically a method for simulating the diffusion behavior of nanoparticles in tissue interstitial spaces using molecular dynamics, such as... Figure 1 As shown, it includes the following steps: Step 1: Constructing a multi-component coarse-grained molecular model: A regular polymer network is generated using MATLAB algorithms to simulate the ECM structure. Nanoparticles (NPs) of the desired configuration are generated as drug carrier models. Freely distributed polymer chains are embedded in the pores of the polymer network to simulate soluble macromolecular components in interstitial fluid. The polymer network, nanoparticles, and polymer chains together constitute the multi-component coarse-grained molecular model. All components in the model are composed of coarse-grained beads (CG beads). Furthermore, the bonding and non-bonding interaction force fields required for the molecular dynamics model are defined.

[0023] Furthermore, in the multi-component coarse-grained molecular model constructed in step 1, the polymer network simulating the ECM is defined as a structure composed of multiple cross-linked polymer chains. It is an scalable topology, and its network size and porosity characteristics can be controlled by adjusting the lattice parameters. The simulation space of the multi-component coarse-grained molecular model adopts a cubic periodic boundary box, the size of which is configured according to the network scale.

[0024] Furthermore, the polymer chain parameters introduced in step 1 satisfy the following proportional relationships: the ratio of the total number of chain CG beads to the total number of network CG beads is adjustable within a range of 0 to 2, and the ratio of the length of a single chain to the size of the network pores is adjustable within a range of 0 to 2. By adjusting the above proportional parameters, it is possible to effectively simulate the interstitial environment of different viscosities from physiological to pathological states.

[0025] Furthermore, the NPs constructed in step 1 are rigid structures composed of several CG beads linked by strong harmonic bonds, and their size and shape can be configured according to simulation requirements; the number of NPs implanted in the model can be controlled and adjusted by adjusting the ratio of NPs to the total volume of the system.

[0026] Furthermore, the bonding interactions in step 1 are described using harmonic potentials, specifically including: (i) Bond stretching interaction, whose potential energy function is defined by equation (1), is used to describe the stretching vibration of the chemical bond between two bonding particles: (1) In the formula, The potential energy is the bond length stretching energy. Let be the bond stretching constant acting between particle i and particle j, with the superscript b indicating bond stretching. The instantaneous distance between the particle pairs is [missing information]. To balance its bond length.

[0027] (ii) Bond angle bending interaction, whose potential energy function is defined by formula (2), is used to describe the bending vibration of the bond angle formed by three consecutive bonding particles: (2) In the formula, The potential energy is the bond angle bending energy. The superscript represents the bending force constant acting on the bond angle formed by the sequentially connected particles i, j, and k. Indicates that the bond angle is bent. This is the instantaneous value of the bond angle. Its equilibrium value.

[0028] Furthermore, the non-bonded interactions in step 1 are described by the Lennard-Jones potential as van der Waals interactions between particles, and its potential energy function is defined by equation (3): (3) In the formula, For Lennard-Jones potential energy, and , respectively, represent the interaction strength parameter between particle i and particle j and the characteristic distance when the potential energy is zero.

[0029] To precisely control physical processes and improve computational efficiency, differentiated cutoff radii are set for the Lennard-Jones potential function. Specific schemes include: (i) For the interactions between nanoparticles and ECMs, and between nanoparticles and macromolecules, a first cutoff radius is used to fully cover their van der Waals attraction and repulsion.

[0030] (ii) For interactions within the ECM, macromolecules, nanoparticles, and between the ECM and macromolecules, a second cutoff radius is adopted, which is set to retain only the purely repulsive portion of the interaction. Preferably, the second cutoff radius is set at the Lennard-Jones potential minimum.

[0031] Furthermore, the interaction potential function and force field setting in step 1 also include applying an elastic restoring force to the polymer network to constrain it to its initial spatial position.

[0032] Step 2: Model Optimization: The conjugate gradient method is used to minimize the energy of the multi-component coarse-grained model, with the energy variation tolerance and maximum force tolerance both set to 1e. -25 The maximum number of iterations is 1,000,000. During the minimization process, the system potential energy, simulation box size, pressure tensor, and total atomic potential energy are output every 10,000 steps. The energy minimization eliminates unreasonable molecular contacts and stresses within the system, thereby obtaining an energy-stable initial configuration.

[0033] Step 3: Molecular dynamics simulation: Molecular dynamics equilibrium simulation was performed on the system after energy minimization under the canonical (NVT) ensemble. A Langevin thermoelectric bath was used to maintain an isothermal environment at a reduction temperature of 1.0ε. The integration step size was set to 0.002τ and the total simulation time was 10,000,000 steps. During the equilibrium simulation, the trajectory file, system thermodynamic parameters and motion trajectory data were monitored and recorded.

[0034] Among them, the thermodynamic parameters of the system are key physical quantities describing the macroscopic physical state of the system during molecular dynamics simulations. They reflect the dynamic equilibrium characteristics of the simulated system in terms of energy, temperature, pressure, etc., and are the core basis for judging whether the simulation has reached a stable state and verifying the rationality of the model. The specific thermodynamic parameters of the system include:

[0035] Energy parameters: total potential energy of the system, total potential energy of atoms (including the contributions of bonding energy, non-bonding energy, etc.).

[0036] Temperature parameters: macroscopic temperature of the system (maintained at the approximation temperature of 1.0ε using a Langevin heat bath).

[0037] Pressure parameters: system pressure, pressure tensor (reflecting the pressure distribution of the system in three-dimensional space).

[0038] System dimensional parameters: The dimensions of the simulated periodic boundary box. In the NVT simulation of this embodiment, the box dimensions are fixed, and the system volume remains constant.

[0039] The main purposes of monitoring the thermodynamic parameters of the system are as follows: 11. Verify the stability of the simulation system: Monitor thermodynamic parameters in real time during the simulation. If the parameters (such as temperature, pressure, and potential energy) fluctuate little and tend to be constant over a long period of time, it indicates that the system has reached thermodynamic equilibrium and the subsequent diffusion behavior data is valid. If the parameters continue to drift (such as the temperature deviating from 1.0ε), the heat bath parameters or simulation settings need to be adjusted to avoid the non-equilibrium state affecting the results.

[0040] 2. Optimize the model and simulation parameters: If the energy parameters are abnormally high, it may be due to unreasonable molecular contact in the initial configuration. You need to return to step 2 (model optimization) and minimize the energy again. If the pressure deviates from the reasonable range, you can adjust the periodic boundary conditions or the cutoff radius of the interaction potential.

[0041] 3. Supplementary analysis of diffusion mechanism: Changes in thermodynamic parameters can help explain the differences in the diffusion coefficient of nanoparticles. For example, an increase in system potential energy may be accompanied by enhanced intermolecular interactions, indirectly confirming the inhibitory effect of increased viscous resistance on nanoparticle diffusion (e.g., in the examples, when the concentration / length of macromolecules increases, the potential energy is positively correlated with viscosity, and the diffusion coefficient decreases).

[0042] Specifically, the motion trajectory data consists of the spatial position information of all particles (polymer network particles, nanoparticles, and polymer chain particles) recorded during the simulation at different time steps. This data serves as the raw data for subsequent calculations of the mean square displacement (MSD) and analysis of nanoparticle diffusion paths. The motion trajectory data specifically includes:

[0043] Time-series coordinates of particles: the three-dimensional spatial coordinates (X, Y, Z) of nanoparticles, polymer network CG beads, and polymer chain CG beads at each time step.

[0044] Particle bonding and topological relationships: the bonding state between particles that changes over time (such as the stretching and bending of polymer chains).

[0045] Unique dynamic information of nanoparticles: the position of the center of mass of the nanoparticle, the direction of motion, the instantaneous velocity, etc. (extracted separately for diffusion behavior analysis).

[0046] The diffusion path of nanoparticles can be visualized by using trajectory files (such as the dump.Lammpstrj file output by LAMMPS), allowing us to observe whether their movement in the pores of the polymer network is hindered by the polymer chains, and to intuitively verify the conclusion that "increased macromolecular concentration / length leads to impeded diffusion".

[0047] Step 4: Plot the MSD-t curve.

[0048] Based on the motion trajectory data, the position coordinates of nanoparticles at different time points are extracted, and then substituted into the mean square displacement formula (4) to calculate the MSD value at different time intervals t, and then the MSD-t curve is plotted.

[0049] (4) In the formula, N The total number of particles, yes The position of particle i at time i This represents the initial position of particle i. The single brackets <> indicate multiple sampling observations (for time intervals). t (Sampling) and average.

[0050] Step 5: Calculate the diffusion coefficient of nanoparticles in the interstitial space model: Perform linear fitting on the MSD-t curve and solve for the slope of the MSD-t curve. k The diffusion coefficient of nanoparticles in the interstitial space was calculated based on Einstein's diffusion equation. D .

[0051] In three-dimensional space, the diffusion coefficient is... D With the slope of the MSD-t curve k The conversion relationship is as follows .

[0052] This invention, through molecular dynamics simulations, overcomes the limitations of the traditional Stokes-Einstein diffusion formula in non-homogeneous media such as interstitial spaces, accurately describing the anomalous diffusion behavior of nanoparticles. Furthermore, while existing studies consider the influence of the solid matrix in isolation, this invention, for the first time, comprehensively considers the synergistic effect of the solid and liquid phases (including the crucial interstitial fluid viscosity) on diffusion behavior in the interstitial space. It reveals the interaction mechanism between nanoparticles and multiphase media at the molecular scale, providing a more precise and reliable microscopic research method for addressing the bottleneck of nanomedicine delivery efficiency in complex tissue barriers.

[0053] Example 1 This embodiment provides a method for simulating the diffusion behavior of nanoparticles in tissue interstitial spaces using molecular dynamics. The specific steps for using LAMMPS software to perform a simulation on a computer or server are as follows: Step 1: Construct a multi-component coarse-grained model, specifically as follows: A regular porous periodic polymer network with a lattice constant of 10σ, containing 4 repeating units, a total of 1840 CG beads, and a size of 40σ×40σ×40σ was constructed. Twenty gold nanoparticles were inserted, and 200 free macromolecular polymer chains with a length of 10σ were randomly generated to construct a multi-component coarse-grained model.

[0054] Step 1.1: Generate the polymer network model This paper describes the automated construction of polymer networks using MATLAB algorithms, ultimately outputting standard GROMACS format coordinate files (GRO files). The implementation process includes three core stages:

[0055] (i) In the node generation stage, by setting the parameter system with period number n=4, lattice constant a=10 and boundary offset p=5, a triple loop is used to traverse the three-dimensional space from 0 to n*a. For each coordinate point (x,y,z), the offset (xp,yp,zp) is calculated and modulo the lattice constant a to obtain (modx,mody,modz). When the lattice condition (modx+mody==0) or (modx+modz==0) or (mody+modz==0) is satisfied, the point is determined as a network node and its coordinates are recorded in matrix A. At the same time, a complete coordinate matrix B containing atomic numbers is generated.

[0056] (ii) In the key identification stage, all node pairs are traversed by a double loop to calculate their Euclidean distance. When the distance is exactly equal to 1 unit length, the key relationship is established and recorded in matrix S.

[0057] (iii) During the file output stage, the coordinate data of matrix B is synchronously written to atom.txt and the general text network10.gro structure file. The gro file strictly follows the format requirements of molecular dynamics software and includes system description, number of atoms, detailed atom information and box size of 40σ×40σ×40σ. At the same time, the bond relationship recorded in matrix S is written to the bond.txt file.

[0058] Step 1.2: Generate nanoparticle model The nanoparticles were constructed using MATLAB algorithms. A three-dimensional mathematical model of gold crystal was established based on the face-centered cubic lattice features, and an initial atomic lattice was generated using a supercell extension algorithm. Subsequently, a geometric screening method was used, with the criterion |x|+|y|+|z|≤L of the truncated octahedron mathematical expression as the standard for spatial screening of atomic coordinates. The final number of atoms was precisely controlled to 147 by adjusting the size parameter L. After the structural screening was completed, a coordinate file NP.gro in GROMACS format was generated through the coordinate output module. At the same time, topological connection information was automatically generated based on the bonding criterion that the interatomic spacing is less than 3.2 Å and exported as a bonding file 2.txt.

[0059] Step 1.3: Construct the initial configuration of the NP-ECM composite simulation system The pre-built nanoparticle coordinate file NP.gro and the polymer network structure file network10.gro were submitted to a high-performance computing cluster. Using GROMACS's insertion function, a geometric insertion command was executed to randomly and non-overlappingly insert 20 NPs into the pores of the polymer network. After the cluster completed the calculation, it output a composite structure file network10-np20.gro containing the coordinates of all atoms. Subsequently, the atomic coordinate information (atom number and X, Y, Z coordinates) from this .gro file was extracted and formatted into a standard four-column data format, saved as the coordinate file 1.txt, and used as one of the core inputs for subsequent MATLAB programs.

[0060] Step 1.4: Integration of multi-component data and insertion of chains Place the following four core data files in the MATLAB working directory: 1.txt contains the coordinates of all atoms in the NP-ECM composite system, 2.txt defines the bonding relationships within the NP, atom.txt records the basic atomic information of the polymer network, and bond.txt defines the bonding topology of the polymer network. In the parameter settings area of ​​the MATLAB integrated program, set the chain number parameter to 200 and the chain length parameter to 10. After running the program, the system first reads atom.txt and bond.txt to construct the polymer network skeleton, then loads the atomic coordinates from 1.txt and the NP bonding information from 2.txt, and forms the coordinate matrix of the NP-polymer network binary system through atomic number remapping. Then, based on the set parameters of 200 chains and 10 CG beads per chain, the program executes a chain generation algorithm within the simulation box: randomly generating starting coordinates in the box space, placing chain CG beads sequentially along a random unit vector direction with a fixed step size, and detecting coordinate conflicts in real time through a spatial repulsion function. Finally, bond connections and bond corner topologies are automatically created for all successfully generated chains, completing the construction of a three-component system containing 200 chains, each with 10 CG beads. In the series of simulations in this embodiment, the chain number parameter is set to multiple different values ​​to generate several comparative models with ratios of 0, 0.25, 0.5, 1, and 2 between the total number of CG beads in the polymer chains and the total number of CG beads in the network.

[0061] Step 1.5: Automated generation of LAMMPS data files The program outputs all integrated model information in data format according to the requirements of LAMMPS data files, named 200chain.data. This file contains the following three parts: The first part is atomic information, which records in detail the atom ID, molecular ID, atom type, and three-dimensional coordinate data of all ECM atoms, NP atoms, and chain CG beads; the second part is topological information, which fully describes all bond relationships of polymer network bonds, NP internal bonds, and chain bonds, and also includes the bond angle information of the chains; the third part is force field parameters and simulation box settings, which clearly gives the specific parameters of atomic mass, bond type, angle type, and the precise size of the simulation box.

[0062] Step 1.6: Define the interaction potential function and force field This invention is based on the Lennard-Jones reduced dimensionless unit system and uses LAMMPS input files to precisely define the interactions in the system: (12) The resonant kinetic potential is used to describe the bond interaction. The bond potential type is defined by bond_style harmonic and three types of bond parameters are set: Type 1 (polymer network bond) force constant 23, equilibrium bond length 1.0, Type 2 (NP internal bond) force constant 200, equilibrium bond length 0.47, and Type 3 (chain bond) force constant 10, equilibrium bond length 1.0.

[0063] (ii) The resonant kinetic potential is used to describe the bond angle interaction. The bond angle potential type is defined by angle_style harmonic, and the bond angle force constant is set to 10 and the equilibrium bond angle is 180 degrees.

[0064] (iii) The Lennard-Jones potential function is used to describe nonbonded interactions. The potential function type and global cutoff radius are defined using pair_style lj / cut2.0, and the interaction parameters between each atom type are configured according to Table 1. At the same time, pair_modify shift yes is enabled to ensure energy continuity. Table 1 Interaction parameters between different atom types Step 2: Model optimization, specifically: This invention minimizes the system energy of a composite system to eliminate local stress and unreasonable contact in the initial configuration: setting the energy convergence tolerance to 1e. -25 The force convergence tolerance is 1e -25 The maximum number of iterations is 500,000. The conjugate gradient algorithm is used for optimization via the min_style command, and the changes in key parameters such as system potential energy, box size, pressure, and total atomic energy are monitored in real time via the thermo and thermo_style commands. The minimization process is executed via the minimize command to ensure that the system reaches the lowest energy stable state.

[0065] Step 3: Molecular dynamics simulation, specifically: Step 31: Assign initial velocities to all atoms using the velocity command, and maintain the system at a constant temperature using a Langevin thermometer to achieve molecular dynamics simulation under a canonical ensemble.

[0066] Step 32: Define the mean square displacement of the centroid of the gold nanoparticle swarm using the compute command, set the thermodynamic output parameters using the thermo_style command, and save the system trajectory file at fixed time intervals using the dump command.

[0067] Step 33: Reset the time step using the reset_timestep command, run the molecular dynamics simulation for the specified number of steps, remove the constraints after the simulation is complete, and save the system restart file.

[0068] Step 34: After completing all the aforementioned simulation parameter settings, this invention submits the three core files prepared in the previous steps—the LAMMPS data file (200chain.data) containing complete architecture information, the input script (in.txt) containing all interaction parameters and simulation procedures, and the cluster job scheduling script (lmp_jobs.lsf)—to the high-performance computing cluster. The job scheduling script specifies the computing resource requirements, including the number of nodes, the number of cores, the runtime, and the LAMMPS execution path. Upon receiving the task, the cluster automatically allocates computing resources and executes energy minimization and molecular dynamics simulations sequentially according to the procedures set in the input script. It then outputs the trajectory files, thermodynamic data, and restart files generated during the simulation to a designated directory, completing the entire computation task.

[0069] Step 4: Plot the MSD-t curve, specifically: After the cluster completes the molecular dynamics simulation, the trajectory file (dump.Lammpstrj) and log file (log.200chain.txt) are downloaded first, and the time step and mean square displacement data of gold nanoparticles, which were recorded in advance by the compute command, are extracted from them. Then, the MSD-t curve is plotted with the simulation time as the x-axis and the mean square displacement as the y-axis.

[0070] Step 5: Calculate the diffusion coefficient of nanoparticles in the interstitial space model, specifically: The slope k was obtained by least-squares fitting of the linear interval of the MSD-t curve, and the diffusion coefficient D of nanoparticles in the interstitial space was calculated according to the Einstein diffusion equation: .

[0071] The diffusion process of nanoparticles in the interstitial spaces was simulated using the above method. The completed model is shown in the diagram below. Figure 2 As shown.

[0072] The diffusion process of nanoparticles in the interstitial space was simulated using the above method. The simulation diagram of the system under equilibrium conditions is shown below. Figure 3 As shown. Figure 3 This is a simulation diagram of the system reaching equilibrium when the chain concentration ratio (total number of chain CG beads / total number of network CG beads) is 1.

[0073] Molecular dynamics simulations were performed on multiple models with different chain concentrations. The MSD-t curves generated after the molecular dynamics simulations are shown below. Figure 5As shown. Calculations showed that, with a polymer chain length of 10σ and a total number of network CG beads of 1840, the nanoparticle diffusion coefficients at different macromolecular chain concentrations (the ratio of the total number of polymer chain CG beads to the total number of network CG beads was 0, 0.25, 0.5, 1, and 2, respectively) were 1.57 × 10⁻⁶. -2 1.43×10 -2 1.27×10 -2 1.23×10 -2 1.07×10 -2 ).

[0074] Figure 7 The diffusion coefficients of nanoparticles (NPs) are shown below under different polymer chain concentrations. Comparative analysis reveals that the diffusion coefficients of NPs decrease significantly with increasing polymer chain concentration. This is mainly because the increased polymer chain concentration strengthens the intermolecular interactions within the system. This enhanced interaction directly manifests as an increase in the macroscopic viscosity of the interstitial fluid, thereby significantly inhibiting the diffusion of nanoparticles through enhanced viscous resistance.

[0075] Example 2 This invention provides a method for simulating the diffusion behavior of nanoparticles in tissue interstitial spaces. The simulation is performed on a computer or server using LAMMPS software. The specific implementation steps are as follows: Step 1: Construct a multi-component coarse-grained model, specifically as follows: A regular porous periodic polymer network with a lattice constant of 10σ, containing 4 repeating units, a total of 1840 CG beads, and a size of 40σ×40σ×40σ was constructed. Twenty gold nanoparticles were inserted, and 200 free polymer chains with a length of 15σ were randomly generated to construct a multi-component coarse-grained model.

[0076] Step 1.1: Generate the polymer network model This paper describes the automated construction of polymer networks using MATLAB algorithms, ultimately outputting standard GROMACS format coordinate files (GRO files). The implementation process includes three core stages:

[0077] (13) In the node generation stage, by setting the parameter system with period number n=4, lattice constant a=10 and boundary offset p=5, the three-dimensional space from 0 to n*a is traversed by a triple loop. For each coordinate point (x,y,z), its offset (xp,yp,zp) is calculated and modulo the lattice constant a to obtain (modx,mody,modz). When the lattice condition (modx+mody==0) or (modx+modz==0) or (mody+modz==0) is satisfied, the point is determined as a network node and its coordinates are recorded to matrix A. At the same time, a complete coordinate matrix B containing atomic numbers is generated.

[0078] (14) In the link identification stage, all node pairs are traversed by double loops to calculate their Euclidean distance. When the distance is exactly equal to 1 unit length, the link relationship is established and recorded in matrix S.

[0079] (15) During the file output stage, the coordinate data of matrix B is synchronously written to atom.txt and the general text network10.gro structure file. The gro file strictly follows the format requirements of molecular dynamics software and includes system description, number of atoms, detailed atom information and box size of 40σ×40σ×40σ. At the same time, the bond connection relationship recorded in matrix S is written to the bond.txt file.

[0080] Step 1.2: Generate nanoparticle model The nanoparticles were constructed using MATLAB algorithms. A three-dimensional mathematical model of gold crystal was established based on the face-centered cubic lattice features, and an initial atomic lattice was generated using a supercell extension algorithm. Subsequently, a geometric screening method was used, with the criterion |x|+|y|+|z|≤L of the truncated octahedron mathematical expression as the standard for spatial screening of atomic coordinates. The final number of atoms was precisely controlled to 147 by adjusting the size parameter L. After the structural screening was completed, a coordinate file NP.gro in GROMACS format was generated through the coordinate output module. At the same time, topological connection information was automatically generated based on the bonding criterion that the interatomic spacing is less than 3.2 Å and exported as a bonding file 2.txt.

[0081] Step 1.3: Construct the initial configuration of the NP-ECM composite simulation system The pre-built nanoparticle coordinate file NP.gro and the polymer network structure file network10.gro were submitted to a high-performance computing cluster. Using GROMACS's insertion function, a geometric insertion command was executed to randomly and non-overlappingly insert 20 NPs into the pores of the polymer network. After the cluster completed the calculation, it output a composite structure file network10-np20.gro containing the coordinates of all atoms. Subsequently, the atomic coordinate information (atom number and X, Y, Z coordinates) from this .gro file was extracted and formatted into a standard four-column data format, saved as the coordinate file 1.txt, and used as one of the core inputs for subsequent MATLAB programs.

[0082] Step 1.4: Integration of multi-component data and insertion of chains Place the following four core data files in the MATLAB working directory: 1.txt contains the coordinates of all atoms in the NP-ECM composite system, 2.txt defines the bonding relationships within the NP, atom.txt records the basic atomic information of the polymer network, and bond.txt defines the bonding topology of the polymer network. In the parameter settings area of ​​the MATLAB integrated program, set the chain number parameter to 200 and the chain length parameter to 15. After running the program, the system first reads atom.txt and bond.txt to construct the polymer network skeleton, then loads the atomic coordinates from 1.txt and the NP bonding information from 2.txt, and forms the coordinate matrix of the NP-polymer network binary system through atomic number remapping. Then, based on the set parameters of 200 chains and 15 CG beads per chain, the program executes a chain generation algorithm within the simulation box: randomly generating starting coordinates in the box space, placing chain CG beads sequentially along a random unit vector direction with a fixed step size, and detecting coordinate conflicts in real time through a spatial repulsion function. Finally, the bonding and bond angle topology relationships are automatically created for all successfully generated chains, completing the construction of a three-component system containing 200 chains, each with 15 CG beads. In the series of simulations in this embodiment, the chain length parameter is set to multiple different values ​​to generate several comparative models with polymer chain length to network pore size ratios of 0, 0.5, 1, 1.5, and 2.

[0083] Step 1.5: Automated generation of LAMMPS data files The program outputs all integrated model information in data format according to the requirements of LAMMPS data files, named 15Lchain.data. This file contains the following three parts: The first part is atomic information, which records in detail the atom ID, molecular ID, atom type, and three-dimensional coordinate data of all ECM atoms, NP atoms, and chain CG beads; the second part is topological information, which fully describes all bond relationships of polymer network bonds, NP internal bonds, and chain bonds, and also includes the bond angle information of the chains; the third part is force field parameters and simulation box settings, which clearly gives the specific parameters of atomic mass, bond type, angle type, and the precise size of the simulation box.

[0084] Step 1.6: Define the interaction potential function and force field This invention is based on the Lennard-Jones reduced dimensionless unit system and uses LAMMPS input files to precisely define the interactions in the system: (i) The bonding interaction is described by resonant kinetic potential. The bond potential type is defined by bond_style harmonic and three types of bond parameters are set: type 1 (polymer network bond) force constant 23, equilibrium bond length 1.0, type 2 (NP internal bond) force constant 200, equilibrium bond length 0.47, and type 3 (chain bond) force constant 10, equilibrium bond length 1.0.

[0085] (ii) The resonant kinetic potential is used to describe the bond angle interaction. The bond angle potential type is defined by angle_style harmonic, and the bond angle force constant is set to 10 and the equilibrium bond angle is 180 degrees.

[0086] (iii) The Lennard-Jones potential function is used to describe nonbonded interactions. The potential function type and global cutoff radius are defined by pair_style lj / cut2.0, and the interaction parameters between each atom type are configured according to Table 2. At the same time, pair_modify shift yes is enabled to ensure energy continuity.

[0087] Table 2 Interaction parameters between different atom types Step 2: Model optimization, specifically: This invention minimizes the system energy of a composite system to eliminate local stress and unreasonable contact in the initial configuration: setting the energy convergence tolerance to 1e. -25 The force convergence tolerance is 1e -25The maximum number of iterations is 500,000. The conjugate gradient algorithm is used for optimization via the min_style command, and the changes in key parameters such as system potential energy, box size, pressure, and total atomic energy are monitored in real time via the thermo and thermo_style commands. The minimization process is executed via the minimize command to ensure that the system reaches the lowest energy stable state.

[0088] Step 3: Molecular dynamics simulation, specifically: Step 31: Assign initial velocities to all atoms using the velocity command, and maintain the system at a constant temperature using a Langevin thermometer to achieve molecular dynamics simulation under a canonical ensemble.

[0089] Step 32: Define the mean square displacement of the centroid of the gold nanoparticle swarm using the compute command, set the thermodynamic output parameters using the thermo_style command, and save the system trajectory file at fixed time intervals using the dump command.

[0090] Step 33: Reset the time step using the reset_timestep command, run the molecular dynamics simulation for the specified number of steps, remove the constraints after the simulation is complete, and save the system restart file.

[0091] Step 34: After completing all the aforementioned simulation parameter settings, this invention submits the three core files prepared in the previous steps—the LAMMPS data file (15Lchain.data) containing complete architecture information, the input script (in.txt) containing all interaction parameters and simulation procedures, and the cluster job scheduling script (lmp_jobs.lsf)—to the high-performance computing cluster. The job scheduling script specifies the computing resource requirements, including the number of nodes, cores, runtime, and LAMMPS execution path. Upon receiving the task, the cluster automatically allocates computing resources and executes energy minimization and molecular dynamics simulations sequentially according to the procedures set in the input script. It then outputs the trajectory files, thermodynamic data, and restart files generated during the simulation to a designated directory, completing the entire computation task.

[0092] Step 4: Plot the MSD-t curve, specifically: After the cluster completes the molecular dynamics simulation, the trajectory file (dump.Lammpstrj) and log file (log.15Lchain.txt) are downloaded first, and the time step and mean square displacement data of gold nanoparticles, which were recorded in advance by the compute command, are extracted from them. Then, the MSD-t curve is plotted with the simulation time as the x-axis and the mean square displacement as the y-axis.

[0093] Step 5: Calculate the diffusion coefficient of nanoparticles in the interstitial space model, specifically: The slope k was obtained by least-squares fitting of the linear interval of the MSD-t curve, and the diffusion coefficient D of nanoparticles in the interstitial space was calculated according to the Einstein diffusion equation. .

[0094] The diffusion process of nanoparticles in the interstitial spaces was simulated using the above method. The completed model is shown in the diagram below. Figure 2 As shown.

[0095] The diffusion process of nanoparticles in the interstitial space was simulated using the above method. The simulation diagram of the system under equilibrium conditions is shown below. Figure 4 As shown. Figure 4 This is a simulation diagram of the system reaching equilibrium when the ratio of polymer chain length to network pore size is 1.5.

[0096] Molecular dynamics simulations were performed on multiple models with different chain lengths. The MSD-t curves generated after the molecular dynamics simulations are shown below. Figure 6 As shown. Calculations showed that, with 200 polymer chains and 1840 CG beads in the network, the nanoparticle diffusion coefficients for different polymer chain lengths (the ratio of polymer chain length to network pore size was 0, 0.5, 1, 1.5, and 2) were 1.57 × 10⁻⁶. -2 1.53×10 -2 1.45×10 -2 1.25×10 -2 1.18×10 -2 )

[0097] Figure 8 The diffusion coefficients of nanoparticles (NPs) are shown for different polymer chain lengths. Comparative analysis reveals that the diffusion coefficients of NPs decrease significantly with increasing polymer chain length. This is mainly because the increased polymer chain length enhances the strength of intermolecular interactions within the system. This enhanced interaction directly manifests as an increase in the macroscopic viscosity of the interstitial fluid, thereby significantly inhibiting the diffusion of nanoparticles through enhanced viscous resistance.

[0098] Based on the same inventive concept, the present invention also provides a simulation system for the diffusion behavior of nanoparticles in tissue interstitial spaces, including a model building module, a model optimization module, and a diffusion coefficient acquisition module.

[0099] Specifically, the model building module is used to construct a multi-component coarse-grained molecular model, which includes a polymer network simulating the extracellular matrix (ECM), nanoparticles serving as drug carriers, and polymer chains embedded in the pores of the polymer network to simulate macromolecular components in interstitial fluid.

[0100] The model optimization module is used to minimize the energy of a multi-component coarse-grained model using the conjugate gradient method, and obtains an energy-stable initial configuration by setting a convergence tolerance.

[0101] The simulation module is used to simulate the molecular dynamics equilibrium of the diffusion behavior of nanoparticles in the interstitial space under constant temperature conditions in a Langevin thermostat under a canonical ensemble, and to record the trajectory data of the nanoparticles.

[0102] The diffusion coefficient acquisition module is used to obtain the mean square displacement (MSD) of nanoparticles based on motion trajectory data. The MSD of nanoparticles is plotted against the simulation time t to obtain the MSD-t curve. The MSD-t curve is linearly fitted to solve for the slope of the MSD-t curve, and the diffusion coefficient of nanoparticles in the tissue interstitial space is determined based on the slope of the MSD-t curve.

[0103] The modules in the simulation system for the diffusion behavior of nanoparticles in interstitial spaces can be implemented entirely or partially through software, hardware, or a combination thereof. These modules can be embedded in or independent of the processor in a computer device, or stored in the memory of a computer device as software, so that the processor can call and execute the corresponding operations of each module.

[0104] The present invention also provides a computer device, including a memory, a processor, and a computer program stored in the memory. The processor executes the computer program to implement the steps in the method embodiment for simulating the diffusion behavior of nanoparticles in the interstitial spaces. Specific implementation methods can be found in the method embodiment, and will not be repeated here.

[0105] Furthermore, the present invention also provides a non-transitory computer-readable storage medium containing instructions on which a computer program is stored. For example, a memory containing instructions that can be executed by a processor of a computer device to perform the above-described method. For example, the non-transitory computer-readable storage medium may be a ROM, random access memory (RAM), CD-ROM, magnetic tape, floppy disk, and optical data storage device, etc. When the computer program is executed by the processor, it can implement the steps in the method embodiment simulating the diffusion behavior of nanoparticles in the interstitial spaces. Specific implementation methods can be found in the method embodiments, which will not be repeated here.

[0106] Those skilled in the art will understand that embodiments of the present invention can provide methods, systems, or computer program products. Therefore, the present invention can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present invention can take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0107] This invention is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and / or block diagrams, as well as combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart. Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0108] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0109] These computer program instructions can also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0110] It should be noted that the specific embodiments described above enable those skilled in the art to more fully understand the present invention, but do not limit the present invention in any way. Therefore, although the present invention has been described in detail in this specification and embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the present invention; and all technical solutions and improvements that do not depart from the spirit and scope of the present invention are covered within the protection scope of the present invention patent. No reference numerals in the claims should be construed as limiting the scope of the claims. Any simple variations or equivalent substitutions of technical solutions that can be readily obtained by those skilled in the art within the scope of the technology disclosed in the present invention are within the protection scope of the present invention.

Claims

1. A method for simulating the diffusion behavior of nanoparticles in tissue interstitial spaces, characterized in that, include: A multi-component coarse-grained molecular model was constructed, comprising a polymer network simulating the extracellular matrix (ECM), nanoparticles serving as drug carriers, and polymer chains embedded in the pores of the polymer network simulating macromolecular components in interstitial fluid. The energy of the multi-component coarse-grained model is minimized using the conjugate gradient method, and an energy-stable initial configuration is obtained by setting a convergence tolerance. Under a canonical ensemble, a Langevin thermostat was used to maintain constant temperature conditions. Molecular dynamics equilibrium simulations of the diffusion behavior of nanoparticles in the interstitial space were performed on the initial configuration, and the motion trajectory data of the nanoparticles were recorded. The mean square displacement (MSD) of the nanoparticles is obtained based on the motion trajectory data. The MSD of the nanoparticles is plotted against the simulation time t to obtain the MSD-t curve. The MSD-t curve is linearly fitted to solve for the slope of the MSD-t curve. The diffusion coefficient of the nanoparticles in the interstitial space is determined based on the slope of the MSD-t curve.

2. The method for simulating the diffusion behavior of nanoparticles in interstitial tissues according to claim 1, characterized in that, All components in the multi-component coarse-grained molecular model are composed of coarse-grained CG beads; The polymer network is an expandable topology structure formed by the cross-linking of multiple polymer chains, and its network size and porosity characteristics can be controlled by adjusting the lattice parameters; By adjusting the ratio of the total number of polymer chain CG beads to the total number of network CG beads within the range of 0 to 2, and adjusting the ratio of single chain length to network pore size within the range of 0 to 2, the interstitial environment of different viscosities from physiological to pathological states can be simulated. The nanoparticles are rigid structures composed of several CG beads linked by strong harmonic bonds, and their size and shape are configured according to the simulation requirements; the number of nanoparticles inserted in the multi-component coarse-grained molecular model is adjusted by adjusting the ratio of nanoparticles to the total volume of the system.

3. The method for simulating the diffusion behavior of nanoparticles in interstitial tissues according to claim 1, characterized in that, After constructing a multi-component coarse-grained molecular model, the process also includes setting a corresponding molecular dynamics interaction field; this field includes both bonded and unbonded interactions; the bonded interactions are described using harmonic potentials, including: The bond stretching interaction has the following potential energy function: In the formula, The potential energy is the bond length stretching energy. Let b be the bond stretching constant acting between particle i and particle j, and let b represent the bond stretching. The instantaneous distance between the particle pairs is [missing information]. To balance its bond length; The potential energy function of the bond-angle bending interaction is: In the formula, The potential energy is the bond angle bending energy. Let $\frac{i}{j}$ be the bending force constant acting on the bond angle formed by the sequentially connected particles $i$, $j$, and $k$. Indicates that the bond angle is bent. This is the instantaneous value of the bond angle. Its equilibrium value; The nonbonded interactions are described using the Lennard-Jones potential, whose potential function is: In the formula, For Lennard-Jones potential energy, and These are the interaction strength parameter between particle i and particle j and the characteristic distance when the potential energy is zero, respectively. The interaction force field sets different cutoff radii for different particle pairs, specifically including two types: Type I cutoff radius: its value is greater than This is used to preserve the full range of attractive forces in the Lennard-Jones potential in order to simulate long-range tunable interactions between particles; The second type of cutoff radius: its value is precisely set at the minimum value of the Lennard-Jones potential energy, and the potential energy is shifted to make the potential value at that point zero, thereby forming a modified potential field that retains only short-range repulsion.

4. The method for simulating the diffusion behavior of nanoparticles in the interstitial space according to claim 1, characterized in that, When using the conjugate gradient method to minimize the energy of a multi-component coarse-grained model, and obtaining an energy-stable initial configuration by setting convergence tolerances, both the energy convergence tolerance and the force convergence tolerance for energy minimization are set to 1e. -25 The maximum number of iterations required to minimize energy is 1,000,000. During the energy minimization process, the potential energy of the system, the size of the simulation box, the pressure tensor, and the total atomic potential energy are output every 10,000 steps.

5. The method for simulating the diffusion behavior of nanoparticles in interstitial tissues according to claim 1, characterized in that, Under a canonical ensemble, a Langevin thermostat was used to maintain constant temperature conditions. Molecular dynamics equilibrium simulations of the diffusion behavior of nanoparticles in the interstitial tissue were performed on the initial configuration, and the trajectory data of the nanoparticles were recorded. Specifically: The multi-component coarse-grained model after energy minimization was subjected to equilibrium simulation under a canonical ensemble; a Langevin heat bath was used to maintain an isothermal environment at a reduction temperature of 1.0ε; the integration step size was set to 0.002τ, and the total simulation time was 10,000,000 steps; the thermodynamic parameters and motion trajectory data of the system were monitored and output.

6. The method for simulating the diffusion behavior of nanoparticles in interstitial tissues according to claim 5, characterized in that, The system's thermodynamic parameters include energy parameters, temperature parameters, pressure parameters, and system size parameters; the motion trajectory data includes the time-series coordinates of the particles, the bonding and topological relationships of the particles, and the specific dynamic information of the nanoparticles.

7. The method for simulating the diffusion behavior of nanoparticles in interstitial tissues according to claim 1, characterized in that, The diffusion coefficient of nanoparticles in the interstitial space was calculated based on Einstein's diffusion equation. D : ;in t represents the slope of the MSD-t curve.

8. A simulation system for the diffusion behavior of nanoparticles in interstitial tissues, characterized in that, include: The model building module is used to build a multi-component coarse-grained molecular model, which includes a polymer network simulating the extracellular matrix (ECM), nanoparticles as drug carriers, and polymer chains embedded in the pores of the polymer network to simulate macromolecular components in interstitial fluid. The model optimization module is used to minimize the energy of a multi-component coarse-grained model using the conjugate gradient method, and obtains an energy-stable initial configuration by setting the convergence tolerance. The simulation module is used to perform molecular dynamics equilibrium simulation of the diffusion behavior of nanoparticles in the interstitial space under constant temperature conditions using a Langevin heat bath under a canonical ensemble, and to record the motion trajectory data of the nanoparticles. The diffusion coefficient acquisition module is used to obtain the mean square displacement (MSD) of the nanoparticles based on the motion trajectory data, plot the MSD of the nanoparticles against the simulation time t to obtain the MSD-t curve, perform linear fitting on the MSD-t curve, solve for the slope of the MSD-t curve, and determine the diffusion coefficient of the nanoparticles in the interstitial space based on the slope of the MSD-t curve.

9. A computer device, comprising a memory, a processor, and a computer program stored in the memory, characterized in that, The processor executes the computer program to implement the steps of the method according to any one of claims 1 to 7.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is loaded by the processor, it is able to perform the steps of the method according to any one of claims 1 to 7.

Images