Design method of structural mechanical gradient regulation and control interbody fusion cage aiming at end plate defect

Abstract

The invention discloses a design method of a structural mechanical gradient regulation and control interbody fusion cage for end plate defects. The design method comprises the following steps of data acquisition and model construction; analyzing the shape of the defective end plate; carrying out interface mechanical adaptation design on the defect end plate and the structural gradient fusion cage; performing multi-objective optimization; the invention relates to a 3D printing defect endplate mechanical and biological adaptive spinal fusion cage. The individual end plate defect model can be established according to the difference of different patients, the local elastic modulus of the fusion cage is dynamically adjusted through morphological analysis and typing, the stress conduction height matching with the defect end plate is realized, the settlement risk is reduced, and the fusion effect is improved; by combining the biological adaptation design of the end plate, the pore structure and material characteristics of the fusion cage are optimized, and bone healing is promoted; through the scheme of image analysis, digital design, intelligent optimization and precise manufacturing, a mechanical and biological dual-adaptive fusion cage system solution is provided for a patient with the defect of the end plate.

Description

Technical Field

[0001] The present invention relates to the field of medical engineering, and in particular to a design method for a structural mechanics gradient-controlled intervertebral fusion device for end plate defects. Background Art

[0002] The changes in the magnetic resonance signals of the bony endplate and adjacent bone marrow under the vertebral cartilage endplate are called Modic changes, which are more common in middle-aged and elderly patients with spinal degeneration. Abnormal lumbar endplate structure and biomechanical changes have been shown to be associated with lumbar instability and degeneration. Aging leads to the loss of proteoglycans in the intervertebral disc, resulting in structural degeneration, which cannot provide buffering and intervertebral mobility. The endplate is subjected to increased mechanical loads, resulting in microfractures. If the pressure on the endplate exceeds its threshold, it will cause fractures of the endplate trabeculae. Imaging shows local collapse of the endplate or discontinuity of the endplate. A small number of clinical studies have shown that the proportion of endplate collapse and fusion device sinking after posterior lumbar fusion in the endplate defect group is higher than that in the group with normal endplate structure. During the fusion process, the large difference in elastic modulus will cause slight sinking of the fusion device. If there is a defect in the endplate where the fusion device is implanted, it may increase the sinking rate of the fusion device, and it will completely sink into the vertebral body before bony fusion, which will secondarily lead to pseudoarthrosis formation, osteophyte formation and displacement of the fusion device, which is not conducive to patient prognosis and rehabilitation, and increases the risk of secondary surgery.

[0003] Xiao Yang and others from the Department of Orthopedics at West China Hospital of Sichuan University believe that fusion failure or fusion device sinking can be prevented or reduced by preoperative evaluation of endplate sclerosis, reduction of iatrogenic endplate injury, fine treatment of intervertebral disc space, management of osteoporosis and selection of appropriate fusion devices. At present, the most widely used implant in clinical practice is polyetheretherketone fusion device, but it is a bioinert material and has a long fixed shape. After implantation in the body, it is easy to cause biomechanical mismatch at the defect of differentiated endplates, resulting in sinking of the fusion device, and fibrous scar hyperplasia of the fusion device-bone interface, causing loosening of the fusion device. In addition, bone fusion in the intervertebral disc also takes longer, causing adjacent segment degeneration and increasing post-fusion surgery. Through an 8-year follow-up cohort study, South Korea and other countries found that the reoperation rate of 60-69 years old reached 13.2%, and the reoperation rate of spinal fusion after degenerative disease in the United States was 9% to 20%, which greatly increased the burden and economic pressure on patients and easily caused conflicts between doctors and patients.

[0004] Different endplate defects in different patients lead to individual and complicated stress conduction between the endplate and the fusion device. The polyetheretherketone fusion device with uniform structure commonly used in clinical practice makes the stress conduction and bone fusion between the differentiated defect endplate and the fusion device poor, which is prone to vertebral collapse and fusion device sinking, affecting the prognosis and rehabilitation of patients. Therefore, providing a spinal fusion device that can match the complex interface stress conduction of patients with endplate defects, reduce the risk of fusion device sinking and promote intervertebral fusion is a technical problem that needs to be solved urgently by technicians in this field. Summary of the invention

[0005] In view of the above problems existing in the prior art, the present application provides a design method for a structural mechanics gradient-regulated intervertebral fusion device for endplate defects, which can design a spinal fusion device that matches the complex interface stress conduction of patients with endplate defects, reduce the risk of fusion device subsidence, and promote intervertebral fusion.

[0006] The technical solution of the present invention is as follows:

[0007] A design method for a structural mechanics gradient-regulated intervertebral fusion device for endplate defects includes the following steps:

[0008] S1. Data acquisition and model construction, the specific steps are as follows:

[0009] S1-1. Collect CT images and MRI images of the patient's lumbar spine;

[0010] S1-2. Import the CT images and MRI images into Mimics software. If the file format is DICOM, use the CT images for bone structure analysis, perform gray threshold segmentation on the bone tissue in the CT images to obtain the cartilage endplate region, extract the data of the cartilage endplate region and save it as an STL file, that is, the three-dimensional model of the endplate; otherwise, use the MRI images for soft tissue lesion analysis, obtain the T1 and T2 signals of the MRI images, and preliminarily determine the Modic lesion regions: Modic I, Modic II, Modic III; Modic I is the edema region, Modic II is the fatty infiltration region, and Modic III is the sclerosis region;

[0011] S2. Use Matlab software to perform morphological analysis on the cartilage endplate region, the specific steps are as follows:

[0012] S2-1. Perform lesion morphological analysis on the cartilage endplate region, the specific steps are as follows:

[0013] S2-1-1. Use Canny edge detection to extract the initial boundary, that is: read the MRI image of the cartilage endplate region, apply the Canny algorithm to detect the edge, and generate a binary endplate edge map;

[0014] S2-1-2. Fill the boundary holes by closing operation, that is: perform morphological closing operation on the edge map to fill the discontinuous regions and obtain the complete endplate contour;

[0015] S2-1-3. Extract the coordinates of the contour points from the complete endplate contour, save them as an ordered point sequence, and obtain the endplate contour point set Edges[];

[0016] S2-1-4. Load the three-dimensional point cloud data and extract the mesh information, the specific steps are as follows:

[0017] S2-1-4-1. Use the stlread function to load the three-dimensional model of the endplate, obtain the vertices and faces of the three-dimensional model of the endplate, and save the two to the vertex array vertices[] and the face array faces[] respectively;

[0018] S2-1-4-2. Determine the number of segments Edge-Sum of the upper surface edge of the endplate contour point set Edges[] in S2-1-3 for loop control of subsequent thickness calculation;

[0019] S2-1-4-3. Call the size function to count the total number of faces Mesh-Sum, that is, Mesh-Sum = size(faces, 1), for traversal of curvature analysis;

[0020] S2-1-5. Thickness distribution analysis, the specific steps are as follows:

[0021] S2-1-5-1. Initialize the thickness calculation loop value i = 0, the total number of loops is Edge-Sum, and define the thickness storage array thickness_values[];

[0022] S2-1-5-2. Obtain the coordinates (x i , y i , z i _up) of the i-th point on the upper surface edge line in the contour point set, and obtain the coordinates (x i , y i , z i _down) of the corresponding position point on the lower surface edge line of this point;

[0023] S2-1-5-3. Calculate the local thickness value, that is: use the Euclidean distance formula to calculate the vertical distance between the upper and lower surface edge lines at the i-th point according to (x i , y i , z i _up) and (x i , y i , z i _down), and save the result to thickness_values[i];

[0024] S2-1-5-4. If i ≤ Edge-Sum, then i = i + 1 and execute S2-1-5-2, otherwise execute 2-1-5-5;

[0025] S2-1-5-5. Generate a thickness distribution map, that is: map thickness_values[] to the surface of the three-dimensional model of the endplate, use the scatter3 function to generate thickness_map, and visualize it;

[0026] S2-1-6. Curvature analysis, the specific steps are as follows:

[0027] S2-1-6-1. Initialize the mesh traversal loop value j = 0, and define the storage arrays for Gaussian curvature K and mean curvature H;

[0028] S2-1-6-2. Extract the three vertex coordinates P1, P2, and P3 of the patch j from faces[];

[0029] S2-1-6-3. Calculate the normal vector N of the patch j j , that is: extract the vertex coordinates of the patches adjacent to the patch j from faces[], calculate the normal vector of the patch j through cross product, and normalize it;

[0030] S2-1-6-4. Calculate the vertex vectors P j and Q j ;

[0031] S2-1-6-5. Fit the local quadratic surface, that is: take the patch j as the center, select the adjacent patches to form a set, extract all vertices from the patch set, and fit to obtain the quadratic surface equation z;

[0032] S2-1-6-6. Calculate the principal curvatures (K1, K2) of the local quadratic surface according to the quadratic surface equation z, and calculate the Gaussian curvature K = K1 * K2 and the mean curvature H = (K1 + K2) / 2 according to the principal curvatures (K1, K2);

