A simulation method for the evolution of cytoskeleton porosity

Abstract

This invention discloses a simulation method for cytoskeleton porosity evolution, belonging to the field of computational biomechanics and micromanipulation simulation technology. The method first constructs a multi-component DPD model including the cell membrane and cytoskeleton; second, it introduces a dual-pathway Bell-Evans dynamics model, accurately characterizing the nonlinear lifetime characteristics of cross-linked proteins in the cytoskeleton under stress—characterized by "first strengthening, then weakening"—by coupling "capture" and "slip" paths; third, it constructs a stress redistribution factor based on atomic virial stress theory to quantify the risk of local stress transfer and cascade failure caused by bond breakage; finally, it employs an Euler-Lagrange coupled background mesh mapping method to solve the three-dimensional porosity field evolution during cytoskeleton deformation. This invention can faithfully reproduce the nonlinear mechanical response of cells under micromanipulation, revealing the cross-scale damage mechanism from microscopic bond breakage to macroscopic structural densification, providing a theoretical basis for path planning and damage assessment in cell microsurgery.

Description

Technical Field

[0001] This invention belongs to the fields of computational biology, biomechanical modeling and micro / nano robotics micromanipulation technology, specifically involving a computational method for simulating the evolution of the internal structure of cells under mechanical load using dissipative particle dynamics (DPD). Background Technology

[0002] Micromanipulation is a core technique in life science research, widely used in somatic cell nuclear transfer (SCNT), intracytoplasmic sperm injection (ICSI), transgenic microinjection, and single-cell drug delivery. With the development of robotics, automated micromanipulation has become a trend, significantly improving throughput and reproducibility. However, regardless of whether it's manual or robotic manipulation, cells undergo significant deformation during physical contact due to the mechanical forces applied by microneedles, optical tweezers, or microfluidic chips.

[0003] Mechanical damage to cells is a key factor limiting the success rate of micromanipulation. Excessive mechanical stress can lead to cell membrane rupture, cytoskeleton depolymerization, or irreversible breakage, thereby inducing apoptosis or necrosis. Therefore, a deep understanding of the mechanical response and damage mechanisms of cells under complex mechanical loads is of great significance for optimizing manipulation parameters and designing novel micro / nano tools.

[0004] To predict the mechanical behavior of cells, researchers have developed a variety of computational models, mainly including continuum mechanics models (such as the finite element method, FEM) and discrete particle models (such as molecular dynamics, MD, and dissipative particle dynamics, DPD).

[0005] Limitations of the continuum model: The FEM model treats the cell as a homogeneous or layered viscoelastic continuum. While it is efficient in calculating macroscopic deformations, it cannot describe the discrete microscopic structure of the cytoskeleton, nor can it simulate discontinuous processes such as chemical bond breaking and molecular rearrangement. For problems involving large deformations and topological changes such as cell puncture and crack propagation, continuum mechanics faces difficulties due to mesh distortion and constitutive equation failure.

[0006] Limitations of traditional discrete particle models: While DPD and coarse-grained MD methods can effectively capture cell membrane fluidity and cytoskeleton network features, existing DPD cell models suffer from significant oversimplification in describing intermolecular interactions, primarily in the following three aspects: 1. Simplified "Slip Bond" Breaking Mechanism: Existing models typically use simple "critical length" or "constant force" as bond breakage criteria. These models assume that the greater the external force, the easier the bond breaks, a phenomenon known as "slip bond" behavior. However, numerous biophysical experiments have confirmed that cell adhesion molecules (such as selectins) and skeletal crosslinking proteins (such as actin-crosslinkers) exhibit counterintuitive "catch bond" behavior: within a certain range of tension, external force can induce conformational changes in proteins, resulting in tighter binding and extended bond lifespan. Ignoring the catch bond mechanism leads to models that severely underestimate the mechanical stability of cells under low and moderate stress, failing to explain the cell's "strengthened by increased stress" mechanical adaptation.

[0007] 2. Lack of quantitative analysis of stress redistribution: In discrete networks, the fracture of a fiber causes the load it bears to be released instantaneously and transferred to neighboring fibers; this process is called "mechanical redistribution." If neighboring fibers cannot withstand the new load, cascading fracture will occur. Existing techniques mostly focus on macroscopic strain contour maps, lacking the tracking and quantification of microscopic stress transfer paths, making it difficult to identify potential damage nucleation sites.

[0008] 3. Neglecting the dynamic evolution of cytoskeleton porosity: The cytoskeleton is a porous medium filled with cytoplasmic fluid. During micromanipulation (such as squeezing and suction), the cytoskeleton network deforms, leading to drastic changes in local porosity. Changes in porosity directly affect the diffusion and transport of intracellular substances (such as drug molecules and organelles). Existing models often simplify the cytoplasm to background friction, failing to establish the connection between mechanical deformation and changes in the structure of material transport channels.

[0009] In summary, existing technologies lack a high-fidelity cell simulation method that can simultaneously integrate nonlinear bonding dynamics, micro-stress redistribution mechanisms, and porous media structure evolution. Summary of the Invention

