Cell stretching tensor calculation method based on dissipative particle dynamics

Abstract

The application belongs to the technical field of micromanipulation biological modeling, and is based on a cell stretching tensor calculation method of dissipative particle dynamics, and the specific steps are as follows: step 1), a cell membrane model of dissipative particle dynamics is established; step 2), a simulation boundary condition and a type of ensemble integrator are determined; step 3), a cell membrane and a cytoskeleton model are aggregated; step 4), a microneedle penetration experiment platform is built in a simulation environment and the microneedle penetration experiment is started; and step 5), a strain of the cell in a penetration direction is calculated, the microneedle penetration simulation data are processed, finite strain theory is used for analysis and calculation, and a normal strain is output. The method can calculate the internal strain of the cell model under different parameters, and the effectiveness of the method is verified by comparison with the results in a real experiment.

Description

Technical Field

[0001] This invention belongs to the field of micromanipulation biological modeling technology, specifically relating to a method for calculating cell stretching tensors based on dissipative particle dynamics. Background Technology

[0002] Micromanipulation is a technique that involves manipulating cells or early embryos using micromanipulators under a high-powered microscope. As an important manipulation method in cell engineering, micromanipulation is widely used in the manipulation of early embryos and cells.

[0003] With the rapid development of computer and robotics technologies, robotic micromanipulation techniques have become increasingly mature. Compared to manual operation, robotic micromanipulation does not require professional personnel and features low contamination, high repeatability, and a significantly improved success rate.

[0004] However, biological cells undergo mechanical deformation during micromanipulation. If this deformation is not mitigated in time, it can lead to structural damage, affecting subsequent development and, in severe cases, cell death. Therefore, improving cell survival rate and embryonic developmental potential after micromanipulation is a key issue in promoting micromanipulation methods.

[0005] Various studies have shown that intracellular strain is a key link between micromanipulation and cell damage. Some existing methods for sensing intracellular strain require physical markers or sensing devices to be attached to or within the cell. Point-tracking-based visual methods can effectively sense strain through two-dimensional images without damaging the cell, but they cannot provide strain in three-dimensional space. In this context, we can obtain the internal strain of cell models and discover patterns through simulation. Therefore, a model that can reflect cell deformation during manipulation is essential. Summary of the Invention

[0006] This invention addresses the problems existing in the prior art and provides a method for calculating cell stretching tensors based on dissipative particle dynamics. This method can calculate the internal strain of a cell model under different parameters. Comparison with results from real-world experiments verifies the effectiveness of this method and provides guidance for practical experimental operations.

[0007] To achieve the above objectives, the present invention provides the following technical solution: a method for calculating the cell stretching tensor based on dissipative particle dynamics, comprising the following steps:

[0008] Step 1) Establish a dissipative particle dynamics cell membrane model:

[0009] Step 2) Determine the simulation boundary conditions and ensemble integrator type: In the dissipative particle dynamics simulation software, determine the boundaries, boundary types, and ensemble of the simulation environment in three dimensions;

[0010] Step 3) The cell membrane and cytoskeleton model are polymerized;

[0011] Step 4) Set up the microneedle insertion experimental platform in the simulation environment and start the microneedle insertion experiment.

[0012] 4.1 Load the microneedle and suction needle module into the simulation environment. The microneedle is a hollow cylindrical structure with radius Ri, inner diameter ri, and length Li; the suction needle is a cylinder with radius Rs, inner diameter rs, and length Ls.

[0013] 4.2 Load the cell model into the simulation environment and place it between the microneedles and the suction needle, with a 0.5 simulation distance between them. At the same time, make the needle openings of the microneedles and the suction needle face each other and place them along the cell diameter.

[0014] 4.3 Add membrane damage conditions to the cell model so that the sides of the triangular facets on the cell membrane break after exceeding a threshold.