[0033] S2-1-6-7. If j ≤ Mesh-Sum, execute S2-1-6-2, otherwise execute S2-1-6-8;

[0034] S2-1-6-8. Curvature visualization, that is: first map the Gaussian curvature K and mean curvature H in S2-1-6-6 to the surface of the endplate three-dimensional model, and then use the scatter3 function to generate surfature_map to obtain the visualized Gaussian curvature distribution;

[0035] S2-1-7. Volume analysis, that is: traverse the patches in faces[] and calculate the volume V according to the following formula:

[0036]

[0037] where N j 、P j and Q j are from S2-1-6-3 and S2-1-6-4;

[0038] S2-1-8. Superimposed display of thickness and curvature distribution, i.e., superimpose thickness_map and curvature_map onto the three-dimensional endplate model, and use different color channels to distinguish thickness and curvature information;

[0039] S2-2. Use Matlab to analyze the morphological differences of cartilage endplate defects and classify them. The specific steps are as follows:

[0040] S2-2-1. Based on the curvature thickness_map and thickness thickness_map in S2-1-8, obtain the defect area of interest to form a mask image: seg_mask = curvature_map > 0.5 & thickness_map < mean(thickness_map) * 0.7;

[0041] And extract the total number of non-zero pixels K-Sum in the mask image: K_Sum = sum(seg_mask(:));

[0042] The mask image is a binary image used to extract the image of the endplate defect area of interest. The pixel value 1 is the area to be retained, that is, the non-zero pixel area; the pixel value 0 is the area to be excluded;

[0043] S2-2-2. Set the count n of the non-zero pixel area in the mask image to 1, and define the noise area threshold as 5;

[0044] S2-2-3. Extract the noise value Noise_n of the nth non-zero pixel area. If Noise_n ≥ 50, then remove this noise area from SegMask, that is, SegMask = SegMask \setminus Noise_n;

[0045] S2-2-4. If n ≤ K-Sum, then n = n + 1, and continue to execute S2-2-2, otherwise execute S2-2-5;

[0046] S2-2-5. Define the feature map size A_Sum = K-Sum / 4, set the initial count a of feature extraction in the mask image to 1, and define the pooling window size as 2 * 2;

[0047] S2-2-6. Apply the maxpool function to sample the mask image to extract high-level features Feature_a = maxpool(FeatureMap, [2, 2], 'Stride', 2), and calculate the following morphological parameters:

[0048] (1) Edge Depression Ed(a) = maximum width of the defect area / maximum width of the entire endplate; the Edge Depression Ed(a) is used to measure the degree of depression at the edge of the defect area and is defined as the ratio of the maximum width of the defect area to the maximum width of the entire endplate;

[0049] (2) Shape Factor SF(a) = 4π * area of the defect area / square of the perimeter of the defect area; the Shape Factor is used to describe the degree to which the defect area approaches a circular shape, and its value range is (0, 1], being 1 for a perfect circle;

[0050] (3) Irregularity Ratio Ra(a) = maximum length of the defect area / maximum width of the defect area; the Irregularity Ratio is used to describe the degree of elongation of the defect area and is defined as the ratio of the length of the major axis of the area to the length of the minor axis;

[0051] (4) Circularity C(a) = average edge curvature / maximum edge curvature; the Circularity is used to evaluate the smoothness of the edge through the curvature distribution and is defined as the ratio of the average edge curvature to the maximum edge curvature:

[0052] S2-2-7. Integrate the above morphological parameters into a multi-dimensional feature vector: MorphoFeatures_a = [Ed(a), SF(a), Ra(a), C(a)], where a on the left side of the equal sign represents the serial number; obtain the feature matrix M based on A_Sum multi-dimensional feature vectors o = [MorphoFeatures_1, …, MorphoFeatures_A_Sum];

[0053] S2-2-8. If a ≤ A - Sum, then a = a + 1 and continue to execute S2-2-6; otherwise, use the normalize function in Matlab to normalize the feature matrix M o to obtain the normalized feature matrix M and execute S2-2-9;

[0054] S2-2-9. Define the defect morphology Labels = ['depressed type', 'triangular type', 'circular type','rectangular type', 'irregular'], construct and train the neural network Softmax classifier. The functions and codes used in Matlab are as follows:

[0055]

[0056] Extract the weight matrix W and the bias term b from net:

[0057] W = net.layers(4).Weights;

[0058] b = net.Layers(4).Bias;

[0059] S2-2-10. Morphological classification prediction, with input parameters: normalized feature matrix M, weight parameter W, and bias term b, and the output result: morphological classification Class_M = ['concave type', 'triangular type', 'round type','rectangular type', 'irregular'];

[0060] S2-2-11. Defect degree prediction, with input parameters: thickness_map, curvature_map, and V; calculating the characteristic parameters: average thickness T and average curvature C; judging the defect degree Type based on the average thickness T and average curvature C, that is: Type1, Type2, Type3;

[0061] S2-2-12. Combining the output results of S2-2-10 and S2-2-11 to obtain the classification label: [morphological classification & defect degree];

[0062] S2-3. Embedding the classification label into Abaqus, and the specific steps are as follows:

[0063] S2-3-1. Obtaining the classification label from S2-2 and matching the predefined material property parameter library Material according to the defect degree Type:

[0064]

[0065] where E represents the elastic modulus and v represents the Poisson's ratio;

[0066] S2-3-2. Three-dimensional mesh data conversion, and the specific steps are as follows:

[0067] S2-3-2-1. Converting the vertex array vertices[] into the node matrix Nodes[] of Abaqus;

[0068] S2-3-2-2. Converting the triangular patches in faces[] into tetrahedrons through Matlab, calling the TetGen tool to generate a high-quality mesh, and saving the mesh vertex coordinates;

[0069] S2-3-2-3. Assigning numerical values to the material properties obtained in S2-3-1;

[0070] S2-3-2-4. Generating an Abaqus input file in inp format, and the Abaqus input file includes: classification label, node matrix Nodes[], mesh vertex coordinates, and material properties;

[0071] S3: Mechanical adaptation design of the defective endplate and the fusion cage. The fusion cage refers to the structural mechanics gradient regulated intervertebral fusion cage. The specific steps are as follows:

[0072] S3-1. Design and modeling of the fusion cage. The specific steps are as follows:

[0073] S3-1-1. Design of the elliptical recessed unit cell structure. The specific steps are as follows:

[0074] S3-1-1-1. Set the number of elliptical nodes i-Sum = 4, and define the major axis length a and minor axis length b of the ellipse. The value ranges are [1.8, 3] and [0.8, 2] respectively; a starts from 1.8 and is incremented by 0.2, with a total of 7; b starts from 0.8 and is incremented by 0.2, with a total of 7.

[0075] S3-1-1-2. The auxetic parameters are: α = 0.4a, β = 0.2b, k = 6.

[0076] S3-1-1-3. Calculate the node coordinates (X i Y i ) of the elliptical recessed unit cell structure according to the following formula:

[0077]

[0078] S3-1-1-4. After cycling through 4 nodes, retain the node coordinates of the elliptical recessed unit cell structure and calculate the aspect ratio AR = a / b of the elliptical recessed unit cell structure.

[0079] S3-1-1-5. Superimpose the chiral offset, that is:

[0080]

[0081] S3-1-2. Stacking of the hybrid structure layers, that is: set L as the layer value, and L = [1, 2, 3, 4]; set the total number of layers L-Sum = 4; stack the unit structures layer by layer, adjust the shape control function f(L) = 1 + 0.1*(L - 1), and output the node coordinates of the hybrid structure:

[0082]

[0083] where w L is the weight coefficient of each layer, used to control the layer contribution ratio;

[0084] S3-1-3. After performing biological constraint screening on the hybrid structure in S3-1-2, output the node coordinates and elastic modulus of the hybrid structure.

[0085] S3-2. Interface matching design and calculation of the defective endplate and the fusion cage. The specific steps are as follows:

[0086] S3-2-1. Perform a preliminary modeling of the elastic modulus distribution of the fusion device, that is: traverse the entire space to calculate the elastic modulus of the fusion device at the position (x, y). where E min is the minimum elastic modulus, representing the flexibility limit of the material; E max is the maximum elastic modulus, representing the rigidity limit of the material; d is the gradient adjustment coefficient, controlling the steepness of the change in elastic modulus; T(x, y) is the endplate thickness and morphological eigenvalue at the corresponding position, from the thickness_map in S2-1-8; T0 is the set key thickness adjustment point, used to define the maximum elastic modulus change region.

[0087] S3-2-2. Interface matching optimization, the specific steps are as follows:

[0088] S3-2-2-1. Set the upper limit value of the loop count Er-Sum.

[0089] S3-2-2-2. Define the elastic modulus of the defective endplate at the position (x, y) as E endplate (x, y); Compare E r (x, y) with E endplate (x, y), and calculate the interface matching function value ΔE(x, y):