[0010] This invention aims to overcome the aforementioned deficiencies of existing technologies and provide a simulation method for cytoskeleton porosity evolution based on the Catch Bond dynamic fracture mechanism and mechanical redistribution. This invention guides dissipative particle dynamics simulations by constructing a sophisticated physical-mathematical model, thereby obtaining analytical results rich in deep biomechanical insights.

[0011] The methodology of this invention comprises three core innovative modules: a Catch-Slip dual-path fracture model based on Kramers theory, a stress redistribution calculation module based on the virial theorem, and a dynamic porosity calculation module based on Euler-Lagrange coupling.

[0012] The specific technical process is as follows: Step 1: Constructing a multi-component heterogeneous cell microscopic model Within the framework of dissipative particle dynamics (DPD), a three-dimensional discrete model is constructed, comprising cell membrane particles, cytoskeleton particles, fluid environment (as a cytoplasmic equivalent), and cross-linked protein particles.

[0013] Cell membrane particles: A triangular mesh network topology is adopted, and the interaction between particles is represented by the Worm-Like Chain (WLC) potential function. The WLC model can accurately describe the entropic elastic behavior of polymer chains under stretching and the strain stiffening characteristics under limiting elongation, which is crucial for simulating the real mechanical response of cell membranes under large deformations caused by microneedle insertion.

[0014] Cytoskeleton particles are modeled as a network composed of randomly cross-linked actin microfilaments and cross-linked proteins. The bonds between cytoskeleton particles are modeled using a harmonic potential function, with the initial configuration optimized through short-time energy minimization or kinetic relaxation. Harmonic potential functions offer high computational efficiency and highly adjustable parameters, making them suitable as the fundamental mechanical framework for complex catch bond breakage mechanisms.

[0015] Fluid environment: Represented by the implicit background friction field, used to provide equivalent intracellular fluid damping and fluid dynamic interactions.

[0016] An initial skeleton network is constructed using a constrained random fiber generation algorithm, and local overlap and residual stress are eliminated by short-time energy minimization or dynamic relaxation.

[0017] Step 2: Define the dynamic breakage mechanism of Catch Bond based on the Bell-Evans model In order to simulate the real mechanochemical behavior of proteins, this invention abandons the deterministic breakage criterion and adopts the classic Bell-Evans theory to model bond lifetime.

[0018] Introducing a modified Bell-Evans two-pathway model to describe bond dissociation rates Bond tension The model includes a Catch component that decreases with force and a Slip component that increases with force; according to the relationship between the changes, the model contains a Catch component that decreases with force and a Slip component that increases with force; Calculate the simulation time step Bond breakage probability within And perform probabilistic key breaking operations in the simulation.

[0019] The Bell-Evans dual-pathway model posits that bond dissociation is a superposition of two independent physical processes. The catch pathway describes the physical process in which external forces induce conformational changes in proteins within a low-stress range, thus hindering bond breaking; while the slip pathway describes the process in which external forces overcome energy barriers within a high-stress range, thus accelerating bond breaking.

[0020] Dissociation rate equation: ; in, The tensile scalar value borne by the key; and These are the zero-force dissociation rate and characteristic force parameters of the Catch path, respectively; therefore, the Catch component... It decreases with increasing force, representing the strengthening stage of the bond; and These represent the zero-force dissociation rate and characteristic force parameters of the Slip path, respectively; therefore, the Slip component... It increases with increasing force, representing the weakening stage of the bond.

[0021] This formula accurately reproduces the non-monotonic behavior of bond lifetime, which first increases with force (Catch-dominated) and then decreases with force (Slip-dominated), as described in literature (e.g., Weak catch bonds make strong networks).

[0022] Random breakage algorithm: based on Calculate the simulation time step Bond breakage probability within The fracture probability The calculation formula is: The Monte Carlo method is used to determine whether a bond breaks: at each step of the simulation, random numbers are generated. ,like If the bond is broken, a probabilistic bond breakage operation will be performed in the simulation.

[0023] Step 3: Construction and dynamic simulation of the micromanipulation simulation environment Introduce parameterized microneedle or indenter models into the simulation space and set their geometric parameters (radius, cone angle) and kinematic parameters (velocity, acceleration).

[0024] Set boundary conditions: Use Lennard-Jones potential energy walls to simulate the constraint effect of the substrate or holding pin.

[0025] Perform micromanipulation simulations, start a DPD integrator (such as Velocity-Verlet), and perform time evolution. During the simulation, update particle positions, velocities, and bond topology in real time, and record the entire particle trajectory.

[0026] Step 4: Calculate the virial stress and mechanical redistribution. Based on the microscopic particle dynamics data output from the simulation, the instantaneous local stress tensor of each skeletal particle is calculated using the Virial Theorem; within the neighborhood of the bond fracture event, the stress redistribution factor (SRF) is calculated to quantify the stress transfer and concentration caused by the fracture event.