[0015] 4.4 Add a force field as a negative pressure at the needle nozzle to fix the position of the cell model; after the microneedle is stabilized, apply a fixed velocity Vi to the microneedle to move towards the center of the cell, set the time step ts and the number of simulation steps, and perform the simulation.

[0016] Step 5) Calculate the strain of the cell in the direction of insertion.

[0017] The simulation data was processed by inserting microneedles, and the finite strain theory was used for analysis and calculation to output the normal strain.

[0018] Furthermore, the specific steps of step 5) are as follows:

[0019] 5.1 First, sample the three-dimensional simulation results from step 4) and project the sampled results onto a two-dimensional plane;

[0020] 5.2 Using the first frame of the results as the original configuration and the frame before cell damage as the current configuration, calculate the distance difference vector d between each particle i and its neighboring particle j in the cell model in both frames. 0 d;

[0021] 5.3 Calculate the local affine transformation matrix J, where ∑ j∈N |d 0 Jd|=0, where matrix J is the deformation gradient tensor in the discrete state;

[0022] 5.4 The transformation matrix J is decomposed into a unique set of orthogonal local rotation tensors R and stretch tensors U using the polar decomposition theorem, where U represents the ratio of the distance difference in the current configuration to the distance difference in the initial configuration. Strain is defined as E = UI, where I is the identity matrix. Then the strain tensor E is the ratio of the elongation of the distance difference in the current configuration to the distance difference in the initial configuration, and its main diagonal elements are the normal strains in the X, Y, and Z directions.

[0023] Furthermore, in step 1), the cell membrane model structure is represented by a series of triangular meshes on a sphere. The vertices of the triangular facets are represented by particles from dissipative particle dynamics, and the vertices are connected by spring edges. When modeling energy constraints, elastic energy, bending energy, and the area constraint energy and volume constraint energy of the mesh are considered.

[0024] Furthermore, in step 2), periodic boundaries are used in all three directions to ensure the conservation of matter within the system, and a microcanonical system is selected for the ensemble to ensure the conservation of energy in the system.

[0025] Furthermore, the polymerization / depolymerization model in step 3) employs a bonding / dissociation rate expression derived from the Boltzmann distribution based on affinity. When the cross-linked protein particles and microfilament particles are sufficiently close, the bonding / dissociation rate can be k... on A bond is formed at a rate of k; conversely, if the length of an existing bond exceeds the breaking distance, it is broken at a rate of k. off fracture.

[0026] Compared with existing technologies, the beneficial effects of this invention are as follows: This invention can be used to calculate the cell stretching tensor in cell models of dissipative particle dynamics. By building a microneedle insertion experimental platform in a simulation environment and conducting microneedle insertion experiments, the internal strain of the cell model under different parameters is calculated to reflect the cell stretching tensor. Comparison with results from visual methods in real-world experiments demonstrates the effectiveness of the calculation method of this invention and its guiding significance for actual experimental operations. Attached Figure Description

[0027] Figure 1 This is a schematic diagram of the cell model of the present invention.

[0028] Figure 2 This is a cloud map showing the results of microneedle insertion experiments and two-dimensional projected strain calculations using a cell model.

[0029] Figure 3 This is a schematic diagram of strain calculation using the visual optical flow method of this invention.

[0030] Figure 4 This is a schematic diagram of the strain distribution in the X direction during a specific experiment in the example.

[0031] Figure 5 This is a schematic diagram of the strain distribution in the X direction during the simulation of the example. Detailed Implementation

[0032] To enable those skilled in the art to better understand the technical solution of the present invention, the present invention will be described in detail below with reference to the accompanying drawings and specific embodiments.

[0033] A method for calculating the cell stretching tensor based on dissipative particle dynamics includes the following steps:

[0034] Step 1) Establish a dissipative particle dynamics cell membrane model

[0035] The cell membrane is composed of various lipids and proteins, and is a viscous liquid substance. Therefore, it has a certain degree of viscoelasticity; at the same time, it can maintain the shape of the cell and also has a certain degree of resistance to compression and bending.