[0090]

[0091] where r is the loop count, and the value range is [1, Er-Sum].[[]]

[0092] S3-2-2-3. If ΔE(x, y) <= 10%, then save the elastic modulus E r (x, y) at the position r, otherwise modify the elastic modulus E r (x, y) of the fusion device at this place as E r (x, y) + γ·ΔE(x, y)·E endplate (x, y), where γ is the adjustment parameter.

[0093] S3-2-2-4. Repeat S3-2-2-2 and S3-2-2-3 to obtain the final E r (x, y);

[0094] S3-2-4. Match the E r (x, y) in S3-2-2 with the node coordinates and elastic modulus of the hybrid structure in S3-1-3 to confirm the structural design of the fusion device.

[0095] S3-3. Abaqus biomechanics and fluid mechanics analysis, the specific steps are as follows:

[0096] Perform operations similar to S2-3 to generate the Abaqus input file 2 in inp format that contains the defective endplate and the fusion device, including all node coordinates, mesh vertex coordinates, and material properties of the hybrid structure;

[0097] Import the Abaqus input file 2 into Abaqus and perform biomechanical and hydrodynamics analysis and calculations, with the output data being "specific surface area, settlement risk, porosity, peak endplate stress, peak fusion device stress";

[0098] S4. Complete multi-objective optimization and solve for the optimal structure solution set in Matlab. The specific steps are as follows:

[0099] S4-1. Determine the optimization objectives and design variables; the optimization objectives are the Abaqus output data of S3-3, namely: maximum specific surface area f1(x), minimize settlement risk f2(x), optimize porosity f3(x), peak endplate stress f4(x), peak fusion device stress f5(x);

[0100] The design variables are the fusion device structure design parameters of S3-1, namely:

[0101] Design variable x = [a / b, h, t, θ, r, p, l, α],

[0102] which includes:

[0103] (1) Parameters of the elliptical recessed unit cell structure: elliptical axis ratio a / b, recessed depth h, wall thickness t;

[0104] (2) Parameters of the chiral structure: twist angle θ, twist radius r, chiral unit spacing p;

[0105] (3) Hybrid unit cell design variables: unit cell size l, hybrid ratio α of the elliptical recessed unit cell structure and the chiral structure;

[0106] S4-1-2. Construct the optimization problem and determine the optimization objective. The objective function expression is:

[0107]

[0108] S4-2. Use the weighted method to search for the Pareto solution set and perform Pareto multi-objective front analysis to ensure the balance between mechanical properties and biological properties;

[0109] S4-3. Use the NSGA-II algorithm to initialize the population; calculate the objective function values of each individual; perform non-dominated sorting and crowding degree calculation; perform crossover, mutation, and selection operations to generate the next generation of the population; iterate until convergence to the Pareto optimal front to obtain the optimal solution; the NSGA-II algorithm is the non-dominated sorting genetic algorithm;

[0110] S5. Derive a 3D model file in STL format from the optimal solution of S4-3, and then use 3D printing technology to fabricate the fusion device.

[0111] Furthermore, in step S1-3, the rules for classifying bone tissues according to the gray threshold are as follows:

[0112] If the gray threshold HU < 350, it is defined as the cartilage endplate;

[0113] If the gray threshold 350 < HU < 850, it is defined as cancellous bone;

[0114] If the gray threshold HU > 850, it is defined as cortical bone.

[0115] Furthermore, the specific steps of S1-4 are as follows:

[0116] S1-4-1. Define the thresholds of three groups of T1 and T2 according to the following rules, which are respectively used to determine three Modic lesion regions. The meanings represented by the letters are: I, II, and III represent serial numbers, lower represents the lower limit value, and upper represents the upper limit value:

[0117] I_T1_lower = 50, I_T1_upper = 150;

[0118] I_T2_lower = 700, I_T2_upper = 1200;

[0119] II_T1_lower = 500, II_T1_upper = 800;

[0120] II_T2_lower = 600, II_T2_upper = 900;

[0121] III_T1_lower = 50, III_T1_upper = 200;

[0122] III_T2_lower = 100, III_T2_upper = 300;

[0123] S1-4-2. Initially determine the lesion regions according to the thresholds of T1 and T2. The rules are as follows:

[0124] If T1 = [I_T1_lower, I_T1_upper] and T2 = [I_T2_lower, I_T2_upper], it is determined as Modic I;

[0125] If T1 = [II_T1_lower, II_T1_upper] and T2 = [II_T2_lower, II_T2_upper], then it is determined as Modic II;

[0126] If T1 = [III_T1_lower, III_T1_upper] and T2 = [III_T2_lower, III_T2_upper], then it is determined as Modic III;

[0127] S1-4-2. Dynamically adjust the thresholds of T1 and T2, and the specific steps are as follows:

[0128] S1-4-2-1. Recalculate the thresholds of T1 and T2 in S1-4-1 according to the following method:

[0129] new_threshold_low = μ - k·σ, new_threshold_high = μ + k·σ,

[0130] where new_threshold_low is the lower limit value of the thresholds of T1 and T2, new_threshold_high is the upper limit value of the thresholds of T1 and T2, μ is the mean of the regional gray level, σ is the standard deviation of the regional gray level, and k is the adjustment coefficient;

[0131] S1-4-2-2. Compare new_threshold_low and new_threshold_high with the thresholds of T1 and T2 in S1-4-1. If there is a change, continue to execute S1-4-2-1; otherwise, execute S1-4-2-3;

[0132] S1-4-2-3. Check the neighborhood pixels of the center point of each Modic lesion area. If the gray level values of the neighborhood pixels are within the threshold range of S1-4-2, merge these neighborhood pixels into the Modic lesion area.

[0133] Furthermore, in step S2-2-11, the rules for judging the defect degree according to the average thickness T and the average curvature C are as follows:

[0134] Type1 = [T > 2mm, -0.5 < C < 0.5, V < 50mm 3 ,

[0135] Type2 = [1 < T < 2mm, C < -0.5, 50 < V < 200mm 3 ,

[0136] Type3 = [T < 1mm, C > 0.5, V > 200mm 3 .

[0137] Further, in step S2-3-2-3, the specific values of the material properties are as follows:

[0138] Type1: E1 = 1000 MPa, v1 = 0.3;

[0139] Type2: E2 = 500 MPa, v2 = 0.35;

[0140] Type3: E3 = 200 MPa, v3 = 0.4.

[0141] Further, in step S3-1-3, the rules for biological screening of the hybrid structure are as follows:

[0142] (1) Specific surface area Specific > 12, and the calculation formula is Specific = Asurface / Vsolid, where Asurface is the surface area of the hybrid structure and Vsolid is the solid volume of the hybrid structure;

[0143] (2) Porosity The calculation formula is where Vtotal is the external volume of the hybrid structure, including the intermediate pores.

[0144] The beneficial technical effects of the present invention are as follows:

[0145] (1) Precise biomechanical adaptation, reducing the risk of subsidence and enhancing the fusion effect: Combining the patient's imaging data, an individualized endplate defect model is established. Through morphological analysis and classification, the local elastic modulus of the fusion device is dynamically adjusted to achieve a high degree of matching of stress conduction with the defective endplate, making the stress distribution of the fusion device on the endplate more uniform, reducing local stress concentration, significantly reducing the subsidence risk compared with traditional PEEK fusion devices, and improving the fusion success rate; targeted optimization for different types of endplate defects is carried out to achieve precise treatment; by calculating the gradient distribution of the elastic modulus, the fusion device has different mechanical properties in different parts, better adapting to the endplate lesion area and improving long-term stability.

[0146] (2) Integrating biological and mechanical factors to promote bone healing: Combining the biological adaptation design of the endplate, optimizing the pore structure and material properties of the fusion device, the combined design of elliptical concave single cells and chiral torsion realizes the negative Poisson's ratio effect and controllable porosity, improves the bone ingrowth ability, and promotes intervertebral fusion. The 3D printing technology realizes high-precision manufacturing, making the microstructure of the fusion device more in line with the biomechanical requirements and improving the adaptability and growth environment of bone tissue.

[0147] (3) The present invention provides a fusion device system solution that is biocompatible in both mechanics and biology for patients with endplate defects through the "imaging analysis - digital design - intelligent optimization - precise manufacturing" solution: by matching the T1 / T2 signal thresholds of the Modic lesion area, it reduces iatrogenic damage to the endplate, decreases the migration rate of the fusion device, and shortens the rehabilitation period. Through the combined modeling of Matlab and Abaqus, it realizes the full-process calculation from image data processing to finite element analysis, improving the design accuracy and efficiency. By adopting the multi-objective optimization method, it systematically optimizes key factors such as the mechanical compatibility, anti-settlement property, and bone-bonding performance of the fusion device, ensuring its applicability to different patients and enhancing its clinical practicality. BRIEF DESCRIPTION OF THE DRAWINGS

[0148] Figure 1 is the overall flowchart of the present invention;

[0149] Figure 2 is the logic block diagram of data acquisition and model construction;

[0150] Figure 3 is the logic block diagram of soft tissue lesion analysis;

[0151] Figure 4 is the logic block diagram of the morphological analysis of cartilage endplate lesions;

[0152] Figure 5 is the logic block diagram of the extraction and classification of the differences in the defective endplate;

[0153] Figure 6 is the flowchart of embedding the classification label into Abaqus;

[0154] Figure 7 is the logic block diagram of the design and modeling of the fusion device;

[0155] Figure 8 is the logic block diagram of the interface matching design and calculation between the defective endplate and the fusion device;

[0156] Figure 9 is the flowchart of solving the optimal structure design solution set in the target space in Matlab. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0157] The present invention will be specifically described below with reference to the accompanying drawings and embodiments. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all of them. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0158] As Figure 1 shown, the embodiment is mainly divided into 5 steps:

[0159] S1. Data collection and model construction;

[0160] S2. Analysis of the morphology of the defective endplate;

[0161] S3. Mechanical adaptability design of the interface between the defective endplate and the structural gradient fusion device;

[0162] S4. Multi-objective optimization;

[0163] S5. 3D printing of a spinal fusion device with mechanical and biological adaptability for the defective endplate.

[0164] I. The specific steps of data collection and model construction are as Figures 2 - 3 shown:

[0165] S1-1. Collect CT images and MRI images of the patient's lumbar spine;

[0166] S1-2. Import the CT images and MRI images into Mimics software. If the file format is DICOM, use the CT images for bone structure analysis to obtain the bony region and proceed to S1-3; otherwise, use the MRI images for soft tissue lesion analysis and proceed to S1-4;

[0167] S1-3. Conduct bony region analysis as Figure 2 shown, that is, use Mimics software to perform gray-scale threshold segmentation on the bone tissue in the CT images to obtain the region of the cartilage endplate; the rules for classifying bone tissue according to the gray-scale threshold are as follows:

[0168] If the gray-scale threshold HU < 350, it is defined as the cartilage endplate;

[0169] If the gray-scale threshold 350 < HU < 850, it is defined as cancellous bone;

[0170] If the gray-scale threshold HU > 850, it is defined as cortical bone;

[0171] Extract the data of the cartilage endplate and save it as an STL file, that is, the three-dimensional model of the endplate;

[0172] S1-4. Conduct soft tissue region analysis as Figure 3 shown, obtain the T1 and T2 signals of the MRI images, and preliminarily determine the Modic lesion regions: Modic I, Modic II, Modic III; Modic I is the edema region, Modic II is the fatty infiltration region, and Modic III is the sclerosis region; the specific steps are as follows:

[0173] S1-4-1. Define the thresholds of three groups of T1 and T2 according to the following rules, which are respectively used to determine three Modic lesion regions. The meanings represented by the letters are as follows: I, II, and III represent serial numbers, lower represents the lower limit value, and upper represents the upper limit value:

[0174] I_T1_lower = 50, I_T1_upper = 150;

[0175] I_T2_lower = 700, I_T2_upper = 1200;

[0176] II_T1_lower = 500, II_T1_upper = 800;

[0177] II_T2_lower = 600, II_T2_upper = 900;

[0178] III_T1_lower = 50, III_T1_upper = 200;

[0179] III_T2_lower = 100, III_T2_upper = 300;

[0180] S1-4-2. Initially determine the lesion region according to the thresholds of T1 and T2. The rules are as follows:

[0181] If T1 = [I_T1_lower, I_T1_upper] and T2 = [I_T2_lower, I_T2_upper], then it is determined as Modic I;

[0182] If T1 = [II_T1_lower, II_T1_upper] and T2 = [II_T2_lower, II_T2_upper], then it is determined as Modic II;

[0183] If T1 = [III_T1_lower, III_T1_upper] and T2 = [III_T2_lower, III_T2_upper], then it is determined as Modic III;

[0184] S1-4-2. Dynamically adjust the thresholds of T1 and T2. The specific steps are as follows:

[0185] S1-4-2-1. Recalculate the thresholds of T1 and T2 in S1-4-1 according to the following method:

[0186] new_threshold_low = μ - k·σ, new_threshold_high = μ + k·σ,

[0187] where new_threshold_low is the lower limit of the T1 and T2 thresholds, new_threshold_high is the upper limit of the T1 and T2 thresholds, μ is the mean of the regional gray level, σ is the standard deviation of the regional gray level, and k is an adjustment coefficient;

[0188] S1-4-2-2. Compare new_threshold_low and new_threshold_high with the thresholds of T1 and T2 in S1-4-1. If there is a change, continue to execute S1-4-2-1; otherwise, execute S1-4-2-3.

[0189] S1-4-2-3. Check the neighborhood pixels of the center point of each Modic lesion area. If the gray values of the neighborhood pixels are within the threshold range of S1-4-2, merge these neighborhood pixels into the Modic lesion area.

[0190] Step 1 combines CT and MRI image data, and through a reasonable segmentation strategy and dynamic optimization algorithm, realizes the accurate identification of the bony structure and Modic lesion area, completes the accurate modeling of the endplate and its lesion area, and provides reliable data support for subsequent biomechanical research and clinical diagnosis. (1) The reliability of CT gray threshold method for bone tissue segmentation: HU value has a clear physical meaning in bone tissue segmentation, can accurately distinguish cartilage endplate, cancellous bone and cortical bone, and provides high-quality geometric data for subsequent finite element analysis. (2) The feasibility of MRI-based Modic lesion classification: T1 and T2 weighted imaging have been widely used for Modic lesion typing, and can better reflect the histological characteristics of edema, fat infiltration and sclerosis area. The traditional fixed threshold method may not be able to adapt to the individual differences of different patients, affecting the accuracy of lesion identification; therefore, MRI-based Modic typing combined with dynamic threshold adjustment improves the robustness of lesion area identification and avoids misjudgment caused by individual differences; further combining the gray information of neighborhood pixels for region merging helps to reduce artifact interference and makes lesion identification more stable.

[0191] II. The specific steps of the defective endplate morphology analysis are as Figures 4 - 6 shown:

[0192] S2-1. As Figure 4 shown, perform lesion morphology analysis on the cartilage endplate area. The specific steps are as follows:

[0193] S2-1-1. Use Canny edge detection to extract the initial boundary, that is: read the MRI image of the cartilage endplate area, apply the Canny algorithm to detect the edge, and generate a binary endplate edge map.

[0194] S2-1-2. Perform closing operation to fill boundary holes, that is: perform morphological closing operation on the edge map to fill discontinuous regions and obtain the complete endplate contour;

[0195] S2-1-3. Extract the coordinates of the contour points from the complete endplate contour, save them as an ordered point sequence, and obtain the endplate contour point set Edges[];

[0196] S2-1-4. Load the 3D point cloud data and extract mesh information. The specific steps are as follows:

[0197] S2-1-4-1. Use the stlread function to load the 3D endplate model, obtain the vertices and faces of the 3D endplate model, and save them to the vertex array vertices[] and the face array faces[] respectively;

[0198] S2-1-4-2. Determine the number of segments Edge-Sum of the upper surface edge of the endplate contour point set Edges[] in S2-1-3 for loop control in subsequent thickness calculation;

[0199] S2-1-4-3. Call the size function to count the total number of faces Mesh-Sum, that is, Mesh-Sum = size(faces,1), for traversal in curvature analysis;

[0200] S2-1-5. Thickness distribution analysis. The specific steps are as follows:

[0201] S2-1-5-1. Initialize the loop value i = 0 for thickness calculation, the total number of loops is Edge-Sum, and define the thickness storage array thickness_values[];

[0202] S2-1-5-2. Obtain the coordinates (x i , y i , z i _up) of the i-th point on the upper surface edge line in the contour point set, and obtain the coordinates (x i , y i , z i _down) of the corresponding point on the lower surface edge line;

[0203] S2-1-5-3. Calculate the local thickness value, that is: use the Euclidean distance formula to calculate the vertical distance between the upper and lower surface edge lines at point i according to (x i , y i , z i _up) and (x i , y i , z i _down), and save the result to thickness_values[i];

[0204] S2-1-5-4: If i ≤ Edge-Sum, then i = i + 1 and execute S2-1-5-2; otherwise, execute 2-1-5-5;

[0205] S2-1-5-5: Generate a thickness distribution map, that is, map thickness_values[] to the surface of the endplate 3D model, use the scatter3 function to generate thickness_map, and visualize it;

[0206] S2-1-6: Curvature analysis, the specific steps are as follows:

[0207] S2-1-6-1: Initialize the mesh traversal loop value j = 0, and define storage arrays for Gaussian curvature K and mean curvature H;

[0208] S2-1-6-2: Extract the three vertex coordinates P1, P2, and P3 of patch j from faces[];

[0209] S2-1-6-3: Calculate the normal vector N of patch j j , that is, extract the vertex coordinates of the patches adjacent to patch j from faces[], calculate the normal vector of patch j through cross product, and normalize it;

[0210] S2-1-6-4: Calculate the vertex vectors P j and Q j ;

[0211] S2-1-6-5: Fit a local quadratic surface, that is, take patch j as the center, select adjacent patches to form a set, extract all vertices from the patch set, and fit to obtain the quadratic surface equation z;

[0212] S2-1-6-6: Calculate the principal curvatures (K1, K2) of the local quadratic surface according to the quadratic surface equation z, and calculate the Gaussian curvature K = K1 * K2 and the mean curvature H = (K1 + K2) / 2 according to the principal curvatures (K1, K2);

[0213] S2-1-6-7: If j ≤ Mesh-Sum, then execute S2-1-6-2; otherwise, execute S2-1-6-8;

[0214] S2-1-6-8: Curvature visualization, that is, first map the Gaussian curvature K and mean curvature H in S2-1-6-6 to the surface of the endplate 3D model, and then use the scatter3 function to generate surfature_map to obtain the visualized Gaussian curvature distribution;

[0215] S2-1-7: Volume analysis, that is, traverse the patches in faces[] and calculate the volume V according to the following formula:

[0216]

[0217] where N j 、P j and Q j are from S2-1-6-3 and S2-1-6-4;

[0218] The superposition display of the thickness and curvature distributions of S2-1-8, that is: superimpose the thickness_map and the thickness_map onto the three-dimensional model of the endplate, and use different color channels to distinguish the thickness and curvature information;

[0219] S2-2, as Figure 5 shown, use Matlab to analyze the morphological differences of the cartilage endplate defects and perform classification. The specific steps are as follows:

[0220] S2-2-1, Based on the curvature thickness_map and the thickness thickness_map of S2-1-8, obtain the defect region of interest and form a mask image: seg_mask = curvature_map>0.5 & thickness_map<mean(thickness_map)*0.7;

[0221] And extract the total number of non-zero pixels K-Sum in the mask image: K_Sum = sum(seg_mask(:));

[0222] The said mask image is a binary image, used to extract the image of the endplate defect region of interest. The pixel value 1 is the region to be retained, that is, the non-zero pixel region; the pixel value 0 is the region to be excluded;

[0223] S2-2-2, Set the count n of the non-zero pixel region of the mask image to 1 and define the noise area threshold 5;

[0224] S2-2-3, Extract the noise value Noise_n of the nth non-zero pixel region. If Noise_n≥50, then remove this noise region from SegMask, that is, SegMask = SegMask\setminus Noise_n;

[0225] S2-2-4, If n≤K-Sum, then n = n + 1 and continue to execute S2-2-2, otherwise execute S2-2-5;

[0226] S2-2-5, Define the feature map size A_Sum = K-Sum / 4, set the initial count a of the feature extraction of the mask image to 1, and define the pooling window size 2*2;

[0227] S2-2-6. Sample the masked image using the maximum pooling function maxpool to extract high-level features Feature_a = maxpool(FeatureMap, [2, 2], 'Stride', 2), and calculate the following morphological parameters:

[0228] (1) Edge depression degree Ed(a) = maximum width of the defect area / maximum width of the entire endplate; the edge depression degree Ed(a), i.e., Edge Depression, is used to measure the depression degree of the defect area at the edge, defined as the ratio of the maximum width of the defect area to the maximum width of the entire endplate;

[0229] (2) Shape factor SF(a) = 4π * area of the defect area / square of the perimeter of the defect area; the shape factor, i.e., Shape Factor, is used to describe the degree to which the defect area approaches a circle, and its value range is (0, 1], and it is 1 for a perfect circle;

[0230] (3) Irregularity ratio Ra(a) = maximum length of the defect area / maximum width of the defect area; the irregularity ratio, i.e., Irregularity Ratio, is used to describe the elongation degree of the defect area, defined as the ratio of the length of the major axis of the area to the length of the minor axis;

[0231] (4) Circularity C(a) = average edge curvature / maximum edge curvature; the circularity, i.e., Circularity, evaluates the edge smoothness through the curvature distribution, defined as the ratio of the average edge curvature to the maximum edge curvature:

[0232] S2-2-7. Integrate the above morphological parameters into a multi-dimensional feature vector: MorphoFeatures_a = [Ed(a), SF(a), Ra(a), C(a)], where a on the left side of the equal sign represents the serial number; obtain the feature matrix M according to A_Sum multi-dimensional feature vectors o = [MorphoFeatures_1, …, MorphoFeatures_A_Sum];

[0233] S2-2-8. If a ≤ A - Sum, then a = a + 1 and continue to execute S2-2-6; otherwise, normalize the feature matrix M using the normalize function in Matlab o to obtain the normalized feature matrix M, and execute S2-2-9;

[0234] S2-2-9. Define the defect morphologies Labels = ['depressed type', 'triangle type', 'round type','rectangle type', 'irregular'], construct and train a neural network Softmax classifier. The functions and codes used in Matlab are as follows:

[0235]

[0236] Extract the weight matrix W and the bias term b from the net:

[0237] W = net.layers(4).Weights;

[0238] b = net.Layers(4).Bias;

[0239] S2-2-10. Morphological classification prediction, with input parameters: the standardized feature matrix M, the weight parameter W, and the bias term b, and the output result: the morphological classification Class_M = ['concave type', 'triangular type', 'circular type','rectangular type', 'irregular'];

[0240] S2-2-11. Defect degree prediction, with input parameters: thickness_map, curvature_map, and V; calculate the characteristic parameters: the average thickness T and the average curvature C; judge the defect degree Type based on the average thickness T and the average curvature C:

[0241] Type1 = [T > 2mm, -0.5 < C < 0.5, V < 50mm 3 ,

[0242] Type2 = [1 < T < 2mm, C < -0.5, 50 < V < 200mm 3 ,

[0243] Type3 = [T < 1mm, C > 0.5, V > 200mm 3 ;

[0244] S2-2-12. Combine the output results of S2-2-10 and S2-2-11 to obtain the classification label: [morphological classification & defect degree];

[0245] S2-3. As Figure 6 shown, embed the classification label into Abaqus, and the specific steps are as follows:

[0246] S2-3-1. Obtain the classification label from S2-2, and match the predefined material property parameter library Material according to the defect degree Type:

[0247]

[0248] where E represents the elastic modulus and v represents the Poisson's ratio;

[0249] S2-3-2. Three-dimensional mesh data conversion, and the specific steps are as follows:

[0250] S2-3-2-1. Convert the vertex array vertices[] into the node matrix Nodes[] of Abaqus;

[0251] S2-3-2-2. Use Matlab to convert the triangular patches in faces[] into tetrahedrons, call the TetGen tool to generate high-quality meshes, and save the mesh vertex coordinates;

[0252] S2-3-2-3. Assign numerical values to the material properties obtained in S2-3-1:

[0253] Type1: (High stiffness) E1 = 1000 MPa, v1 = 0.3;

[0254] Type2: (Medium stiffness) E2 = 500 MPa, v2 = 0.35;

[0255] Type3: (Low stiffness) E3 = 200 MPa, v3 = 0.4.

[0256] S2-3-2-4. Generate an Abaqus input file in inp format, and the Abaqus input file includes: classification labels, the node matrix Nodes[], mesh vertex coordinates, and material properties.

[0257] Step 2 Systematically combine image analysis, morphological modeling, and finite element simulation, which not only improves the recognition accuracy of endplate lesions but also enhances the practicality of biomechanical analysis, providing strong technical support for the precise treatment of spinal degenerative diseases. (1) Endplate defect morphology analysis and classification: Canny edge detection has good anti-noise performance and can accurately obtain endplate edge information. Calculate the endplate surface morphology using Gaussian curvature and mean curvature to comprehensively evaluate the depression degree and smoothness of the endplate damage area, thereby supporting subsequent lesion typing. Use Matlab to extract features of the defect morphology, including edge depression degree, shape factor, irregularity, and circularity, etc., and construct a multi-dimensional feature matrix to comprehensively describe the endplate defect morphology and avoid classification bias that may be caused by a single index. (2) Endplate lesion type and biomechanical modeling: Through the data conversion interface between Matlab and Abaqus, convert CT / MRI data into the mesh data required for finite element analysis to form an imaging analysis and biomechanical simulation process.

[0258] III. The specific steps of the interface mechanical adaptation design between the defective endplate and the structural gradient fusion device are as Figure 7 、 8 shown:

[0259] S3-1. Fuser design and modeling ( Figure 7 ), and the specific steps are as follows:

[0260] S3-1-1. Design of elliptical concave single-cell structure, the specific steps are as follows:

[0261] S3-1-1-1. Set the number of elliptical nodes i-Sum = 4, define the major axis length a and minor axis length b of the ellipse, and the value ranges are [1.8, 3] and [0.8, 2] respectively; a starts from 1.8 and is incremented by 0.2, with a total of 7; b starts from 0.8 and is incremented by 0.2, with a total of 7;

[0262] S3-1-1-2. The auxetic parameters are: α = 0.4a, β = 0.2b, k = 6;

[0263] S3-1-1-3. Calculate the node coordinates (X i Y i ) of the elliptical concave single-cell structure according to the following formula:

[0264]

[0265] S3-1-1-4. After cycling through 4 nodes, retain the node coordinates of the elliptical concave single-cell structure, and calculate the aspect ratio AR = a / b of the elliptical concave single-cell structure;

[0266] S3-1-1-5. Superimpose the chiral offset, that is:

[0267]

[0268] S3-1-2. Stacking of the hybrid structure layers, that is: set L as the layer value, and L = [1, 2, 3, 4]; set the total number of layers L-Sum = 4; stack the unit structures layer by layer, adjust the shape control function f(L) = 1 + 0.1*(L - 1), and output the node coordinates of the hybrid structure:

[0269]

[0270] where w L is the weight coefficient of each layer, used to control the layer contribution ratio;

[0271] S3-1-3. Screen the hybrid structure of S3-1-2 according to the following constraint conditions for biological constraints:

[0272] (1) Specific surface area Specific > 12, and the calculation formula is Specific = Asurface / Vsolid, where Asurface is the surface area of the hybrid structure and Vsolid is the solid volume of the hybrid structure;

[0273] (2) Porosity The calculation formula is where Vtotal is the external volume of the hybrid structure, including the intermediate pores;

[0274] Then output the node coordinates and elastic modulus of the hybrid structure;

[0275] S3-2. Interface matching design and calculation between the defective endplate and the fusion cage ( Figure 7 ), the specific steps are as follows:

[0276] S3-2-1. Perform a preliminary modeling of the elastic modulus distribution of the fusion cage, that is: traverse the entire space to calculate the elastic modulus of the fusion cage at the position (x, y) where E min is the minimum elastic modulus, representing the flexibility limit of the material; E max is the maximum elastic modulus, representing the rigidity limit of the material; d is the gradient adjustment coefficient, controlling the steepness of the change in elastic modulus; T(x, y) is the endplate thickness and morphological characteristic value at the corresponding position, from the thickness_map of S2-1-8; T0 is the set key thickness adjustment point, used to define the maximum elastic modulus change region;

[0277] S3-2-2. Interface matching optimization ( Figure 8 ), the specific steps are as follows:

[0278] S3-2-2-1. Set the upper limit value of the loop count Er-Sum;

[0279] S3-2-2-2. Define the elastic modulus of the defective endplate at the position (x, y) as E endplate (x, y); compare E r (x, y) with E endplate (x, y), and calculate the interface matching function value ΔE(x, y):

[0280]

[0281] where r is the loop count, and the value range is [1, Er-Sum];

[0282] S3-2-2-3. If ΔE(x, y) <= 10%, then save the elastic modulus Er(x, y) at the r position, otherwise modify the elastic modulus E of this place r (x, y) = E r (x, y) + γ·ΔE(x, y)·E endplate (x, y), where γ is the adjustment parameter;

[0283] S3-2-2-4. Repeat S3-2-2-2 and S3-2-2-3 to obtain the final E r (x, y);

[0284] S3-2-4. The E of S3-2-2 r(x, y) is matched with the node coordinates and elastic modulus of the hybrid structure in S3-1-3 to confirm the structural design of the fusion device;

[0285] S3-3. Abaqus biomechanical and fluid mechanics analysis, the specific steps are as follows:

[0286] Take operations similar to those in S2-3 to generate the Abaqus input file 2 in inp format containing the defective endplate and the fusion device, including all node coordinates, mesh vertex coordinates, and material properties of the hybrid structure;

[0287] Import the Abaqus input file 2 into Abaqus and perform biomechanical and fluid mechanics analysis calculations, and the output data is "specific surface area, settlement risk, porosity, peak stress of the endplate, peak stress of the fusion device".

[0288] In Step 3, through the mechanical adaptation optimization of the structural gradient fusion device and the defective endplate, the fusion device can better adapt to the endplate defect area, reduce stress concentration, and lower the settlement risk; by reasonably adjusting the specific surface area and porosity of the fusion device, sufficient bone tissue growth space is provided to promote bone fusion; combined with Abaqus biomechanical and fluid mechanics analysis, the bearing capacity and adaptability of the fusion device are evaluated to provide reliable data support for clinical applications.

[0289] IV. The specific steps of multi-objective optimization are as Figure 9 shown:

[0290] S4-1. Determine the optimization objectives and design variables; the optimization objectives are the Abaqus output data in S3-3, namely: maximum specific surface area f1(x), minimize settlement risk f2(x), optimize porosity f3(x), peak stress of the endplate f4(x), peak stress of the fusion device f5(x);

[0291] The design variables are the fusion device structural design parameters in S3-1, namely:

[0292] Design variable x = [a / b, h, t, θ, r, p, l, α],

[0293] wherein it includes:

[0294] (1) Parameters of the elliptical recessed unit cell structure: elliptical axis length ratio a / b, recessed depth h, wall thickness t;

[0295] (2) Parameters of the chiral structure: twist angle θ, twist radius r, chiral unit spacing p;

[0296] (3) Hybrid unit cell design variables: unit cell size l, hybrid ratio α of the elliptical recessed unit cell structure and the chiral structure;

[0297] S4-1-2. Construct an optimization problem and determine the optimization objective. The objective function expression is as follows:

[0298]

[0299] S4-2. Use the weight method to search for the Pareto solution set and conduct Pareto multi-objective frontier analysis to ensure the balance between mechanical properties and biological properties.

[0300] S4-3. Use the NSGA-II algorithm to initialize the population; calculate the objective function values of each individual; conduct non-dominated sorting and crowding degree calculation; perform crossover, mutation, and selection operations to generate the next generation of the population; iterate until convergence to the Pareto optimal frontier to obtain the optimal solution. The NSGA-II algorithm is a non-dominated sorting genetic algorithm.

[0301] In Step 4, through a multi-objective optimization strategy, five core performance indicators calculated by Abaqus are selected as the optimization objectives to construct a multi-objective optimization problem, ensuring that the fusion device not only has good mechanical support capabilities but also takes into account performance requirements such as bone ingrowth and biocompatibility. The design variables cover multiple key parameters such as the unit cell geometry, chiral structure, and hybrid structure of the fusion device, so as to adjust the structural parameters during the optimization process to meet the biomechanical requirements. Finally, the optimal design of the structural parameters of the fusion device is achieved, enabling it to have high bone integration ability, good mechanical adaptability, and low subsidence risk, laying an important foundation for the future design of personalized spinal fusion devices.

[0302] V. Derive a 3D model file in STL format according to the optimal solution of S4-3, and then use 3D printing technology to fabricate the fusion device.

[0303] Although the embodiments of the present invention are disclosed as above, they are not limited to only the applications listed in the specification and embodiments. It can be fully applied to various fields suitable for the present invention. For those skilled in the art, for ordinary technical personnel in the art, various changes, modifications, substitutions, and variations can be made to these embodiments without departing from the principles and spirit of the present invention. Therefore, without departing from the general concept defined by the claims and their equivalents, the present invention is not limited to specific details.

Claims

1. A design method of a structural mechanics gradient-regulated intervertebral fusion device for endplate defects, characterized in that It includes the following steps: S1. Data collection and model construction, and the specific steps are as follows: S1-1. Collect CT images and MRI images of the patient's lumbar spine; S1-2. Import the CT images and MRI images into Mimics software. If the file format is DICOM, use the CT images for bone structure analysis, perform gray-scale threshold segmentation on the bone tissue in the CT images to obtain the cartilage endplate region, extract the data of the cartilage endplate region and save it as an STL file, that is, the three-dimensional model of the endplate; otherwise, use the MRI images for soft tissue lesion analysis, obtain the T1 and T2 signals of the MRI images, and initially determine the Modic lesion regions: Modic I, Modic II, Modic III; Modic I is the edema region, Modic II is the fatty infiltration region, and Modic III is the sclerosis region; S2. Use Matlab software to perform morphological analysis on the cartilage endplate region, and the specific steps are as follows: S2-1. Perform lesion morphological analysis on the cartilage endplate region, and the specific steps are as follows: S2-1-1. Use Canny edge detection to extract the initial boundary, that is: read the MRI image of the cartilage endplate region, apply the Canny algorithm to detect the edge, and generate a binary endplate edge map; S2-1-2. Perform closing operation to fill the boundary holes, that is: perform morphological closing operation on the edge map to fill the discontinuous regions and obtain the complete endplate contour; S2-1-3. Extract the coordinates of the contour points from the complete endplate contour, save them as an ordered point sequence, and obtain the endplate contour point set Edges[]; S2-1-4. Load the three-dimensional point cloud data and extract the mesh information, and the specific steps are as follows: S2-1-4-1. Use the stlread function to load the three-dimensional model of the endplate, obtain the vertices and faces of the three-dimensional model of the endplate, and save them to the vertex array vertices[] and the face array faces[] respectively; S2-1-4-2. Determine the segmentation number Edge-Sum of the upper surface edge of the endplate contour point set Edges[] in S2-1-3 for loop control of subsequent thickness calculation; S2-1-4-3. Call the size function to count the total number of faces Mesh-Sum, that is, Mesh-Sum = size(faces,1), for traversal of curvature analysis; S2-1-5. Thickness distribution analysis, and the specific steps are as follows: S2-1-5-1. Initialize the thickness calculation loop value i = 0, the total number of loops is Edge-Sum, and define the thickness storage array thickness_values[]; S2-1-5-2. Obtain the coordinates (x i , y i , z i _up) of the i-th point on the upper surface edge line in the contour point set, and obtain the coordinates (x i , y i , z i _down) of the position point corresponding to this point on the lower surface edge line; S2-1-5-3. Calculate the local thickness value, i.e., use the Euclidean distance formula to calculate the vertical distance between the upper and lower surface edge lines at point i based on (x i , y i , z i _up) and (x i , y i , z i _down), and save the result to thickness_values[i]; S2-1-5-4. If i ≤ Edge-Sum, then i = i + 1 and execute S2-1-5-2, otherwise execute 2-1-5-5; S2-1-5-5. Generate a thickness distribution map, that is: map thickness_values[] to the surface of the three-dimensional model of the endplate, use the scatter3 function to generate thickness_map, and visualize it; S2-1-6. Curvature analysis, and the specific steps are as follows: S2-1-6-1. Initialize the grid traversal loop value j = 0, and define the storage arrays for Gaussian curvature K and mean curvature H; S2-1-6-2. Extract the three vertex coordinates P1, P2, and P3 of patch j from faces[]; S2-1-6-3. Calculate the normal vector N of patch j j , that is: extract the vertex coordinates of the patches adjacent to patch j from faces[], calculate the normal vector of patch j through cross product, and normalize it; S2-1-6-4. Calculate the vertex vectors P j and Q j ; S2-1-6-5. Fit the local quadratic surface, that is: taking patch j as the center, select adjacent patches to form a set, extract all vertices from the patch set, and fit to obtain the quadratic surface equation z; S2-1-6-6. Calculate the principal curvatures (K1, K2) of the local quadratic surface according to the quadratic surface equation z, and calculate the Gaussian curvature K = K1 * K2 and the mean curvature H = (K1 + K2) / 2 according to the principal curvatures (K1, K2); S2-1-6-7. If j ≤ Mesh-Sum, execute S2-1-6-2, otherwise execute S2-1-6-8; S2-1-6-8. Curvature visualization, that is: first map the Gaussian curvature K and mean curvature H in S2-1-6-6 to the surface of the endplate three-dimensional model, and then use the scatter3 function to generate surfature_map to obtain the visualized Gaussian curvature distribution; S2-1-7. Volume analysis, that is: traverse the patches of faces[] and calculate the volume V according to the following formula: where N j , P j and Q j are from S2-1-6-3 and S2-1-6-4; S2-1-8. Superimposed display of thickness and curvature distribution, that is: superimpose thickness_map and thickness_map on the endplate three-dimensional model, and use different color channels to distinguish thickness and curvature information; S2-2. Use Matlab to analyze the morphological differences of cartilage endplate defects and perform classification. The specific steps are as follows: S2-2-1. Based on the curvature thickness_map and thickness thickness_map in S2-1-8, obtain the defect area of interest and form a mask image: seg_mask = curvature_map>0.5 & thickness_map<mean(thickness_map)*0.7; Extract the total number of non-zero pixels K-Sum in the mask image: K_Sum = sum(seg_mask(:)); The mask image is a binary image used to extract the image of the endplate defect area of interest. The pixel value 1 is the area to be retained, that is, the non-zero pixel area, and the pixel value 0 is the area to be excluded; S2-2-2. Set the non-zero pixel area count n = 1 of the mask image and define the noise area threshold 5; S2-2-3. Extract the noise value Noise_n of the nth non-zero pixel area. If Noise_n ≥ 50, remove this noise area from SegMask, that is, SegMask = SegMask\setminus Noise_n; S2-2-4. If n ≤ K-Sum, then n = n + 1 and continue to execute S2-2-2, otherwise execute S2-2-5; S2-2-5. Define the feature map size A_Sum = K-Sum / 4, set the initial count a = 1 for feature extraction of the mask image, and define the pooling window size 2*2; S2-2-6. Apply the maxpool function, maxpool, to sample the masked image, extract the high-level feature Feature_a = maxpool(FeatureMap, [2, 2], 'Stride', 2), and calculate the following morphological parameters: (1) Edge Depression Ed(a) = maximum width of the defect area / maximum width of the entire endplate; the Edge Depression Ed(a), i.e., Edge Depression, is used to measure the degree of depression at the edge of the defect area and is defined as the ratio of the maximum width of the defect area to the maximum width of the entire endplate; (2) Shape Factor SF(a) = 4π * area of the defect area / square of the perimeter of the defect area; the Shape Factor is used to describe the degree to which the defect area approaches a circular shape, and its value range is (0, 1], being 1 for a perfect circle; (3) Irregularity Ratio Ra(a) = maximum length of the defect area / maximum width of the defect area; the Irregularity Ratio is used to describe the degree of elongation of the defect area and is defined as the ratio of the length of the major axis to the length of the minor axis of the area; (4) Circularity C(a) = average edge curvature / maximum edge curvature; the Circularity is used to evaluate the smoothness of the edge through the curvature distribution and is defined as the ratio of the average edge curvature to the maximum edge curvature: S2-2-7. Integrate the above morphological parameters into a multi-dimensional feature vector: MorphoFeatures_a = [Ed(a), SF(a), Ra(a), C(a)], where a on the left side of the equal sign represents the serial number; obtain the feature matrix M according to A_Sum multi-dimensional feature vectors o = [MorphoFeatures_1, …, MorphoFeatures_A_Sum]; S2-2-8. If a ≤ A-Sum, then a = a + 1 and continue to execute S2-2-6; otherwise, use the normalize function in Matlab to normalize the feature matrix M o to obtain the normalized feature matrix M, and execute S2-2-9; S2-2-9. Define the defect morphologies Labels = ['concave type', 'triangular type', 'circular type','rectangular type', 'irregular'], construct and train the neural network Softmax classifier. The functions and code used in Matlab are as follows: layers = featureInputLayer(size(M, 2), 'Name', 'input' fullyConnectedLayer(64, 'Name', 'fc1') reluLayer('Name','relu') fullyConnectedLayer(5, 'Name', 'output') softmaxLayer('Name','softmax') classificationLayer('Name', 'class')]; net = trainNetwork(M, labels, layers); Extract the weight matrix W and the bias term b from net: W = net.layers(4).Weights; b = net.Layers(4).Bias; S2-2-10. Morphological classification prediction. The input parameters are: the standardized feature matrix M, the weight parameter W, and the bias term b. The output result is: the morphological classification Class_M = ['concave type', 'triangular type', 'circular type','rectangular type', 'irregular']; S2-2-11. Defect degree prediction. The input parameters are: thickness_map, curvature_map, and V. The calculated characteristic parameters are: average thickness T and average curvature C. The defect degree Type is determined according to the average thickness T and average curvature C, i.e., Type1, Type2, Type3. S2-2-12. Combine the output results of S2-2-10 and S2-2-11 to obtain the classification label: [morphological classification & defect degree]. S2-3. Embed the classification label into Abaqus. The specific steps are as follows: S2-3-1. Obtain the classification label from S2-2 and match the predefined material property parameter library Material according to the defect degree Type: where E represents the elastic modulus and v represents the Poisson's ratio. S2-3-2. Three-dimensional mesh data conversion. The specific steps are as follows: S2-3-2-1. Convert the vertex array vertices[] into the node matrix Nodes[] of Abaqus. S2-3-2-2. Convert the triangular patches in faces[] into tetrahedrons through Matlab, call the TetGen tool to generate high-quality meshes, and save the mesh vertex coordinates. S2-3-2-3. Assign numerical values to the material properties obtained in S2-3-1. S2-3-2-4. Generate an Abaqus input file in inp format. The Abaqus input file contains: classification label, node matrix Nodes[], mesh vertex coordinates, and material properties. S3: Mechanical adaptation design of the interface between the defective endplate and the fusion device. The fusion device refers to the structural mechanics gradient-regulated intervertebral fusion device. The specific steps are as follows: S3-1. Fusion device design and modeling. The specific steps are as follows: S3-1-1. Design of the elliptical concave single-cell structure. The specific steps are as follows: S3-1-1-1. Set the number of elliptical nodes i-Sum = 4, define the major axis length a and minor axis length b of the ellipse, and the value ranges are [1.8, 3] and [0.8, 2] respectively; a starts from 1.8 and is incremented by 0.2, with a total of 7; b starts from 0.8 and is incremented by 0.2, with a total of 7. S3-1-1-2. The auxetic parameters are: α = 0.4a, β = 0.2b, k = 6. S3-1-1-3. Calculate the node coordinates (X i Y i ) of the elliptical recessed unit cell structure according to the following formula: S3-1-1-4. After cycling through 4 nodes, retain the node coordinates of the elliptical concave single-cell structure and calculate the aspect ratio AR = a / b of the elliptical concave single-cell structure. S3-1-1-5. Superimpose the chiral offset, i.e.: S3-1-2. Stacking of the hybrid structure layers, i.e.: Set L as the layer value, and L = [1, 2, 3, 4]; set the total number of layers L-Sum = 4; stack the unit structures layer by layer, adjust the shape control function f(L) = 1 + 0.1*(L - 1), and output the node coordinates of the hybrid structure: where w L is the weight coefficient of each layer, used to control the contribution ratio of each layer; S3-1-3. After performing biological constraint screening on the hybrid structure in S3-1-2, output the node coordinates and elastic modulus of the hybrid structure. S3-2. Interface matching design and calculation between the defective endplate and the fusion device. The specific steps are as follows: S3-2-1. Conduct a preliminary modeling of the elastic modulus distribution of the fusion device, that is, traverse the entire space to calculate the elastic modulus of the fusion device at the position (x, y). where E min is the minimum elastic modulus, representing the flexibility limit of the material; E max is the maximum elastic modulus, representing the rigidity limit of the material. d is the gradient adjustment coefficient, which controls the steepness of the change in elastic modulus; T(x, y) is the endplate thickness and morphological eigenvalue at the corresponding position, obtained from the thickness_map in S2-1-8; T0 is the set key thickness adjustment point, used to define the maximum elastic modulus change region; S3-2-2. Interface matching optimization, the specific steps are as follows: S3-2-2-1. Set the upper limit value of the loop count Er-Sum; S3-2-2-2. Define the elastic modulus of the defective endplate at the position (x, y) as E endplate (x, y); compare E r (x, y) with E endplate (x, y), and calculate the interface matching function value ΔE(x, y): where r is the loop count, and the value range is [1, Er-Sum]; S3-2-2-3. If ΔE(x,y) <= 10%, then save the elastic modulus E at the r position r (x,y), otherwise modify the elastic modulus E of the fuser at this location r (x,y) = E r (x,y) + γ·ΔE(x,y)·E endplate (x,y), where γ is an adjustment parameter; S3-2-2-4. Repeat the execution of S3-2-2-2 and S3-2-2-3 to obtain the final E r (x, y); S3-2-4. Match the E in S3-2-2 r (x, y) with the node coordinates and elastic modulus of the hybrid structure in S3-1-3 to confirm the structural design of the fusion device; S3-3. Abaqus biomechanics and fluid mechanics analysis, the specific steps are as follows: Take operations similar to those in S2-3 to generate the Abaqus input file 2 in inp format containing the defective endplate and the fusion device, including all node coordinates, mesh vertex coordinates, and material properties of the hybrid structure; Import the Abaqus input file 2 into Abaqus and perform biomechanics and fluid mechanics analysis calculations, and the output data is "specific surface area, settlement risk, porosity, endplate peak stress, fusion device peak stress"; S4. Complete multi-objective optimization and solve the optimal structure solution set in Matlab, the specific steps are as follows: S4-1. Determine the optimization objectives and design variables; the optimization objectives are the Abaqus output data in S3-3, namely: maximum specific surface area f1(x), minimum settlement risk f2(x), optimized porosity f3(x), endplate peak stress f4(x), fusion device peak stress f5(x); The design variables are the fusion device structure design parameters in S3-1, namely: Design variable x = [a / b, h, t, θ, r, p, l, α], which includes: (1) Parameters of the elliptical recessed unit cell structure: elliptical axis ratio a / b, recessed depth h, wall thickness t; (2) Parameters of the chiral structure: twist angle θ, twist radius r, chiral unit spacing p; (3) Hybrid unit cell design variables: unit cell size l, hybrid ratio α of the elliptical recessed unit cell structure and the chiral structure; S4-1-2. Construct the optimization problem and determine the optimization objective, and the objective function expression is: S4-2. Use the weighted method to search for the Pareto solution set and perform Pareto multi-objective frontier analysis to ensure the balance of mechanical and biological properties; S4-3. Use the NSGA-II algorithm to initialize the population; calculate the objective function values of each individual; perform non-dominated sorting and crowding degree calculation; perform crossover, mutation, and selection operations to generate the next generation population; iterate until convergence to the Pareto optimal frontier to obtain the optimal solution; the NSGA-II algorithm is the non-dominated sorting genetic algorithm; S5. Derive the STL format three-dimensional model file according to the optimal solution in S4-3, and then use 3D printing technology to fabricate the fusion device.

2. The design method of the structural mechanics gradient-regulated intervertebral fusion device for endplate defect according to claim 1, characterized in that: In step S1-3, the rule for classifying bone tissue according to the gray threshold is as follows: If the gray threshold HU < 350, it is defined as the cartilage endplate; If the gray threshold 350 < HU < 850, it is defined as cancellous bone; If the gray-scale threshold HU > 850, it is defined as cortical bone.

3. The design method of the structural mechanics gradient-regulated intervertebral fusion device for endplate defects according to claim 1, characterized in that: The specific steps of S1-4 are as follows: S1-4-1: Define the thresholds of 3 groups of T1 and T2 according to the following rules, which are respectively used to determine 3 Modic lesion regions. The meanings represented by the letters are: I, II, and III represent serial numbers, lower represents the lower limit value, and upper represents the upper limit value: I_T1_lower = 50, I_T1_upper = 150; I_T2_lower = 700, I_T2_upper = 1200; II_T1_lower = 500, II_T1_upper = 800; II_T2_lower = 600, II_T2_upper = 900; III_T1_lower = 50, III_T1_upper = 200; III_T2_lower = 100, III_T2_upper = 300; S1-4-2: Initially determine the lesion region according to the thresholds of T1 and T2. The rules are as follows: If T1 = [I_T1_lower, I_T1_upper] and T2 = [I_T2_lower, I_T2_upper], it is determined as Modic I; If T1 = [II_T1_lower, II_T1_upper] and T2 = [II_T2_lower, II_T2_upper], it is determined as Modic II; If T1 = [III_T1_lower, III_T1_upper] and T2 = [III_T2_lower, III_T2_upper], it is determined as Modic III; S1-4-2: Dynamically adjust the thresholds of T1 and T2. The specific steps are as follows: S1-4-2-1: Recalculate the thresholds of T1 and T2 in S1-4-1 according to the following method: new_threshold_low = μ - k·σ, new_threshold_high = μ + k·σ, where new_threshold_low is the lower limit value of the thresholds of T1 and T2, new_threshold_high is the upper limit value of the thresholds of T1 and T2, μ is the mean of the regional gray scale, σ is the standard deviation of the regional gray scale, and k is the adjustment coefficient; S1-4-2-2: Compare new_threshold_low and new_threshold_high with the thresholds of T1 and T2 in S1-4-1. If there is a change, continue to execute S1-4-2-1; otherwise, execute S1-4-2-3; S1-4-2-3: Check the neighborhood pixels of the center point of each Modic lesion region. If the gray scale values of the neighborhood pixels are within the threshold range of S1-4-2, merge these neighborhood pixels into the Modic lesion region.

4. The design method of the structural mechanics gradient-regulated intervertebral fusion device for endplate defect according to claim 1, characterized in that: In step S2-2-11, the rules for judging the defect degree according to the average thickness T and the average curvature C are as follows: Type1 = [T > 2mm, -0.5 < C < 0.5, V < 50mm 3 , Type2 = [1 < T < 2mm, C < -0.5, 50 < V < 200mm 3 , Type3 = [T < 1mm, C > 0.5, V > 200mm 3 .

5. The design method of the structural mechanics gradient-regulated intervertebral fusion device for endplate defect according to claim 1, characterized in that: In step S2-3-2-3, the specific numerical values of the material properties are as follows: Type1: E1 = 1000MPa, v1 = 0.3; Type2: E2 = 500MPa, v2 = 0.35; Type3: E3 = 200MPa, v3 = 0.

4.

6. The design method of the structural mechanics gradient-regulated intervertebral fusion device for endplate defect according to claim 1, characterized in that: In step S3-1-3, the rules for biological screening of the hybrid structure are as follows: (1) Specific surface area Specific > 12, and the calculation formula is Specific = Asurface / Vsolid, where Asurface is the surface area of the hybrid structure and Vsolid is the solid volume of the hybrid structure; (2) Porosity The calculation formula is where Vtotal is the external volume of the hybrid structure, including the intermediate pores.

Images