[0027] Virial stress calculation: Calculating the atomic stress tensor of each particle using the virial theorem. It includes both kinetic and potential energy contributions. Local virial stress tensor The calculation formula is: ; in, For particles The representative volume; For particles The quality; and Particles In Cartesian coordinates , Velocity component in the direction; For particles With particles The relative position vector in Components in direction; For particles Acting on particles The interaction forces in Components in direction; subscript , representing the direction of the Cartesian coordinate system .

[0028] Redistribution factor (SRF): Tracks the bonds that break and calculates the rate of change of von Mises stress in neighboring particles before and after the breakage.

[0029] The formula for calculating the stress redistribution factor (SRF) is as follows: ; in, Represents neighborhood particles Stress redistribution factor; Represents particles At the moment of fracture The von Mises equivalent effect; Represents particles After the fracture occurs, a relaxation time window is passed. The von Mises equivalent effect at that time; The moment of fracture; This refers to the relaxation time window after fracture. Number the particles adjacent to the broken bond.

[0030] high The value region indicates that the stress of the particles in the vicinity increases significantly after fracture, reflecting that stress is concentrated in a local area, which is a potential area of ​​secondary fracture risk.

[0031] Step 5: Calculation of skeleton porosity evolution Background grid mapping is employed to divide the simulation space into regular Eulerian grids. A kernel function maps the volume of skeletal particles to grid nodes, calculating the volume ratio of each grid cell. The evolution of the three-dimensional porosity field during skeletal deformation is then calculated in real-time or based on trajectory post-processing, yielding the three-dimensional spatial distribution of intracellular porosity and its evolution over time. Specifically, this includes: 5.1 Spatial Discretization: Establish a regular 3D Eulerian background mesh covering the bounding box of the cell model within the simulation area, with a mesh step size of [missing information]. Set as average particle spacing 0.5-1.0 times; 5.2 Volume Projection: Defining the Kernel Function (such as Gaussian kernels or spline kernels), to structure particles volume Grid nodes assigned to its influence domain Above, compute grid nodes m Solid volume fraction at the point : ; in, For grid nodes Spatial position vector, skeletal particles Spatial position vector, Represents grid nodes With skeletal particles The Euclidean distance between them skeletal particles The representative volume, This represents the set of skeleton particles located within the support domain of the kernel function. Represents grid nodes The volume fraction of solids at that location.

[0032] 5.3 Calculating the porosity of mesh elements: The local porosity at mesh element m is... ,in Indicates the volume fraction of solids. The smaller of the two values, 0 and 1, limits the local porosity to the range of 0 to 1.

[0033] 5.4 For all meshes Interpolation and smoothing processes are performed to generate a porosity scalar field. By analyzing the three-dimensional cloud map or cross-sectional heat map of porosity changes over time, the evolution characteristics of the compacted zone (densification) and the loose zone (expansion) are identified.

[0034] Step 6: Multiphysics Coupling Assessment This study couples and evaluates the local damage evolution and overall mechanical state of cells during micromanipulation by integrating the stress field, damage field (bond breakage distribution), porosity field, and risk dissociation rate extracted based on bond breakage event sequences. By correlating bond breakage locations, high virial stress concentration zones, and porosity abrupt change zones, the influence of manipulation parameters such as microneedle insertion velocity and needle tip shape on damage to intracellular structures is analyzed. Outputs include: virial stress contour maps, bond breakage distribution maps, porosity evolution heatmaps or curves, risk dissociation rate curves, and time-series variation curves of the global damage index and stress concentration coefficient, thus generating a comprehensive mechanical analysis report.

[0035] Step 1, in constructing the cytoskeleton model, abandons the traditional method of random point scattering combined with Monte Carlo relaxation, and instead adopts a geometrically constrained, restricted random fiber generation algorithm. The construction process specifically includes: first, randomly selecting a starting coordinate point within the cell membrane space and randomly setting a growth direction vector; then generating a virtual line segment based on a preset fiber length; subsequently, performing geometric constraint checks on this virtual line segment. The passing criteria are: all sampling points of the virtual line segment are located inside the cell membrane envelope, and the minimum distance between the virtual line segment and the generated cytoskeleton fibers is not less than a preset safety distance; the failing criteria are: any sampling point of the virtual line segment extends beyond the cell membrane boundary, or spatially intersects or overlaps with the generated cytoskeleton fibers, or the minimum distance between them is less than the preset safety distance. If the detection fails, the line segment is discarded and the test is repeated. If the detection passes, the line segment is retained, and cytoskeleton particles are generated uniformly along its direction at a fixed discretization step size. Finally, chemical bonds are established between adjacent cytoskeleton particles, and their potential function is set as a harmonic potential energy. Harmonic angular potential energy is applied between three consecutive cytoskeleton particles on the same straight line to maintain fiber rigidity. The above steps are repeated until the total number of cytoskeleton particles reaches a preset threshold.

[0036] In step 3, the micromanipulation simulation supports both force loading and displacement loading modes; in displacement loading mode, the microneedle operates at a constant speed. Press in until the pressing depth reaches the set value. Then hold or pull back; during the hold phase, count the delayed fracture phenomenon under the Catch Bond mechanism.