[0036] The constructed cell model uses a pig oocyte as an example. Its membrane structure is represented by a series of triangular meshes on a sphere. Particles in dissipative particle dynamics are used to represent the vertices of the triangular facets, and the vertices are connected by spring edges. To reflect the mechanical properties of the cell membrane, elastic energy, bending energy, and the area and volume constraints of the mesh are considered in the energy constraint modeling.

[0037] Step 2) Determine the simulation boundary conditions and the type of ensemble integrator.

[0038] In the dissipative particle dynamics simulation software, the boundaries, boundary types, and ensembles of the simulation environment in three dimensions are determined respectively. The boundary type is a periodic boundary or a non-fixed boundary, and the ensemble is a canonical ensemble (NVT), a microcanonical ensemble (NVE), an isothermal-isobaric ensemble (NPT), or an isobaric-isohanic ensemble (NPH).

[0039] To approximate the simulation system, this invention uses periodic boundaries in all three directions to ensure the conservation of matter within the system, and the ensemble uses a microcanonical system to ensure the conservation of energy in the system.

[0040] Step 3) The cell membrane and cytoskeleton model are polymerized;

[0041] The polymerization / depolymerization model uses an affinity-based Boltzmann distribution to derive the bonding / dissociation rate expression. When the cross-linked protein particles and microfilament particles are sufficiently close, the bonding / dissociation rate can be k... on A bond is formed at a rate of k. Conversely, if the length of an existing bond exceeds the breaking distance, a bond is formed at a rate of k. off fracture.

[0042] Step 4) Set up the microneedle insertion experimental platform in the simulation environment and start the microneedle insertion experiment.

[0043] 4.1 Load the microneedle and suction needle module into the simulation environment. The microneedle is a hollow cylindrical structure with radius Ri, inner diameter ri, and length Li; the suction needle is a cylinder with radius Rs, inner diameter rs, and length Ls.

[0044] 4.2 Load the cell model into the simulation environment and place it between the microneedles and the suction needle, with a 0.5 simulation distance between them. At the same time, make the needle openings of the microneedles and the suction needle face each other and place them along the cell diameter.

[0045] 4.3 Add membrane damage conditions to the cell model so that the sides of the triangular facets on the cell membrane break after exceeding a threshold.

[0046] 4.4 Add a force field as a negative pressure at the needle nozzle to fix the position of the cell model; after the microneedle is stabilized, apply a fixed velocity Vi to the microneedle to move towards the center of the cell, set the time step ts and the number of simulation steps, and perform the simulation.

[0047] Step 5) Calculate the strain of the cell in the insertion direction. Process the simulation data of the microneedle insertion, analyze and calculate using finite strain theory, and output the normal strain. The specific steps are as follows:

[0048] 5.1 First, sample the three-dimensional simulation results from step 4) and project the sampled results onto a two-dimensional plane;

[0049] 5.2 Using the first frame of the results as the original configuration and the frame before cell damage as the current configuration, calculate the distance difference vector d between each particle i and its neighboring particle j in the cell model in both frames. 0 d;

[0050] 5.3 Calculate the local affine transformation matrix J, where ∑ j∈N |d 0 When Jd|=0, the matrix J is regarded as the deformation gradient tensor in the discrete state;

[0051] 5.4 The transformation matrix J is decomposed into a unique set of orthogonal local rotation tensors R and stretch tensors U using the polar decomposition theorem. The physical meaning of U can be expressed as the ratio of the distance difference in the current configuration to the distance difference in the initial configuration. Therefore, strain is defined as E = UI, where I is the identity matrix. The physical meaning of the strain tensor E is the ratio of the elongation of the distance difference in the current configuration to the distance difference in the initial configuration, with its main diagonal elements representing the normal strain in the X, Y, and Z directions.

