From MRI to Model: In Silico Medicine with Wolfram Language

October 7, 2025

Alessandro Mastrofini, Bioengineering Application Developer Intern, Algorithms R&D

In silico medicine, particularly through the use of finite element simulations, is revolutionizing patient-specific healthcare, especially in the musculoskeletal system. By creating detailed computational models of individual anatomy, finite element analysis enables precise simulation of biomechanical behavior under various conditions. This approach allows for personalized treatment strategies, offering predictive insights into how specific interventions—such as surgeries, prosthetics or rehabilitation plans—might affect the patient’s musculoskeletal health. Unlike traditional experiments such as in vivo or in vitro, in silico medicine refers to simulations where the experimental environment is recreated within the processor.

The term in silico comes from the main component of computers: the processor, which is made of silicon. With these simulation techniques, such as finite element simulations, it is possible to reduce the need for invasive procedures and physical testing by providing a noninvasive, cost-effective and highly detailed method to explore the complex interactions within structures like tendons, joints and bones. This innovation is paving the way for more effective, individualized treatments and a deeper understanding of biomechanical health.

In this example, we will explore a detailed demonstration of a finite element uniaxial test experiment applied to a tendon structure. The geometry, derived from medical imaging, highlights patient-specific outcomes. For simplicity, we will start with the anatomical database available in Wolfram Language:

m[x_]:=AnatomyData[x,"Graphics3D"][[1]]
muscles=Join[AnatomyData[(Entity["right calcaneal tendon", Entity["AnatomicalStructure","RightCalcanealTendon"], "Entity[\"AnatomicalStructure\", \"RightCalcanealTendon\"]", "anatomical structure"])]["AttachedMuscles"],AnatomyData[(Entity["right calcaneal tendon", Entity["AnatomicalStructure","RightCalcanealTendon"], "Entity[\"AnatomicalStructure\", \"RightCalcanealTendon\"]", "anatomical structure"])]["InsertedMuscles"]];
(*Graphics3D[{
m[(Entity["right calcaneal tendon", Entity["AnatomicalStructure","RightCalcanealTendon"], "Entity[\"AnatomicalStructure\", \"RightCalcanealTendon\"]", "anatomical structure"])],
Opacity[0.1],m[#]&/@muscles,
m[(Entity["skeleton of right free lower limb", Entity["AnatomicalStructure","SkeletonOfRightFreeLowerLimb"], "Entity[\"AnatomicalStructure\", \"SkeletonOfRightFreeLowerLimb\"]", "anatomical structure"])],
m[#]&/@{(Entity["set of lumbricals of right foot", Entity["AnatomicalStructure","SetOfLumbricalsOfRightFoot"], "Entity[\"AnatomicalStructure\", \"SetOfLumbricalsOfRightFoot\"]", "anatomical structure"]),(Entity["right flexor accessorius", Entity["AnatomicalStructure","RightFlexorAccessorius"], "Entity[\"AnatomicalStructure\", \"RightFlexorAccessorius\"]", "anatomical structure"]),(Entity["right flexor hallucis brevis", Entity["AnatomicalStructure","RightFlexorHallucisBrevis"], "Entity[\"AnatomicalStructure\", \"RightFlexorHallucisBrevis\"]", "anatomical structure"]),(Entity["right adductor hallucis", Entity["AnatomicalStructure","RightAdductorHallucis"], "Entity[\"AnatomicalStructure\", \"RightAdductorHallucis\"]", "anatomical structure"]),(Entity["flexor digiti minimi brevis of right foot", Entity["AnatomicalStructure","FlexorDigitiMinimiBrevisOfRightFoot"], "Entity[\"AnatomicalStructure\", \"FlexorDigitiMinimiBrevisOfRightFoot\"]", "anatomical structure"]),(Entity["plantar interosseous of right foot", Entity["AnatomicalStructure","PlantarInterosseousOfRightFoot"], "Entity[\"AnatomicalStructure\", \"PlantarInterosseousOfRightFoot\"]", "anatomical structure"])},
Opacity[0.05],m[#]&/@{(Entity["skin of right foot", Entity["AnatomicalStructure","SkinOfRightFoot"], "Entity[\"AnatomicalStructure\", \"SkinOfRightFoot\"]", "anatomical structure"]),(Entity["skin of right leg", Entity["AnatomicalStructure","SkinOfRightLeg"], "Entity[\"AnatomicalStructure\", \"SkinOfRightLeg\"]", "anatomical structure"])}},Axes->False,PlotRange->{All,All,{All,500}},Boxed->False,ViewPoint->{-100,100,30}]*)

The process begins with geometry analysis, accommodating both STL-like meshes and volumetric images commonly produced by clinical diagnostic imaging techniques. Following this, we will assign material properties to the tendon structure and simulate a load-bearing experimental scenario, focusing on its behavior during a uniaxial tensile test. This approach demonstrates the workflow for integrating patient-specific data into biomechanical simulations.

Calcaneal Tendon and Virtual Patient

The calcaneal tendon, or Achilles tendon, is the strongest and largest tendon in the human body, connecting the gastrocnemius and soleus muscles of the calf to the calcaneus (heel bone). It plays a vital role in enabling plantarflexion of the foot, a movement essential for walking, running and jumping. This tendon is capable of withstanding tensile loads of up to 12.5 times body weight during high-impact activities, such as running or jumping, and has an elastic modulus ranging from 1.2 to 2.0 GPa, allowing it to store and release energy efficiently. Typically, the tendon can tolerate strains of up to 8–10% before failure.

The Achilles tendon is prone to pathologies caused by overuse, acute trauma and age-related degeneration. Chronic overuse often leads to Achilles tendinopathy, characterized by pain, swelling and impaired function due to repetitive microtrauma. Sudden, forceful movements in sports can result in tendon ruptures, particularly in middle-aged individuals with diminished elasticity. Insertional tendinitis, inflammation near the tendon’s attachment to the heel, is another common issue, exacerbated by biomechanical stressors like Haglund’s deformity or improper footwear. Aging further reduces tendon elasticity, increasing susceptibility to microtears and ruptures, while overloading and biomechanical imbalances, such as overpronation, unevenly distribute stress, heightening the risk of injury.

In this context, in silico medicine plays a crucial role in clinical scenarios. By using computational models tailored to individual patients, clinicians can simulate the biomechanical behavior of the Achilles tendon under various conditions. This approach enables the prediction of injury risk, optimization of treatment strategies and design of personalized rehabilitation plans, all while reducing reliance on invasive procedures. The integration of in silico methods into clinical practice enhances our ability to address the complex challenges of Achilles tendon pathologies and improve patient outcomes.

The following example is structured into multiple sections. First, the “Imaging” section outlines the transformation of medical images into a virtual representation of the patient, utilizing several anatomical marker points. Next, the image is converted into a mesh, a structure suitable for finite element analysis (FEA). Finally, multiple FEAs are performed to examine the mechanical behavior trends and the impact of the pathology.

Imaging

The standard approach to patient-specific biomechanical assessment begins with medical imaging:

tendonMesh=AnatomyData[AnatomyData[(Entity[calcaneal tendon, Entity["AnatomicalStructure","CalcanealTendon"], "Entity[\"AnatomicalStructure\", \"CalcanealTendon\"]", "anatomical structure"]),"BilateralComponents"][[1]],"MeshRegion"]

While handling the various medical image formats is beyond the scope of this text, all formats can be converted into grayscale volumetric images, where a specific threshold can be set to highlight the relevant tissue regions. To provide a broadly applicable example, we will start with the tendon morphology obtained from the AnatomyData of the anatomical structure entity in Wolfram Language:

tendonImage=ImagePad[RegionImage[Rotate[tendonMesh,-5Degree,{1,0,0}]],10,0];

Image3D[tendonImage]

tendonTrimImage=ImageTrim[tendonImage,{{0,0,60},{180,180,180}}]

This image showcases the volume rendering of the tendon, closely resembling the images produced by medical scans such as CT or MRI.

Next, we can convert this volumetric representation into a mesh suitable for FEA. The goal is to generate a mesh that is both sufficiently accurate and computationally efficient, balancing detail with a reasonable number of elements.

To achieve this, we can use the ImageMesh function with different methods. However, in all cases, the resulting mesh retains characteristics tied to the resolution of the original volumetric image. This is a common challenge when converting clinical images from medical equipment into digital models suitable for numerical computations:

(tendonRemeshed=ImageMesh[Closing[MeanFilter[tendonTrimImage,2],10],Method->#,CornerNeighbors->True])&/@{Automatic, "MarchingSquares"}

A more effective approach involves analyzing the mesh structure in terms of its points. These points, now serving as markers on the tendon surface, can then be refined using smoothing algorithms, alternative meshing techniques or, given the tendon’s shape, a surface lofting approach to create a more structured representation.

Anatomical Markers

At this stage, it becomes clear that we can identify multiple contour lines at the same height, which can be utilized to construct a lofted surface:

pts=MeshCoordinates[tendonRemeshed];

Show[tendonRemeshed,Graphics3D[Point[pts]]]

The approach follows these key steps:

Define a fixed number of slices: Establish evenly spaced cross sections along the tendon’s height.

Extract points for each slice: For each cross section, collect all points at the corresponding height.

Generate contour lines: Create a line that interpolates through the extracted points for each slice.

Construct the loft: Connect the sequence of contour lines to form a smooth 3D lofted surface.

This method enables a structured and refined surface representation, improving the accuracy and smoothness of the final mesh while maintaining anatomical fidelity.

For convenience, let’s start with a rotation of the tendon to align it with a more intuitive coordinate system:

actualWidth=Max[pts[[All,2]]]-Min[pts[[All,2]]];
tendonWidth=0.006;

scaleFactor=tendonWidth/actualWidth;
scaling=ScalingTransform[{scaleFactor,scaleFactor,scaleFactor}]; 
pts=scaling/@pts;

pts=RotationTransform[Pi/2,{1,0,0}]/@pts;
pts=RotationTransform[Pi/2,{0,0,1}]/@pts;
pts=RotationTransform[Pi/4,{1,0,0}]/@pts;

tol=scaleFactor*0.5;
markerCoords=scaleFactor*{2.1,25.1,50.9,75.1,100.1,120.1};
filtered=Table[pts[[Flatten[Position[Abs[pts[[All,1]]-markerCoords[[i]]],_?(#<=tol&)]]]],{i,1,Length[markerCoords]}];

Graphics3D[{Point[pts],Red,Point[#]&/@filtered},Axes->True,AxesLabel->{"X","Y","Z"}]

From Imaging to a Computational Mesh

Before creating the loft, we need to clean up the data. First, let’s ensure that each slice is precisely aligned along the loft’s direction (x coordinate):

Do[
filtered[[i,All,1]]=Mean[filtered[[i,All,1]]];,
{i,1,Length[filtered]}
]

As the second step, we need to ensure that each slice forms a closed loop, with points ordered sequentially to maintain a consistent structure:

ordered={};
Do[
	current=filtered[[i]];
	tmp=current[[FindShortestTour[current][[2]]]];
	center=Mean[tmp];
	angles=ArcTan@@(Transpose[(#[[{2,3}]]-center[[{2,3}]])&/@current]);
	clockwisePoints=current[[Ordering[angles,All,Greater]]];
	AppendTo[ordered,clockwisePoints];
,{i,1,Length[filtered]}]

Now, we aim to extract a reduced set of points to simplify the final mesh while maintaining an illustrative example. To achieve this, we will select a small number of marker points, ensuring a relatively simplified yet structured representation.

Additionally, we must ensure alignment of extracted points across all slices. To do this, we:

Start from the first slice and select points based on uniform angular distribution.

Translate the slice to the origin, which can facilitate alignment and ensure consistency.

Use a dynamic visualization element to interactively inspect and refine point selection.

This approach helps maintain geometric coherence while effectively reducing computational complexity:

extractionIndexes={1,6,14,25,35,40,46,52,58,68,80,90,100} ;

sliceImgs=Table[
pts0=(#- (Mean[ordered[[1,All,2;;3]]]))&/@ordered[[1,All,2;;3]];
ptsI=(#- Mean[ordered[[i,All,2;;3]]])&/@ordered[[i,All,2;;3]];

GraphicsGrid[
{{Text["Original"],Text["Aligned"]},
{
Graphics[{Black,Point[c1=Mean[ordered[[i,All,2;;3]]]],
Table[Line[{c1,ordered[[1,All,2;;3]][[extractionIndexes[[l]]]]}],{l,1,Length[extractionIndexes]}],
Red,Point[ordered[[1,All,2;;3]]],Green,Point[ordered[[1,All,2;;3]][[extractionIndexes]]],Blue,Point[ordered[[i,All,2;;3]]]}]

,Graphics[{Black,Point[c0=Mean[pts0]],
Table[Line[{c0,pts0[[extractionIndexes[[l]]]]}],{l,1,Length[extractionIndexes]}],
Red,Point[pts0],Green,Point[pts0[[extractionIndexes]]],Blue,Point[ptsI]}]},
{SwatchLegend[{Red,Green,Blue,Black},{"Slice 1","Extracted","Slice "<>ToString[i],"Angles"}]
}}]
,{i,2,Length[ordered],1}
];
ListAnimate[sliceImgs]

This allows us to extract points that are aligned relative to the first slice:

refined={};
AppendTo[refined,ordered[[1]][[extractionIndexes]]];
extracted0=ordered[[1]][[extractionIndexes]][[All,2;;3]];
centered0=(#-Mean[extracted0])&/@extracted0;
angles0=ArcTan@@@centered0;
Do[
current=ordered[[i]];
extractedI=current[[All,2;;3]];
centeredI=(#-Mean[extractedI])&/@extractedI;
anglesI=ArcTan@@@centeredI;
tmp=Nearest[anglesI,#]&/@angles0;
idx=Flatten[Position[anglesI,#[[1]]]&/@tmp];
AppendTo[refined,current[[idx]]];
,{i,2,Length[ordered]}];

Now, we have two sets of lines forming a grid that covers the tendon. This structured representation can be used to generate a mesh using OpenCascade, ensuring a well-defined and smooth surface reconstruction:

lines=Table[Line[refined[[i]]],{i,1,Length[refined]}];

Graphics3D[{Point[pts],Red,lines,Green,Table[Line[refined[[All,i]]],{i,1,Length[refined[[1]]]}]}]

Needs["OpenCascadeLink`"]
Needs["NDSolve`FEM`"]
(*Needs["MUMPSLink`"]*)

Let’s now ensure that each line forms a closed loop, maintaining continuity and consistency across the structure. This step is essential for generating a well-defined surface mesh:

refined=Append[#,First[#]]&/@refined;
lines=Line[#]&/@refined;

We can proceed to generate the lofted surface using OpenCascadeShapeLoft:

shapes=OpenCascadeShape/@lines;
loft=OpenCascadeShapeLoft[shapes,"BuildSolid"->True];

Finally, we can mesh the lofted surface into a solid tetrahedral mesh, ensuring a volumetric representation suitable for FEA. This step converts the structured loft into a watertight 3D model, enabling accurate biomechanical simulations and structural assessments:

bmesh=OpenCascadeShapeSurfaceMeshToBoundaryMesh[loft];
mesh=ToElementMesh[bmesh,"MaxCellMeasure"->(10^(-10))]

Show[mesh["Wireframe"["MeshElementStyle"->Directive[EdgeForm[],FaceForm[RGBColor[0.76,0.74,0.69]]]]]]

Finite Element Analysis

We aim to compute the stress and strain distribution within the material by simulating the mechanical stretching of a tendon. In this scenario, one end is fixed, while a load is applied to the opposite end. Physiologically, a tendon can experience forces exceeding 10 times the body weight.

The framework for structural FEA makes use of the SolidMechanicsPDEComponent, which is well described in the Solid Mechanics monograph. Let’s extract the coordinates of the two boundary sides, along with the load direction, which can be reasonably assumed to align with the tendon’s principal axis:

base=Min[mesh["Coordinates"][[All,1]]];
end=Max[mesh["Coordinates"][[All,1]]];

direction=Normalize[Mean[refined[[-1]]]-Mean[refined[[1]]]];

endPoints=mesh["Coordinates"][[Flatten@Position[mesh["Coordinates"][[All,1]],end]]];
basePoints=mesh["Coordinates"][[Flatten@Position[mesh["Coordinates"][[All,1]],base]]];

GraphicsRow[{Show[
{Graphics3D[{Red,PointSize[0.03],Arrow[{#,#+0.005*direction}]&/@endPoints[[1;;-1;;2]],Blue,Point[basePoints]}],mesh["Wireframe"["MeshElementStyle"->Directive[FaceForm[LightBrown],EdgeForm[Brown]]]]}],
SwatchLegend[{Red,Brown,Blue},{"Load","Tendon","Fixed"}]
}]

The typical loading scenario for the tendon involves a tensile load applied along its longitudinal direction, which can be represented by fixing an extremity and pulling the other one.

Tendons exhibit relatively high stiffness, which is often modeled using a Neo-Hookean hyperelastic model with a shear modulus with a typical range of 0.1–1 GPa:

μ=0.5*(10^9);
load=(Quantity[9000, "N", newtons, "Newtons"]);
force=load*direction;

pars=<|"SolidMechanicsMaterialModel"->"NeoHookean","ShearModulus"->μ|>;

The mechanical problem can be conveniently formulated using the SolidMechanicsPDEComponent function. Given the high nonlinearity of the problem and the large strains involved, we set up a parametric solver that depends on a load factor k (ranging from 0 to 1). This allows for a gradual increase in the applied load, ensuring better convergence while solving the problem:

vars={{u[x,y,z],v[x,y,z],w[x,y,z]},{x,y,z}};

pde={SolidMechanicsPDEComponent[vars,pars]==SolidBoundaryLoadValue[x==end,vars,pars,<|"Force"->k*force|>],SolidFixedCondition[x==base,vars,pars]};

numericalMethod={"PDEDiscretization"->{"FiniteElement","PDESolveOptions"->{"LinearSolver"->"MUMPS"}}};

pfun=ParametricNDSolveValue[pde,vars[[1]],{x,y,z}∈mesh,k,Method->numericalMethod]

To minimize the nonlinearities of the problem, a parametric solver can be configured to gradually increase the load. This approach is thoroughly explained in the Hyperelasticity monograph:

solver[pfun_,vars_,nsteps_]:=
Module[{displ,k,time},
displ=ConstantArray[0,Length[vars[[1]]]];
         time=AbsoluteTiming[Monitor[
Do[
k=(((step/nsteps))^2);	displ=pfun[k,"InitialSeeding"->Thread[Equal[vars[[1]], displ]]];
,{step, 1, nsteps, 1}],
step];];
Print["Solved in ", time[[1]]," s"];
displ
]

displacement=solver[pfun,vars,20]

As shown in the following figure, the tendon is stretched along its principal direction, exhibiting a displacement with a stretch ratio of up to 200%:

strain=SolidMechanicsStrain[vars,pars,displacement];
equivStrain=PDEModels`EquivalentStrain[vars,pars,strain];

VectorDisplacementPlot3D[{displacement,equivStrain},{x,y,z}∈MeshRegion[mesh],VectorSizes->Full,Boxed->False,Axes->False]

This image illustrates the final deformed configuration resulting from the applied load. The gray-shaded geometry represents the undeformed tendon. Additionally, the deformed geometry is color-coded based on the equivalent strain, which is highest in the central region and gradually decreases toward the boundaries:

Print["Maximum stretch of ",1+Max[Head[equivStrain]["ValuesOnGrid"]]," times the sample size"]

But what happens in the case of pathology? How does the mechanical response change when structural integrity is compromised?

Before investigating pathological conditions, we must account for the fact that tendons are highly transversely isotropic, with a significant load-bearing contribution from collagen fibers (which constitute up to 80% of the tendon structure) aligned along the loading direction.

Typically, this collagen fiber network exhibits a stiffness up to 20 times higher than the surrounding matrix, which is composed of elastin, cartilage, proteoglycans, inorganic components and other extracellular matrix elements. This strong anisotropic behavior plays a crucial role in the tendon's mechanical response and load distribution:

PacletInstall["WolframChemistry/ProteinVisualization"];
Needs["WolframChemistry`ProteinVisualization`"];

GraphicsRow[{ProteinBackboneAtomPlot["1BKV",SphericalRegion->True],BiomoleculePlot3D["1BKV",PlotLegends->False,ViewVertical->{0, 0, 1}]}]

Tendons are primarily composed of type I collagen fibers, which are highly organized and aligned along the loading direction. This fiber-reinforced structure enables tendons to efficiently transmit forces from muscles to bones while withstanding high mechanical loads—often exceeding 10 times body weight during intense physical activity.

The collagen matrix works in conjunction with proteoglycans, elastin and other extracellular components, contributing to the tendon's viscoelastic behavior, allowing it to store and release energy efficiently. In pathological conditions, such as tendinopathy, collagen degradation and disorganization can lead to weakened mechanical properties, increasing the risk of injury or rupture:

StreamPlot3D[direction,{x,y,z}∈mesh,StreamColorFunction->None,StreamMarkers->"Arrow",StreamScale->0.08,AspectRatio->Automatic,BoxRatios->Automatic,StreamColorFunctionScaling->False,StreamPoints->Automatic,RegionBoundaryStyle->Gray,Boxed->False,Axes->False,PlotLabel->"Load bearing fibers alignment"]

Geometry Simplification

Before further investigating the influence of collagen fibers, let’s simplify the geometry to enhance computational efficiency. We can achieve this by extracting a 2D representation of the tendon by taking a cross-sectional slice from the 3D model:

halfSpaceClip[reg_,___]/;!MeshRegionQ[reg]&&!BoundaryMeshRegionQ[reg]&&RegionEmbeddingDimension[reg]!=3=$Failed;
halfSpaceClip[mr_,h_HalfSpace]/;RegionWithin[h,mr]:=mr
halfSpaceClip[mr_,HalfSpace[n_,p_]]/;RegionWithin[HalfSpace[-n,p],mr]=EmptyRegion[3];
halfSpaceClip[mr_,HalfSpace[n_,p_]]:=Block[{rt,rot,bds,clip},rt=RotationTransform[{n,{0,0,-1}},p];
rot=TransformedRegion[mr,rt];
bds=RegionBounds[rot];
clip=DiscretizeRegion[rot,{#1+5{-1,1},#2+5{-1,1},{p[[3]],#3[[2]]+5}}]&@@bds;
InverseTransformedRegion[clip,rt]]
halfSpaceClip[___]=$Failed;

mm=MeshRegion[MeshOrderAlteration[mesh,1]];

dir=Cross[direction,{0,1,0.3}];
hs=HalfSpace[dir,Mean[mesh["Coordinates"]]];
Show[Graphics3D[{Opacity[0.8],hs}],mm,Boxed->False]

The figure displays the 3D structure of the tendon along with a box highlighting its bottom half. This setup allows for a precise mid-section slice, enabling a clearer view of the internal structure:

cut=halfSpaceClip[RegionBoundary[mm],hs];
holes=FindMeshDefects[cut,"HoleEdges","Cell"]["HoleEdges"];
slice=MeshRegion[MeshCoordinates[cut],holes,PlotTheme->"Lines",MeshCellStyle->{1->Black}];

Note that the intersection points are not evenly distributed, which may affect mesh quality. To ensure the best possible mesh, we need to remove points that are too close to each other, preventing overly small elements and improving numerical stability:

coords=MeshCoordinates[slice][[All,1;;2]];

thickness = First[Differences[MinMax[mesh["Coordinates"][[All,3]]]]];
threshold=0.09*thickness;
distance[i_,j_]:=EuclideanDistance[coords[[i]],coords[[j]]]

Before meshing, it is advisable to standardize the points. Specifically, by applying a threshold, we can filter out nearby points, ensuring a higher-quality mesh:

lines=MeshCells[slice,1];
newlines={First[lines[[1]]][[1]]};
skip=False;
currentID=1;
Do[
line=First[lines[[i]]][[2]];
If[distance[newlines[[currentID]],line]>=threshold,
	AppendTo[newlines,line];
	currentID+=1;
];
,{i,1,Length[lines]-1}
]
AppendTo[newlines,newlines[[1]]];

Graphics[{Line[coords[[#]]&/@newlines],
Point[coords[[newlines]]]}]

p=Polygon[coords[[newlines]]];
pulledSideLength=Max[coords[[All,1]]];
pulledSidePoints=Select[MeshCoordinates[p],#[[1]]>=pulledSideLength-(10^(-5))&];
pulledSidePoints=SortBy[pulledSidePoints,#[[2]]&];

clampSidePoints={#[[1]]1.1,#[[2]]}&/@pulledSidePoints;

To better apply the pulling force, let’s consider a general uniaxial strain experiment, where the pulled side is clamped to a rigid support (i.e. a material with significantly higher stiffness compared to the tendon). This setup ensures a well-defined boundary condition, mimicking realistic loading conditions. Geometrically this can be identified as a triangle. So we can mesh the geometry:

un=RegionUnion[p,MeshRegion[Line[Join[{pulledSidePoints[[1]]},clampSidePoints,{pulledSidePoints[[-1]]}]]]];
(bmesh2D=ToBoundaryMesh[un])["Wireframe"]

mesh2D=ToElementMesh[bmesh2D,"MaxCellMeasure"->0.0000008,"MaxBoundaryCellMeasure"->1,"RegionMarker"->{{p1={0.001,0.006},1},{p2={0.029,0.006},2}}]
mesh2D["Wireframe"["MeshElementStyle"->FaceForm/@{Red,Gray}]]
meshNoClamp=ToElementMesh[p,"MaxCellMeasure"->0.0000008];

As a result, we obtain a planar mesh with good quality. The mesh consists of two distinct domains: the tendon (red) and the pulling clamp (gray), each assigned different material properties.

Collagen Fibers’ Influence

Tendons can be viewed as fiber-reinforced materials, where the primary load-bearing components are the fibers, which exhibit significantly higher stiffness compared to the surrounding matrix:

μmatrix=0.1*(10^5);
stiffnessRatio = 20;
pulled=Max[mesh2D["Coordinates"][[All,1]]];
force2D=load*thickness*Normalize[direction[[1;;2]]];

pars2D=<|"SolidMechanicsModelForm"->"PlaneStress","SolidMechanicsMaterialModel"->"NeoHookean","ShearModulus"->If[ElementMarker==1,μmatrix,μmatrix*(10^5)],"Thickness"->thickness|>;
pars2D["ReinforcementStiffnessFactor"]=stiffnessRatio;
pars2D["ReinforcementDirectionVector"]=direction[[1;;2]];

vars2D={{u[x,y],v[x,y]},{x,y}};

pde2D={SolidMechanicsPDEComponent[vars2D,pars2D]==SolidBoundaryLoadValue[x>=pulled-(10^(-3)),vars2D,pars2D,<|"Force"->k*force2D|>],SolidFixedCondition[x==base,vars2D,pars2D]};

pfun=ParametricNDSolveValue[pde2D,vars2D[[1]],{x,y}∈mesh2D,k,Method->numericalMethod]

displ2D=solver[pfun,vars2D,10];

strain2D=SolidMechanicsStrain[vars2D,pars2D,displ2D];
stress2D=SolidMechanicsStress[vars2D,pars2D,strain2D,displ2D];

vonMisesStress = VonMisesStress[vars2D, parsPath, stress2D];

VectorDisplacementPlot[{displ2D,vonMisesStress},{x,y}∈MeshRegion[meshNoClamp],VectorSizes->Full,Frame->False]

This figure displays the deformed mesh, color-coded based on first principal stress. As observed, stress is particularly high in the central bottom region and gradually decreases towards the extremities. This distribution is influenced by both the thinning of the central region and the higher curvature in the bottom area, which, under the horizontal pulling load, experiences significant mechanical stress.

Pathology and Collagen Damage

From a mechanical perspective, tendon pathology is often characterized by a reduced ability of collagen fibers (the main load-bearing component) to transmit forces effectively. This degradation can lead to conditions such as tendinosis, tendinitis or even rupture, depending on the severity of the lesion, as detailed in studies from Arya and Kulig, 2010; Yin et al., 2021; and Freedman et al., 2014:

coords=MeshCoordinates[tendonMesh];
cells=MeshCells[tendonMesh,2];
validCells1=Select[cells,-200<=Max[coords[[#[[1]]]][[All,3]]]<=20&];
validCells2=Select[cells,60<=Max[coords[[#[[1]]]][[All,3]]]<=200&];
validCells3=Select[cells,20<=Max[coords[[#[[1]]]][[All,3]]]<=60&];
clippedMesh1=MeshRegion[coords,validCells1];
clippedMesh2=MeshRegion[coords,validCells2];
internalMesh=MeshRegion[coords,validCells3];
m[x_]:=AnatomyData[x,"Graphics3D"][[1]]
muscles=Join[AnatomyData[(Entity["right calcaneal tendon", Entity["AnatomicalStructure","RightCalcanealTendon"], "Entity[\"AnatomicalStructure\", \"RightCalcanealTendon\"]", "anatomical structure"])]["AttachedMuscles"],AnatomyData[(Entity["right calcaneal tendon", Entity["AnatomicalStructure","RightCalcanealTendon"], "Entity[\"AnatomicalStructure\", \"RightCalcanealTendon\"]", "anatomical structure"])]["InsertedMuscles"]];
G=Graphics3D[{
Opacity[0.05],m[#]&/@muscles,
m[(Entity["skeleton of right free lower limb", Entity["AnatomicalStructure","SkeletonOfRightFreeLowerLimb"], "Entity[\"AnatomicalStructure\", \"SkeletonOfRightFreeLowerLimb\"]", "anatomical structure"])],
m[#]&/@{(Entity["set of lumbricals of right foot", Entity["AnatomicalStructure","SetOfLumbricalsOfRightFoot"], "Entity[\"AnatomicalStructure\", \"SetOfLumbricalsOfRightFoot\"]", "anatomical structure"]),(Entity["right flexor accessorius", Entity["AnatomicalStructure","RightFlexorAccessorius"], "Entity[\"AnatomicalStructure\", \"RightFlexorAccessorius\"]", "anatomical structure"]),(Entity["right flexor hallucis brevis", Entity["AnatomicalStructure","RightFlexorHallucisBrevis"], "Entity[\"AnatomicalStructure\", \"RightFlexorHallucisBrevis\"]", "anatomical structure"]),(Entity["right adductor hallucis", Entity["AnatomicalStructure","RightAdductorHallucis"], "Entity[\"AnatomicalStructure\", \"RightAdductorHallucis\"]", "anatomical structure"]),(Entity["flexor digiti minimi brevis of right foot", Entity["AnatomicalStructure","FlexorDigitiMinimiBrevisOfRightFoot"], "Entity[\"AnatomicalStructure\", \"FlexorDigitiMinimiBrevisOfRightFoot\"]", "anatomical structure"]),(Entity["plantar interosseous of right foot", Entity["AnatomicalStructure","PlantarInterosseousOfRightFoot"], "Entity[\"AnatomicalStructure\", \"PlantarInterosseousOfRightFoot\"]", "anatomical structure"])},
Opacity[0.05],m[#]&/@{(Entity["skin of right foot", Entity["AnatomicalStructure","SkinOfRightFoot"], "Entity[\"AnatomicalStructure\", \"SkinOfRightFoot\"]", "anatomical structure"]),(Entity["skin of right leg", Entity["AnatomicalStructure","SkinOfRightLeg"], "Entity[\"AnatomicalStructure\", \"SkinOfRightLeg\"]", "anatomical structure"])}},Axes->False,PlotRange->{All,All,{All,200}},Boxed->False,ViewPoint->{-100,100,30}];
GraphicsGrid[{
{Text["Healthy"],Text["Tendinosis"],Text["Rupture"]},
{Show[G,Graphics3D[{FaceForm[Brown],EdgeForm[],clippedMesh1,clippedMesh2,internalMesh},Boxed->False]],
Show[G,Graphics3D[{FaceForm[Brown],EdgeForm[],clippedMesh1,clippedMesh2,FaceForm[Red],EdgeForm[LightRed],internalMesh},Boxed->False]],
Show[G,Graphics3D[{FaceForm[Brown],EdgeForm[],clippedMesh1,clippedMesh2},Boxed->False]]
}}]

Mathematically, this can be represented as a localized alteration of the mechanical properties within a specific region of the tendon. This may include changes in stiffness, elasticity or failure thresholds, affecting the tendon’s overall ability to sustain loads:

c={First[Min[mesh["Coordinates"][[All,1]]]+Differences[MinMax[mesh["Coordinates"][[All,1]]]]/2],First[Min[mesh["Coordinates"][[All,2]]]+Differences[MinMax[mesh["Coordinates"][[All,2]]]]/5],
First[Min[mesh["Coordinates"][[All,3]]]+Differences[MinMax[mesh["Coordinates"][[All,3]]]]/2]};
r =2First[Differences[MinMax[mesh2D["Coordinates"][[All,2]]]]/5];

damage[x_,y_,z_]:=If[(((x-c[[1]]))^2)+(((y-c[[2]]))^2)+(((z-c[[3]]))^2)<=(r^2),0.05+0.95*(((((Norm[c-{x,y,z}])/r)))^3),1]
damage[x_,y_]:=damage[x,y,c[[3]]]

DensityPlot3D[damage[x,y,z],{x,y,z}∈mesh,PlotLegends->False,AspectRatio->Automatic,PlotRange->All,Boxed->False,BoxRatios->Automatic,Axes->False,PlotLabel->"Pathology: fibers degradation"]

As previously mentioned, fibers play a crucial role in the stiffness and load-bearing capacity of tendons due to their widespread distribution and significantly higher stiffness (up to 20 times greater than other components). Now, let’s investigate their influence on the tendon’s mechanical behavior using the simplified 2D mesh:

parsPath=pars2D;

parsPath["ReinforcementStiffnessFactor"]=20*damage[x,y];

GraphicsRow[{StreamPlot[direction[[1;;2]],{x,y}∈meshNoClamp,StreamColorFunction->(ColorData[{"Rainbow","Reverse"}][damage[#1,#2]]&),StreamMarkers->"Arrow",StreamScale->0.02,AspectRatio->Automatic,BoxRatios->Automatic,StreamColorFunctionScaling->False,StreamPoints->(10^6),Axes->False,Frame->False,RegionFillingStyle->LightBrown,RegionBoundaryStyle->Black,ImageSize->Large]
,BarLegend[{"Rainbow",{0,1}},LegendLabel->"Damage"]}]

The plot illustrates the undeformed planar section of the tendon along with the fiber distribution, which is predominantly aligned along the longitudinal direction. The rainbow colors represent the damage properties, with the highest concentration (red) located in the central bottom region. This directly correlates with a reduction in mechanical properties, affecting the tendon’s structural integrity:

pdePath={SolidMechanicsPDEComponent[vars2D,parsPath]==SolidBoundaryLoadValue[x>=pulled-(10^(-3)),vars2D,parsPath,<|"Force"->k*force2D|>],SolidFixedCondition[x==base,vars2D,parsPath]};

pfun=ParametricNDSolveValue[pdePath,vars2D[[1]],{x,y}∈mesh2D,k,Method->numericalMethod]

displPathology=solver[pfun,vars2D,10];

GraphicsRow[{Show[
mesh2D["Wireframe"["MeshElementStyle"->Directive[FaceForm[LightGray],EdgeForm[LightGray]]]],
def=ElementMeshDeformation[mesh2D,displPathology,"ScalingFactor"->1]["Wireframe"["MeshElement"->"BoundaryElements","MeshElementStyle"->Red]],
ElementMeshDeformation[mesh2D,displ2D,"ScalingFactor"->1]["Wireframe"["MeshElement"->"BoundaryElements","MeshElementStyle"->Black]]
]
,SwatchLegend[{Gray,Red,Black},{"Undeformed","Pathology","Healthy"}]}]

The relaxation of pathological fibers is clearly observed in the displacement field, with the pulled side exhibiting noticeable rotation. But how does this affect strain and stress distribution?

Both stress and strain play a critical role in injury risk and the progression of inflammation, making their evaluation essential for understanding the mechanical implications of tendon pathology:

strainPath=SolidMechanicsStrain[vars2D,parsPath,displPathology];
stressPath=SolidMechanicsStress[vars2D,parsPath,strainPath,displPathology];

vonMisesStressPath = VonMisesStress[vars2D, parsPath, stressPath];

VectorDisplacementPlot[{displPathology,strainPath[[1,1]]},{x,y}∈p,VectorSizes->Full,Frame->False, PlotLabel->"Deformation colored by uniaxial strain",ColorFunction->(ColorData["Rainbow"][Rescale[#5,{0,20}]]&),ColorFunctionScaling->False]

GraphicsRow[{DensityPlot[vonMisesStress,{x,y}∈meshNoClamp,AspectRatio->Automatic,PlotLegends->False,PlotRange->{All,All,{0,(10^7)}},ColorFunction->(ColorData["Rainbow"][Rescale[#,{0,3*(10^6)}]]&),ColorFunctionScaling->False,Frame->False,PlotLabel->"Healthy",BoundaryStyle->Black],
DensityPlot[vonMisesStressPath,{x,y}∈meshNoClamp,AspectRatio->Automatic,PlotLegends->False,PlotRange->{All,All,{0,(10^7)}},ColorFunction->(ColorData["Rainbow"][Rescale[#,{0,3*(10^6)}]]&),ColorFunctionScaling->False,Frame->False,PlotLabel->"Pathologic",BoundaryStyle->Black],
BarLegend[{"Rainbow",{0,2}},LegendLabel->"σ [MPa]"]
}]

This figure illustrates how, due to the pathology, the stress is significantly higher in the bottom region. Additionally, the final deformation on the right side is more pronounced compared to the healthy case. Furthermore, it allows for a direct comparison of stress distribution between the healthy and pathological cases, highlighting the impact of the condition on the tendon’s mechanical behavior:

avσ=NIntegrate[vonMisesStressPath,{x,y}∈mesh2D];
DensityPlot[(vonMisesStressPath-vonMisesStress)/avσ,{x,y}∈meshNoClamp,AspectRatio->Automatic,PlotRange->{All,All,{-35,35}},PlotLegends->Automatic,ColorFunction->(ColorData["Rainbow"][Rescale[#,{-20,20}]]&),ColorFunctionScaling->False,Frame->False,BoundaryStyle->Black,PlotLabel->"Relative stress increase"]

As we can see, stress slightly increases in the upper region, while it slightly decreases on the bottom side. Note also the increase near the damaged region. These effects warrant further investigation, particularly in relation to pathology progression.

Curious about the technical details? You can find all the relevant information on FEM and structural modeling on the PDEModels Overview reference page.

However, for a deeper analysis, you’ll have to stay tuned for the next post!

References for BML

Viceconti, M., & Emili, L. - Toward good simulation practice: Best practices for the use of computational modelling and simulation in the regulatory process of biomedical products. Cham: Springer Nature Switzerland. https://doi.org/10.1007/978-3-031-48284-7

Gray, Henry. 2015. Gray’s Anatomy: With Original Illustrations by Henry Carter. London, England: Arcturus Publishing Ltd.

Khayyeri H., Longo G., Gustafsson A., Isaksson H. - Comparison of structural anisotropic soft tissue models for simulating Achilles tendon tensile behaviour. 
Journal of the Mechanical Behavior of Biomedical Materials, Vol 61, 2016, https://doi.org/10.1016/j.jmbbm.2016.04.007.

Bajuri, M.N., Isaksson, H., Eliasson, P. et al. A hyperelastic fibre-reinforced continuum model of healing tendons with distributed collagen fibre orientations. Biomech Model Mechanobiol 15, 1457–1466 (2016). https://doi.org/10.1007/s10237-016-0774-5

Shruti A., Kornelia K. - Tendinopathy alters mechanical and material properties of the Achilles tendon - Journal of Applied Physiology 2010 108:3, 670-675 https://doi.org/10.1152/japplphysiol.00259.2009

Freedman BR, Gordon JA, Soslowsky LJ. The Achilles tendon: fundamental properties and mechanisms governing healing. Muscles Ligaments Tendons J. 2014 Jul 14;4(2):245-55. PMID: 25332943; PMCID: PMC4187594.

Nai-Hao Yin Paul Fromme Ian McCarthy Helen L Birch (2021) Individual variation in Achilles tendon morphology and geometry changes susceptibility to injury. https://doi.org/10.7554/eLife.63204