[0037] Step 3 requires the force value at the Catch-Slip inflection point. To determine the safe range of applied forces during micromanipulation, the specific method is as follows: Let The critical force value at which the bond lifetime is longest (most stable) is obtained.

[0038] The background mesh is activated only in the internal region of the cell membrane. By identifying the envelope of cell membrane particles, invalid meshes in the external fluid region are eliminated, thereby improving the accuracy of porosity calculation.

[0039] The risk dissociation rate is obtained by performing Savitzky-Golay local filtering and analytical differentiation on the cumulative number of broken bonds to obtain the macroscopic breakage rate, and then dividing the macroscopic breakage rate by the current number of remaining intact bonds. It is used to characterize the dissociation risk of the unit remaining bonds of the cross-linked protein under dynamic stress conditions.

[0040] The simulation software is a secondary development based on LAMMPS (Large-scale Atomic / Molecular Massively Parallel Simulator). It implements the probabilistic bond breaking logic in step 2 and the real-time output of virial stress in step 4 by writing user-defined fix programs.

[0041] The method of the present invention also includes the calculation of a global evaluation index for the degree of cell damage: defining a global damage index. The ratio of the current number of fractured bonds to the initial number of bonds; the stress concentration factor is defined. The ratio of the maximum von Mises stress to the mean stress within the cell is used to calculate the macroscopic fracture rate based on the curve of the cumulative number of broken bonds over time. The risk dissociation rate per unit of remaining bond is calculated in combination with the current number of intact bonds. The global damage index curve, stress concentration coefficient curve, and risk dissociation rate curve are output to analyze the temporal correlation between damage evolution, stress concentration, and bond dissociation risk.

[0042] Compared with the prior art, the beneficial effects of this invention are as follows: To address the shortcomings of existing dissipative particle dynamics cell models, such as a simplistic bond breakage mechanism, a lack of quantification methods for stress redistribution, and an inability to describe the evolution of internal pore structures, this invention proposes a multiphysics coupling simulation scheme. This invention specifically focuses on the nonlinear breakage behavior of protein bonds in the cytoskeleton network, the microscopic stress redistribution caused by breakage, and the spatiotemporal dynamic changes of intracellular micropore structures. It is applicable to scenarios such as cell microinjection, mechanical property measurement, and single-cell damage assessment. Attached Figure Description

[0043] Figure 1 This is a flowchart illustrating the overall calculation process of the method of the present invention; Figure 2 This is a schematic diagram comparing the dissociation rates of the Catch Bond and Slip Bond of the present invention with force. Figure 3 This is a schematic diagram illustrating the porosity calculation principle based on the background grid method of the present invention; Figure 4 This is a schematic diagram showing the distribution of locally high-stress bond segments inside the cytoskeleton during the pressing process of the indenter of this invention; Figure 5 This is a heatmap showing the porosity distribution at the center cross-section of the cell in this invention. Figure 6 This is a comparison curve of the unit residual bond risk dissociation rate between the Catch Bond model and the Slip Bond model of this invention under the same loading conditions. Detailed Implementation

[0044] Exemplary embodiments of the present invention will now be described in more detail with reference to the accompanying drawings. While exemplary embodiments of the present invention are shown in the drawings, it should be understood that the present invention can be implemented in various forms and should not be limited to the embodiments set forth herein. Rather, these embodiments are provided to enable a more thorough understanding of the present invention and to fully convey the scope of the invention to those skilled in the art. It should be noted that, unless otherwise specified, the embodiments and features described herein can be combined with each other. The present invention will now be described in detail with reference to the accompanying drawings and embodiments.

[0045] See Figures 1 to 6 The embodiments of the present invention provide a method for calculating the internal strain of a multi-component cell model based on dissipative particle dynamics (DPD). The method uses Python scripts for parametric modeling and data post-processing, and LAMMPS (Large-scale Atomic / Molecular Massively Parallel Simulator) for dynamic simulation and implementation of the Catch Bond mechanism. Figure 1 It demonstrates the entire process from model building to post-processing analysis, specifically including the following steps: Step 1: Use a Python script to build a multi-component cell geometry model and generate data files. Unlike traditional methods that rely on MATLAB toolboxes, this embodiment uses Python's open-source geometry processing library for modeling to improve model generation efficiency and topology quality.

[0046] Write a Python script to use the trimesh and numpy libraries to generate the initial configuration file (.data format) for the cell model: Cell membrane particle modeling: The radius is generated using the trimesh.creation.icosphere(subdivisions=3, radius=10.0) function. A spherical triangular mesh (simulation unit). Extract the vertices of the mesh as cell membrane particles (type 1), and extract the edges of the mesh as topological connections between membrane particles.

[0047] Cytoskeleton particle modeling: Within the cell membrane space, a geometrically constrained, restricted random fiber generation algorithm (Type 2) is used to generate cytoskeleton particles. Specifically, numpy.random is used to randomly assign the starting coordinates and growth direction of the fibers, and virtual line segments are formed according to the preset fiber length. The virtual line segments are then tested for in-membrane integrity, intersection / overlap, and minimum spacing. If the tests pass, cytoskeleton particles are generated along the line segments at a fixed spacing. If the tests fail, the particles are discarded and resampled.