[0052] The above describes the calculation of cell strain in the insertion direction by establishing a cell model based on dissipative particle dynamics, building a microneedle insertion experimental platform in a simulation environment, and initiating a microneedle insertion experiment. To verify the effectiveness of the simulation experiment, it is necessary to compare and verify the results with those of actual microneedle insertion experiments.

[0053] The experimental strain results were plotted as a histogram and normalized, and treated as distribution P0; the simulation results were plotted as a histogram and normalized, and treated as distribution P1.

[0054] The JS divergence is used to measure the similarity between two distributions, while the two-sample KS test and the chi-square test are used to measure the similarity between two distributions.

[0055] JS divergence is a statistical measure based on KL divergence, which measures the similarity between two distributions. Generally, JS divergence is symmetric, with values ​​between 0 and 1, where 0 represents identical distributions and 1 represents identical opposite distributions. It is defined as:

[0056]

[0057]

[0058] Where P1 and P2 are two distribution vectors.

[0059] Wasserstein distance is a method for assessing the similarity length between two distributions. It is defined as the minimum cost required to transform a P distribution into a Q distribution, that is, the minimum work required to transform the shape of one distribution into the shape of another, and can reflect its morphological characteristics to a certain extent.

[0060] The two-sample KS test is generally used to test whether two data distributions are consistent, and is expressed as:

[0061] D m,n =sup|F 1,n (x)-F 2,m (x)|

[0062] Where D m,n F refers to the maximum margin between the cumulative distribution functions of two samples. 1,n (x) refers to the empirical distribution function, sup represents taking the maximum difference, and n and m refer to the size of the two sample data. Its calculation p... value The larger the value of p, the greater the probability that the two belong to the same distribution. It is generally believed that p... value The threshold is 0.05.

[0063] The above provides a general overview of the main process and principle of the cell stretching tensor calculation method based on dissipative particle dynamics of this invention, and also presents the corresponding verification method. The effectiveness of the method is illustrated below through a specific implementation process.

[0064] This embodiment uses the dissipative particle dynamics simulation software Lammps (Large-scale Atomic / Molecular Massively Parallel Simulator) for cell modeling and micromanipulation, Ovito Basic for data computation and post-processing, and Matlab for data processing and visualization. The specific process is as follows:

[0065] Step 1) Write a MATLAB script to model the cell membrane.

[0066] The dsphere(p,0,0,0,1) function is used to determine the geometry of the sphere, calculate the signed distance of the sphere centered at coordinates 0,0,0 with a radius of 1, and returns fd.

[0067] Use [p,t] = distmeshsurface(fd,@huniform,0.2,1.1*[-1,-1,-1;1,1,1]) as the spherical triangulation function, where:

[0068] The first parameter is the fd mentioned above;

[0069] The second parameter is an anonymous function that returns the desired side length;

[0070] The third parameter, h0, is the initial distribution's mesh step size, which determines the mesh fineness.

[0071] The fourth parameter is two arrays, given the bounding boxes of the region: [xmin, ymin, zmin: xmax, ymax, zmax]

[0072] The function's first return value is a set p containing all points after triangulation on the sphere, including the coordinates (x, y, z) of each point in space. The second return value is a set t containing all triangular facets, including the three edges of each facet. Each edge is considered a bond, and each bond is represented by the id of two particles.

[0073] Perform a linear transformation on all points, and all points in p will be mapped to a sphere with radius r of 4.

[0074] Iterate through all the edges of the triangular facets and add all the edges to the key set bdata. Sort the particle IDs of each edge in descending order, use the unique(bdata) function to remove duplicates, and calculate the following attributes:

[0075] The length d of each side;

[0076] The area s enclosed by each triangular facet;

[0077] The volume v of the tetrahedron formed by each triangular facet and its center point (0, 0, 0);

[0078] The dihedral angle formed by two triangular facets;

[0079] The calculation results are stored in the result set.

[0080] The output results are in accordance with the Lammps model file format and include:

[0081] Number of particles, bonds, angles, and dihedral angles

[0082] The number of types of particles, bonds, angles, and dihedral angles

[0083] Output all parameters of the bond on a single line, including id, bond type (bondtype), bond type name, energy unit (temp), ratio of initial bond length to maximum bond length (r0), maximum bond length (rmax), shear modulus (μ0), Upow coefficient (qp), and dissipative force parameters (gamc and gamt). The parameters are as follows: bond uniqueness id, 1, wlc, 0.0036, 0.45, 0.5, 300, 2, 30, 90.

[0084] Output all parameters of the key angle on a single line, including id, angle type (angletype), angle type name, hydrostatic elastic energy coefficient (Cq), hydrostatic elastic energy exponent (q), global area constraint coefficient (Ka), initial global area of ​​the cell (Atot0), volume constraint coefficient (kv), initial volume of the cell (Vtot), local area constraint coefficient (kd), and initial local area of ​​the cell (A0). The parameters are as follows: unique angle id, 1, rbc, 0.0, 1, 7500, 200.92, 7500, 267, 300, 0.24.

[0085] Output the parameters for each dihedral angle, including id, the initial angle θ0 of the adjacent triangle, and the cell bending energy coefficient kb, where id starts from 1, the cell bending energy coefficient is 65, and the angle is the result of the dihedral angle dihe calculated above.

[0086] Output the attributes of each particle, including id, particle tag, particle type, and 3D coordinates x, y, z. The id starts from 1, the particle tag is 1, the particle type is 1, and the 3D coordinates are the results from set p in the above calculations.

[0087] Output the attributes of each key, including the id of particles i and j.

[0088] Output the attributes of each bond corner, including the ID of particles i, j, and k.

[0089] Output the attributes of each dihedral angle, including the ID of particles i, j, k, and l.

[0090] Step 2) Determine the size of the simulation box and boundary conditions.

[0091] Due to computational limitations, in order to approximate the simulation system, the X, Y, and Z axis boundaries are set as periodic boundaries to ensure the conservation of particle number, momentum, and energy in the simulation, and the simulation box size is set from (-8,-8,-8) to (8,8,8).

[0092] In molecular dynamics simulations, all particles within the simulated system constitute an ensemble. Due to the different requirements of the simulation and the different objects of study, the ensembles selected are not the same. In dissipative particle dynamics simulations, the NVE (microcanonical ensemble) is commonly used, which can ensure that the number of particles, system volume, and energy within the system remain constant. The NVE ensemble is set using the command fix 1all nve.

[0093] Step 3) Combining the cytoskeleton model with the cell membrane model

[0094] Define the interaction type between particles in the system, that is, determine the type of potential function that governs the interactions between particles. This simulation system uses dissipative particle dynamics, a type of particle potential used for mesoscopic fluid modeling and biological system modeling. The potential function is defined using the command: `pair_styledpd 1.0 2 1`.

[0095] Because this cell model has multiple components, mixed potential energy is needed to describe the relationships between bonding particles. Two-body potentials for bonds, three-body potentials for angles, and four-body potentials for dihedral angles are set for the particles connecting the cell membrane and cytoskeleton. Constraints are applied using the following commands:

[0096] bond_style hybrid wlc harmonic,

[0097] angle_style hybrid rbc harmonic

[0098] dihedral_style bend

[0099] Set the dpd parameters for particles of different components, including conservative force parameters and dissipative force parameters. Use the command to set the conservative force and dissipative force parameters of the membrane components to: 100, 45, 0.5. The command is:

[0100] pair_coeff 1 1 100 45 0.5;

[0101] Use the command to set the parameters of the skeletal particles to 100, 65, 0.5. The command is:

[0102] pair_coeff 2 2 100 65 0.5;

[0103] The parameters between the two components are set to 100, 45, and 0.5, and the command is:

[0104] pair_coeff 1 2 100 45 1.