[0048] File Output: The function follows the LAMMPS atom_style hybrid bond format to write the particle coordinates, bond connection information, and the angle and dihedral topology information derived from the mesh facets into the cell_model.data file. This file declares the number and type of particles, bonds, angles, and dihedrals.

[0049] Step 2: Set up the simulation environment and particle interaction potential in LAMMPS Write the LAMMPS input script.

[0050] Environment Setup: Use the command `units lj` to set the unit system. To simulate a microscopic manipulation environment, set the X-axis to a fixed boundary to prevent particle penetration; set the Y and Z axes to periodic boundaries to ensure lateral momentum and energy conservation. Use the command `read_data cell_model.data` to read the model file generated in step 1.

[0051] Potential function setting: Nonbonded interactions: Use the command `pair_style dpd 1.0 1.0 12345` to set the dissipative particle dynamic potential. Set the repulsive parameters (Conservative Force) between different components: membrane-membrane parameter is `pair_coeff 1 1 10045 0.5`; framework-framework parameter is `pair_coeff 2 2 100 65 0.5`; membrane-framework parameter is `pair_coeff 12 100 45 1.0` (simulating differences in hydrophilicity).

[0052] Bonding interactions: Use the command `bond_style hybrid wlc harmonic`.

[0053] For the bonds between cell membrane particles, a worm-like chain potential is used to simulate the nonlinear entropy elasticity of the cell membrane. The command is bond_coeff 1 wlc 1.0 1.0 0.5.

[0054] For the bonds between cytoskeleton particles, a simple harmonic potential is used as the basic carrier of the Catch Bond mechanism. The command is bond_coeff 2 harmonic 50.0 0.5.

[0055] Step 3: Implement skeleton aggregation and Bell-Evans dynamic fracture mechanism in LAMMPS Initial network construction: The restricted random fiber backbone generated in step 1 is imported into LAMMPS as the initial main chain network; during the simulation relaxation phase, fix bond / create or equivalent custom bonding logic is used to create cross-links only between candidate cross-linked protein-backbone particles that meet the distance threshold, type constraints and maximum bond number limit, thereby forming a cross-linked backbone network consistent with the restricted random fiber topology.

[0056] Introducing the Bell-Evans dual-pathway fracture model: To simulate realistic biomechanical responses, this invention abandons the single critical force fracture criterion and adopts the Bell-Evans dual-pathway model to describe bond dissociation rates. With bond tension Relationship:

[0057] in, and These are the zero-force dissociation rates for the Catch path and the Slip path, respectively. and These represent the characteristic force scales under two different paths. In LAMMPS, this mechanism is implemented by writing custom repair programs bond_bell_break.cpp and bond_bell_break.h: at each time step, the tension is calculated based on the current bond length. Substitute into the above formula to calculate the fracture probability. The mechanism enables the skeleton network to exhibit dynamic reorganization (high dissociation rate) under low stress and extended lifetime (Catch effect) under moderate stress, thus achieving intelligent stress redistribution.

[0058] Step 4: Set up the micromanipulation experimental environment and run the simulation. Loading tool: Define a spherical indenter on the right side of the simulation box. Use the command... Define the motion trajectory of the pressure head.

[0059] Fix the base: Set up a repulsive plate (Wall) on the left side of the cell (negative X-axis direction) and use the command fix wall / reflect xlo... to prevent the cell from translating as a whole.

[0060] Simulation run: Set timestep to 0.0005. Under the NVE (micro-regular) ensemble, use the fixindent command to drive the pressure head at a speed of... Insert the needle into the cell along the negative X-axis, pressing it to a depth equal to the cell radius. .

[0061] Data output: The dump custom command outputs a full particle trajectory file dump.lammpstrj every 1000 steps. The output includes: particle ID, type, and three-dimensional coordinates (x, y, z).

[0062] Step 5: Multiphysics Post-processing – Strain Field and Porosity Field Calculation After the simulation, a strategy combining "Ovito visualization analysis" and "Python deep computing" was adopted to calculate the mechanical strain and microporosity evolution of the cells.

[0063] 5.1 Finite Strain Calculations Using Ovito Data import and configuration: Open Ovito software and load the dump.lammpstrj file. Add the Atomic Strain modifier to the modifier list.

[0064] Parameter settings: Set the cutoff radius to 2.0 to define the computational neighborhood of the local deformation gradient.

[0065] Select frame 0 in the Reference configuration.

[0066] Check the "Output strain tensors" option.

[0067] Visualizing the slice: Add a Slice modifier to cut the cell along the XZ plane, add color coding, and map the attribute to Strain Tensor.XX (normal strain along the insertion direction). At this point, the concentrated area of ​​compressive strain in front of the indenter can be visually observed, which helps to analyze the spatial distribution characteristics of local deformation and damage concentration during micromanipulation.