[0105] The cell membrane and cells use the affinity-based Boltzmann derivation of the bond dissociation rate expression as the polymerization standard, and employ two commands for polymerization and depolymerization:

[0106] fix createBond all bond / create 1 1 2 0.5 881prob 0.5 235825I param 11jparam 1 2

[0107] fix breakBond all bond / break 1 881 0.5prob 0.5 235825

[0108] k on Set to 10e-4, σ on Set to 10e-4, k s Set to 8000, l0 to 0.5.

[0109] Set the timestep to 0.001 and perform a total of 2,000,000 timesteps.

[0110] The cell model after aggregation is output as the simulation result.

[0111] Step 4) Set up a micromanipulation experimental environment and conduct microneedle insertion experiments.

[0112] 4.1 Load the microneedle and suction needle modules into the simulation environment. The microneedle is a hollow cylinder with a radius of Ri = 0.5, an inner diameter of ri = 0.3, and a length of Li = 4; the suction needle is a hollow cylinder with a radius of Rs = 2.5, an inner diameter of rs = 2.35, and a length of Ls = 2.

[0113] 4.2 Import the cell model into the simulation environment and place it between the microneedle and the holding needle, ensuring the nozzles of the microneedle and the holding needle are facing each other. Add a force field as the holding pressure at the nozzle of the holding needle using the command:

[0114] fix p field1 addforce 0 0 -5

[0115] The position of the fixed cell model is defined, where field1 is the group id of the needle insertion site.

[0116] The holding process consists of 50,000 time steps.

[0117] After the system stabilizes, apply a fixed velocity Vi to the microneedles and move them toward the center of the cells using the following command:

[0118] velocity injector set 0 0 0-${Vi}

[0119] Set the simulation steps, where Vi is the speed, set to 2, 4, 6, 8, 10, and perform 180,000, 90,000, 60,000, 45,000, and 36,000 time steps respectively.

[0120] By command:

[0121] fix membBreak membrane bond / break 10 1 0.5prob 1 12456

[0122] A condition for cell membrane rupture is given: a length threshold is set for all bonds on the cell membrane, meaning that bonds exceeding 0.5 will break.

[0123] The cell model is output as the simulation result.

[0124] Step 5) Calculate intracellular strain

[0125] Since the experiment uses two-dimensional images and optical flow to calculate intracellular strain, this step requires projecting the three-dimensional simulation results onto a two-dimensional plane before performing strain calculations.

[0126] Input the output from the previous step into Ovito Basic.

[0127] In the LAMMPS dump reader, select Edit column mapping, and change the Property of y in File column to Mass or another property that does not affect the calculation result.

[0128] In the pipeline's Data source section, select Simulation cell and set Dimensionality to 2D.

[0129] At this point, the three-dimensional simulation results will be mapped onto the two-dimensional XZ plane.

[0130] Add an Atomic strain modifier, set the cutoff redius to 2, and check the output stretchtensors option.

[0131] The current configuration is defined as the moment just before the cell membrane breaks, and the first frame is defined as the initial configuration.

[0132] Add a compute property modifier, set the Output property to Eeng, and enter StretchTensor.YY-1 in Expression.

[0133] Add a Histogram modifier, select Eeng for Property, and set the Number of histogram bins to 20.

[0134] Check the "fix x-range" box to set the range to between -1 and 0.

[0135] Output the bar chart result, such as Figure 5 As shown.

[0136] The following Python script is written for verification, calculating the distance.

[0137] You need to import the NumPy and SciPy Python packages. For simulation data, denote it as `dataSim`, and for experimental data, denote it as `dataRel`. Then use the methods in SciPy:

[0138] To calculate the JavaScript divergence, use the function: scipy.sptial.distance.jensenshannon(dataSim,dataRel);

[0139] To calculate the Wasserstein distance, use the function: scipy.stats.wasserstein_distance(dataSim,dataRel);

[0140] Use the function to perform a KS test on the data: scipy.stats.ks_2samp(dataSim,dataRel).

[0141] The above metrics were used to measure the similarity between the degree of deformation of simulated cells and that of cells in actual experiments. The calculated JS divergence, Wasserstein distance, and KS test pValue are shown in the table below.

[0142]

[0143]

[0144] Verification has shown that the calculation method of this invention is effective and has certain guiding significance for actual experimental operations.

[0145] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A method for cell stretch tensor calculation based on dissipative particle dynamics, characterized in that, According to the following steps: Step 1) Establish a dissipative particle dynamics cell membrane model; The cell membrane model structure is represented by a series of triangular meshes on a spherical surface, and the vertices of the triangular meshes are represented by particles using dissipative particle dynamics. The vertices are connected by spring edges, and the energy constraints include elastic energy, bending energy, and area and volume constraint energies of the mesh; Step 2) Determine the simulation boundary conditions and the type of ensemble integrator: determine the boundary, boundary type, and ensemble in the three-dimensional direction of the simulation environment in the dissipative particle dynamics simulation software; Step 3) Aggregation and disaggregation of the cell membrane model and the cytoskeleton model; Using the bond formation / dissociation rate expressions derived from the affinity-based Boltzmann distribution, a bond can form at a rate of when the cross-linking protein particle is close enough to the microfilament particle, and conversely, a bond can break at a rate of if the length of the existing bond exceeds the breaking distance. Step 4) Build a microneedle penetration experiment platform in the simulation environment and start the microneedle penetration experiment; 4.1 Load the microneedle and holding needle module into the simulation environment. The microneedle is a hollow cylindrical structure with a radius of Ri, an inner diameter of ri, and a length of Li. The holding needle is a cylindrical structure with a radius of Rs, an inner diameter of rs, and a length of Ls. 4.2 Load the cell model into the simulation environment and place it between the microneedle and the holding needle with a mutual spacing of 0.5 simulation distance. At the same time, make the needle openings of the microneedle and the holding needle opposite and place them along the diameter of the cell. 4.3 Add a membrane damage condition to the cell model, which causes the edges of the triangular facets on the cell membrane to break when their lengths exceed a threshold value. 4.4 Add a force field to the holding needle opening as a holding negative pressure to fix the position of the cell model. After the holding is stable, give the microneedle a fixed size of velocity Vi and move it towards the center of the cell. Set the time step ts and the number of simulation steps to perform the simulation. Step 5) Calculate the strain of the cell in the penetration direction; Process the microneedle penetration simulation data, analyze and calculate using finite strain theory, and output the normal strain.

2. The dissipative particle dynamics based cell stretch tensor calculation method of claim 1, wherein: The specific steps of step 5) are as follows: 5.1 First, sample the three-dimensional simulation results in step 4) and project the sampling results onto a two-dimensional plane; 5.2 Take the first frame in the result as the original configuration, and the frame before the cell membrane is broken as the current configuration, and calculate the distance difference vector of each particle of the cell model and its surrounding neighbor particles in the two frames respectively , , , ; 5.3 Computing the local affine transformation matrix where , the matrix is the deformation gradient tensor in discrete state; 5.4 Decomposition of the transformation matrix using polar decomposition theorem into a unique set of orthogonal local rotation tensors and stretch tensors where is expressed as the ratio of the distance difference in the current configuration to the distance difference in the initial configuration, the strain is defined as where is the identity matrix, the strain tensor is the ratio of the elongation of the distance difference in the current configuration to the distance difference in the initial configuration, and its principal diagonal elements are the normal strains in the X, Y, Z directions.

3. The dissipative particle dynamics based cell stretch tensor calculation method of claim 2, wherein: In step 2), periodic boundaries are used in all three directions to ensure mass conservation in the system, and the microcanonical ensemble is used to ensure energy conservation of the system energy.

Images