[0068] 5.2 Calculating the dynamic porosity field using Python scripts Figure 3 This demonstrates the mapping process from the volume of skeletal particles to background mesh nodes. Specifically, on a regular 3D Eulerian background mesh, using skeletal particles as Lagrangian discrete points, their representative volumes are distributed to neighboring mesh nodes using a kernel function with distance weights; mesh nodes closer to the particles receive larger volume allocations, while nodes farther away receive smaller allocations. This allows the construction of a continuous solid-phase volume fraction field on the Eulerian background mesh, and further calculation of local porosity distribution.

[0069] To quantify the impact of skeletal deformation on intracellular transport channels during micromanipulation, a Python script was written to calculate porosity distribution using the Eulerian-Lagrange coupled background grid mapping method. The specific process is as follows: 5.2.1 Mesh Construction and Spatial Discretization: A regular 3D Eulerian mesh is constructed within the simulation region containing the maximum bounding box of cell deformation using numpy.meshgrid. The mesh step size is then adjusted. The average spacing between skeletal particles is set to 0.5 to 1.0 times (in this embodiment). Simulation units are used to ensure sufficient spatial resolution.

[0070] 5.2.2 Boundary Recognition and Mask Generation: To avoid interference from the extracellular fluid region on porosity calculations, a closed envelope surface is first constructed based on the cell membrane particle coordinates (e.g., using the Convex Hull algorithm or the Alpha Shape algorithm). A binary mask matrix is ​​then defined. For any grid node If it is located inside the cell envelope, then it is marked. Otherwise, mark .

[0071] 5.2.3 Kernel Mapping: For each skeletal particle (Location ,volume The volume is distributed to the mesh nodes within its supporting domain using a cubic spline kernel function in Smooth Particle Hydrodynamics (SPH). Above. Solid volume fraction at grid nodes. The calculation formula is:

[0072] in, For smooth length, kernel function This ensures volume conservation and spatial continuity. Nearest neighbor search is performed using scipy.spatial.cKDTree to accelerate the mapping process.

[0073] 5.2.4 Calculation of local porosity: Define local porosity This represents the proportion of fluid within a unit volume. The calculation formula is:

[0074] This step converts the solid volume fraction into a porosity field and restricts the value range to the interval [0,1].

[0075] 6. Results Output and Visualization: Contour plotting: Use matplotlib.pyplot.contourf to plot the cross-sectional heat map of porosity. In the plotting settings, force the background color to pure white (RGB:255,255,255), remove the default gray background and grid lines, and use the "Viridis" color level to display dense areas with high contrast.

[0076] Results Analysis: The analysis results show that the skeleton is severely compressed around the indenter's insertion path and in front of the needle tip, resulting in a significant decrease in porosity. This forms a dense region that hinders the diffusion of matter; this is highly consistent with the spatial distribution of the hardening phenomenon in the high-stress region caused by the Catch Bond mechanism, verifying the coupling relationship between mechanical response and structural evolution.

[0077] 6.1 Macroscopic Fracture Dynamics and Risk Dissociation Rate Extraction Based on Bond Breaking Event Sequence This section pertains to the statistical post-processing of the probabilistic bond breaking results in step 2: step 2 is responsible for calculating the breakage probability and performing bond breaking within the simulation time step, while this section extracts the macroscopic breakage rate and hazard dissociation rate based on the recorded bond breaking time series. In microscopic molecular dynamics (MD) simulations, the breaking of cross-linked bonds essentially follows a discrete Poisson-like stochastic process. To eliminate numerical noise caused by high-frequency thermal fluctuations in the system and effectively remove the interference of residual network stress on the dynamic evolution, this invention designs a post-processing algorithm based on Savitzky-Golay (SG) local filtering and hazard dissociation rate calculation to accurately extract the true fracture dynamics characteristics under macroscopic mechanical response. The specific steps are as follows: 6.1.1 Simulation Time-Domain Mapping and Cardinality Statistics Extract the cumulative number of breakages of the slip bond and catch bond at each time step from the simulation terminal. The total number of initial keys recorded based on the initial network topology. Calculate at any time Number of remaining cross-linked bonds in good condition:

[0078] 6.1.2 Dynamic Smoothing and Analytical Differentiation Based on Savitzky-Golay Algorithm Because directly performing finite-difference differentiation on discrete cumulative data produces a severe "noise amplification effect," and conventional global polynomial fitting is prone to inducing Runge oscillations at the boundaries, this step introduces a large-window Savitzky-Golay (SG) filter. This algorithm uses a low-order polynomial (such as second-order) for local least-squares fitting within the sliding data window and directly outputs the smoothed analytical first derivative, thereby obtaining the instantaneous macroscopic fracture rate.

[0079] This operation forcibly flattens out the random thermal fluctuations at the microscale, extracting the fracture rate trend that is entirely driven by macroscopic compressive strain.

[0080] 6.1.3 Hazard Rate Reconstruction In the later stages of loading, the macroscopic fracture rate decreases due to the sharp reduction in the number of intact cross-linked bonds. A spurious decrease may occur due to the "depletion effect." This is to restore the probability of single-bond dissociation in the Bell-Evans theoretical model. This step calculates the risk dissociation rate of the unit remaining bond. :

[0081] This index precisely maps the true microscopic probability of cross-linked proteins dissociating in a dynamic time-varying stress field.

[0082] 6.1.4 Initial Relaxation Shielding and Dynamic Peak Extraction In the initial stage of the simulation, the presence of local residual "pre-strain" in the randomly cross-linked network easily leads to early and extensive breakage of weak bonds, forming "pseudo-peaks" not driven by external loads. To accurately quantify the mechanical failure dominated by external extrusion, a time threshold is introduced in this step. To shield the initial energy relaxation region, and in Within the forced deformation range, extract the maximum peak value of the risk dissociation rate. and their corresponding peak times :

[0083] 6.2 Comparison of Dynamic Characteristics and Analysis of Results Figure 4 This is a schematic diagram of the distribution of locally high-stress bond segments inside the cytoskeleton during the indentation process. The red area represents the high-stress concentration zone. Figure 5 Heatmap of porosity distribution in the central section of a cell; heatmap of porosity evolution in the central section of a cell. Figure 6This paper compares the unit residual bond risk dissociation rates of the Catch Bond model and the Slip Bond model under the same loading conditions. The results show that, facing the same initial network prestress and macroscopic extrusion deformation, the Slip bond exhibits an extremely high early dissociation rate (fragility characteristic); while the Catch bond shows significant dissociation rate stagnation and a downward curve in the early stages of extrusion. Quantitative data confirm that the catastrophic yield peak time of the Catch bond (…) Lags behind the Slip key ( This dynamic evolution curve provides conclusive evidence from a statistical physics perspective: Catch Bond's "force-strengthening" inverse locking mechanism successfully suppressed initial network damage and significantly improved the topological stability and toughness of the cytoskeleton when subjected to large deformation mechanical damage.

[0084] The present invention has been described in detail above through embodiments, but the content described is only an exemplary embodiment of the present invention and should not be considered as limiting the scope of the present invention. The scope of protection of the present invention is defined by the claims. Any technical solutions designed by those skilled in the art using the technical solutions described in the present invention, or similar technical solutions designed by those skilled in the art under the inspiration of the technical solutions of the present invention, within the substance and scope of protection of the present invention, to achieve the above-mentioned technical effects, or equivalent changes and improvements made to the scope of the application, should still fall within the patent protection scope of the present invention. It should be noted that, for clarity, descriptions of some components and processes that are not directly and obviously related to the scope of protection of the present invention but are known to those skilled in the art have been omitted in the description of the present invention.

Claims

1. A method for simulating the evolution of cytoskeleton porosity, characterized in that, Includes the following steps: Step 1: Within the framework of dissipative particle dynamics, a three-dimensional cell model is established, including cell membrane particles, cytoskeleton particles, fluid environment, and cross-linked protein particles. The cell membrane particles adopt a triangular mesh network topology, and the inter-particle connections adopt a worm-like chain potential energy model. The cytoskeleton particles adopt a randomly cross-linked network topology, and the inter-particle connections adopt a harmonic potential energy model. The initial cytoskeleton network is constructed using a constrained random fiber generation algorithm, and local overlap and residual stress are eliminated through short-time energy minimization or dynamic relaxation. Step 2: Introduce a modified Bell-Evans two-pathway model to describe the bond dissociation rate. Bond tension The model includes a Catch component that decreases with force and a Slip component that increases with force; according to the relationship between the changes, the model contains a Catch component that decreases with force and a Slip component that increases with force; Calculate the simulation time step Bond breakage probability within And perform probabilistic key breaking operations in the simulation; Step 3: Load the microneedle or indenter module into the simulation environment, set the parameters and boundary conditions, and perform micromanipulation simulation; during the simulation, update the particle position, velocity and bond connection topology in real time, and record the entire particle trajectory; Step 4: Based on microscopic particle dynamics data, calculate the instantaneous local stress tensor of each skeletal particle using the virial theorem; calculate the stress redistribution factor in the neighborhood where the bond breakage event occurs, and quantify the stress transfer and concentration caused by the breakage event. Step 5: Using the background mesh mapping method, the simulation space is divided into regular Eulerian meshes. The volume of the skeletal particles is mapped to the mesh nodes through the kernel function. The volume ratio of each mesh unit is calculated, thereby obtaining the three-dimensional spatial distribution of the porosity inside the cell and its evolution cloud map over time. Step 6: Correlate the bond fracture location, high-dimensional stress zone, and porosity abrupt change zone, and extract the macroscopic fracture rate and risk dissociation rate based on the fracture bond time series, outputting a comprehensive mechanical analysis report including stress cloud map, damage distribution map, porosity evolution curve, and risk dissociation rate curve.

2. The method according to claim 1, characterized in that, The construction process in step 1 is as follows: First, a starting coordinate point is randomly selected in the internal space of the cell membrane, and a growth direction vector is randomly set. A virtual line segment is generated according to the preset fiber length. Then, the virtual line segment is subjected to geometric constraint detection. The passing criteria for the geometric constraint detection are: all sampling points of the virtual line segment are located inside the cell membrane envelope, and the minimum distance between the virtual line segment and the generated skeleton fiber is not less than the preset safety distance. The failing criteria are: any sampling point of the virtual line segment penetrates the cell membrane boundary, or spatially intersects or overlaps with the generated skeleton fiber, or the minimum distance between the two is less than the preset safety distance. If the detection fails, the line segment is discarded and retried. If the detection passes, the line segment is retained, and skeleton particles are generated uniformly along its direction at a fixed discretization step size. Finally, chemical bonds are established between adjacent skeleton particles, and their potential function is set as harmonic potential energy. Harmonic angular potential energy is applied between three consecutive skeleton particles on the same straight line to maintain fiber rigidity. The above steps are repeated until the total number of cytoskeleton particles reaches the preset threshold.

3. The method according to claim 1, characterized in that, In step 2, the bond dissociation rate The mathematical model adopts the two-path Bell-Evans model, and its expression is: ； in, The tensile scalar value borne by the key; and These represent the zero-force dissociation rate and characteristic force parameters of the Catch path, respectively. and These represent the zero-force dissociation rate and characteristic force parameters of the Slip path, respectively. The probability of fracture The calculation formula is: In each step of the simulation, random numbers are generated. If the bond is broken, a probabilistic bond breakage operation will be performed in the simulation.

4. The method according to claim 1, characterized in that, In step 3, the micromanipulation simulation supports both force loading and displacement loading modes; in displacement loading mode, the microneedle operates at a constant speed. Press in until the pressing depth reaches the set value. Then hold or pull back; during the hold phase, statistical analysis of delayed breakage phenomena under the Catch Bond mechanism is performed.

5. The method according to claim 3, characterized in that, Step 3 requires the force value at the Catch-Slip inflection point. To determine the safe range of applied forces during micromanipulation, the specific method is as follows: Let The critical force value at which the bond lifetime is longest is obtained.

6. The method according to claim 1, characterized in that, In step 4, the local virial stress tensor The calculation formula is: ； in, For particles The representative volume, For the mass of the particle, and Particles exist , velocity components in the direction, For particles With particles The relative position of Components in direction, For particles and The interaction force between them in the b direction; subscript Represents the direction of the Cartesian coordinate system ; The formula for calculating the stress redistribution factor is: ； in, Represents neighborhood particles Stress redistribution factor; Represents particles At the moment of fracture The von Mises equivalent effect; Represents particles After the fracture occurs, a relaxation time window is passed. The von Mises equivalent effect at that time; The moment of fracture; This refers to the relaxation time window after fracture. Number the particles adjacent to the broken bond.

7. The method according to claim 1, characterized in that, In step 5, the calculation of the dynamic porosity field is specifically... include: 5.1 Spatial Discretization: Establish a regular 3D Eulerian background mesh covering the bounding box of the cell model within the simulation area, with a mesh step size of [missing information]. Set as average particle spacing 0.5-1.0 times; 5.2 Volume Projection: Defining the Kernel Function , scaffold particles volume Grid nodes assigned to its influence domain Above, compute grid nodes m Solid volume fraction at the point : ； in, For grid nodes Spatial position vector, skeletal particles Spatial position vector, Represents grid nodes With skeletal particles The Euclidean distance between them skeletal particles The representative volume, This represents the set of skeleton particles located within the support domain of the kernel function. Represents grid nodes The volume fraction of solids at that location; 5.3 Calculating the porosity of mesh elements: The local porosity at mesh element m is... ; 5.4 For all meshes Interpolation and smoothing processes are performed to generate a porosity scalar field. By analyzing the three-dimensional cloud map or cross-sectional heat map of porosity changes over time, the evolution characteristics of compacted and loose zones can be identified.

8. The method according to claim 7, characterized in that, The background mesh is activated only in the internal region of the cell membrane. By identifying the envelope of cell membrane particles, invalid meshes in the external fluid region are eliminated, thereby improving the accuracy of porosity calculation.

9. The method according to claim 1, characterized in that, The risk dissociation rate is obtained by performing Savitzky-Golay local filtering and analytical differentiation on the cumulative number of broken bonds to obtain the macroscopic breakage rate, and then dividing the macroscopic breakage rate by the current number of remaining intact bonds. It is used to characterize the dissociation risk of the unit remaining bonds of the cross-linked protein under dynamic stress conditions.

10. The method according to claim 1, characterized in that, The method also includes the calculation of a global evaluation index for the degree of cell damage: the global damage index is defined as the ratio of the total number of broken bonds to the total number of initial bonds at the current moment; the stress concentration factor is defined as the ratio of the maximum von Mises stress to the mean stress in the cell; The macroscopic fracture rate is calculated based on the curve of the cumulative number of broken bonds over time, and the risk dissociation rate per unit of remaining bonds is calculated in combination with the current number of intact bonds. The damage index curve, stress concentration factor curve and risk dissociation rate curve are output to analyze the temporal correlation between damage evolution, stress concentration and bond dissociation risk.

Images