Probabilistic finite element analysis of a craniofacial finite element model

12 Pages·2012·1.31 MB·English

by  Michael A. Berthaume

JournalofTheoreticalBiology300(2012)242–253
ContentslistsavailableatSciVerseScienceDirect
Journal of Theoretical Biology
journal homepage: www.elsevier.com/locate/yjtbi
Probabilistic finite element analysis of a craniofacial finite element model
Michael A. Berthaumea,b, Paul C. Dechowc, Jose Iriarte-Diazd, Callum F. Rossd, David S. Straite,
Qian Wangf, Ian R. Grossea,n
aDepartmentofMechanical&IndustrialEngineering,UniversityofMassachusetts,Amherst,USA
bDepartmentofAnthropology,UniversityofMassachusetts,Amherst,USA
cDepartmentofBiomedicalSciences,TexasA&MHealthScienceCenterBaylorCollegeofDentistry,Dallas,TX,USA
dOrganismalBiologyandAnatomy,UniversityofChicago,Chicago,IL,USA
eDepartmentofAnthropology,UniversityatAlbany,Albany,NY,USA
fDivisionofBasicMedicalSciences,MercerUniversitySchoolofMedicine,Macon,GA,USA
a r t i c l e i n f o a b s t r a c t
Articlehistory: We employed a probabilistic finite element analysis (FEA) method to determine how variability in
Received7July2011 material property values affects stress and strain values in a finite model of a Macaca fascicularis
Receivedinrevisedform cranium.Thematerialbehaviorofcorticalbonevariedinthreeways:isotropichomogeneous,isotropic
1December2011 non-homogeneous,andorthotropicnon-homogeneous.Thematerialbehaviorofthetrabecularbone
Accepted18January2012
and teeth was always treated as isotropic and homogeneous. All material property values for the
Availableonline27January2012
craniumwererandomizedwithaGaussiandistributionwitheithercoefficientsofvariation(CVs)of
Keywords: 0.2orwithCVscalculatedfromempiricaldata.Latinhypercubesamplingwasusedtodeterminethe
Probabilisticanalysis valuesofthematerialpropertiesusedinthefiniteelementmodels.Intotal,fourhundredandtwenty
Finiteelementanalysis
sixseparatedeterministicFEsimulationswereexecuted.
Modelvalidation
Wetestedfourhypothesesinthisstudy:(1)uncertaintyinmaterialpropertyvalueswillhavean
Isotropicandorthotropicmaterial
insignificanteffectonhighstressesandasignificanteffectonhighstrainsforhomogeneousisotropic
properties
Macacafascicularis models;(2)theeffectofvariabilityinmaterialpropertyvaluesonthestressstatewillincreaseasnon-
homogeneityandanisotropyincrease;(3)variationintheinvivoshearstrainvaluesreportedbyStrait
et al. (2005) and Ross et al. (2011) is not only due to variations in muscle forces and cranial
morphology,but also due to variation in material property values; (4) the assumption of a uniform
coefficientofvariationforthematerialpropertyvalueswillresultinthesametrendinhowmoderate-
to-highstressesandmoderate-to-highstrainsvarywithrespecttothedegreeofnon-homogeneityand
anisotropy as the trend found when the coefficients of variation for material property values are
calculatedfromempiricaldata.Ourresultssupportedthefirstthreehypothesesandfalsifiedthefourth.
When material properties were varied with a constant CV, as non-homogeneity and anisotropy
increasedthelevelofvariabilityinthemoderate-to-highstrainsdecreasedwhilethelevelofvariability
inthemoderate-to-highstressesincreased.However,thisisnotthepatternobservedwhenCVscalculated
fromempiricaldatawereappliedtothematerialpropertieswherethelowestlevelofvariabilityinboth
stressesandstrainsoccurredwhenthecraniumwasmodeledwithalowlevelofnon-homogeneityand
anisotropy.Therefore, whenconstant materialproperty variability is assumed, inaccurate trends inthe
levelofvariabilitypresentinmodest-to-highmagnitudestressesandstrainsareproduced.
When the cranium is modeled with the highest level of accuracy (high non-homogeneity and
anisotropy)andwhenrandomnessinthematerialpropertiesiscalculatedfromempiricaldata,thereisa
largelevelofvariabilityinthesignificantstrains(CV¼0.369)andalowlevelofvariabilityinthemodest-to-
high magnitude stresses (CV¼0.150). This result may have important implications with regard to the
mechanicalsignalsdrivingboneremodelingandadaptationthroughnaturalselection.
&2012ElsevierLtd.Allrightsreserved.
1. Introduction
1.1. Background:FEAinbiology
Over the past ten years many researchers in comparative
nCorrespondingauthor.Tel.:þ14135451350. morphology and physical anthropology have adopted finite ele-
E-mailaddress:[e-mail](I.R.Grosse). ment analysis (FEA) as a tool to simulate feeding biomechanics
0022-5193/$-seefrontmatter&2012ElsevierLtd.Allrightsreserved.
doi:10.1016/j.jtbi.2012.01.031
M.A.Berthaumeetal./JournalofTheoreticalBiology300(2012)242–253 243
(Straitetal.,2005;Rossetal.,2011;Dumontetal.,2011;Wood orthogonal to each other (Malvern, 1969). Spatial non-homoge-
etal.,2011;O’Higginsetal.,2011;Porroetal.,2011;Davisetal., neityrequiresspatiallydetailedmeasuresofmaterialproperties,
2011; Berthaume et al., 2010; Rayfield, 2007; Dumont et al., which can be accounted for by assigning different material
2005). Finite element models of craniofacial skeletal structures properties to separate regions in a FEM (Strait et al., 2005;
enableresearcherstopredictthecompletestressandstrainfields Porroetal.,2011;WangandDechow,2006).
throughout the cranium as a result of loading conditions that Anisotropy and spatial non-homogeneity in bone material
simulate orofacial functions, such as mastication or ingestion of properties can produce an increase in spatial variation in stress
mechanically resistant foods, and to test hypotheses regarding and strain values across a single model. When combined with
craniofacial shape, feeding performance, and behavior (Strait variation in material properties of other connective tissues (e.g.,
etal.,2005;Rossetal.,2011;Rossetal.,2005). sutures (Reed et al., 2011; Wang et al., 2010; Rayfield, 2005)),
FEAisacomputer-basedtechniqueusedtosolvethegovern- constraining conditions, and external forces, variation in bone
ing field equations in continuum mechanics, such as the equa- material properties can produce significant variation in model
tions of elasticity. It is widely used by engineers for structural behavior(Porroetal.,2011).Theimplicationsofthisvariationfor
analysis, thermal analysis, fluid flow analysis, and electromag- our understanding of biological morphology through natural
netic analysis. Initial applications of FEA to biological systems selectionareonlyjustbeginningtobeexplored.
involved relatively simple geometric idealizations of human
skeletal structures, such as the femur (Brekelmans et al., 1972; 1.2. Viewingmaterialpropertiesas‘‘random’’variables
Rybicki et al., 1972). Subsequent improvements in technology
enabled geometrically complex tissue structures to be digitally A parameter whose values are not known precisely due to
reconstructedfromCTscans,openingthedoortotheapplication uncertaintyinmeasurementand/orotherfactorsisconsideredto
ofFEAtoawidevarietyofbiologicalsystems,suchasmammalian bearandomorstochasticvariable.Ifthevariableiscontinuous,it
skulls (Dumont et al., 2005), dinosaur skulls (Rayfield, 2004), may be mathematically modeled by a probability density func-
alligatorskulls(Porroetal.,2011;Metzgeretal.,2005;Reedetal., tion, which mathematically defines the probability that the
2011), shark jaws (Wroe et al., 2008), seahorse snouts (Leysen variablewill have a specific value. It should be emphasized that
etal.,2010),andbluecrablegs(Hechtetal.,2010). randomness in a variable does not require that the variable has
One of the few finite element models of a biological system equal probability to take on any possible value between two
validated by in vivo data is the macaque cranium model devel- extremes. Such a random variable is defined by a uniform
oped by Strait et al. (2005) and Ross et al. (2005) to simulate distribution—which would be a rather poor representation of a
feeding biomechanics. (This model was validated again with material property value. Instead, distributions such as Gaussian,
additional in vivo data gathered from multiple individuals in lognormal,orWeibulldistributionhavebeenusedtomodelthe
Ross et al. (2011). Other FE models validated with in vivo data natural variability observed in mechanical property values for
include a rat ulna (Kotha et al., 2004) and an alligator cranium bothbiologicalandengineeredmaterials(McElhaneyetal.,1970;
(Metzgeretal.,2005).)Itshouldbenotedthattherewasalarge Huesteetal.,2004).Comparatively,engineeredmaterialstendto
levelofvariationintheinvivodatathatwasusedtovalidatethe have small coefficients of variation due to high levels of control
macaque cranium model, and understanding from where that during the production process whereas biological materials can
levelofvariationoriginatedisofparticularinterestinthisstudy. havefairlylargecoefficientsofvariationsforthereasonsalluded
In their FE model the researchers varied the type of material toearlier.
behavior (isotropic versus orthotropic) and the degree to which We also note that randomness does not necessarily imply a
material property values vary spatially, i.e. non-homogeneity. lackoforderorstructure.Twoormorerandomvariablesmaybe
Thiswasassessedbydividingthecraniumintonumerousregions, correlated,andthiscorrelationimposesameasureoforderonthe
each of which could be assigned different material property randomness. For example, consider the Youngs moduli of an
values with different levels of anisotropy. Results from this set orthotropicmaterial: E , E , and E . These moduli are defined as
3 2 1
of models were then compared to the range of strain values the maximum, intermediate, and minimum stiffnesses, respec-
measured in vivo on different individuals and collected and tively,correspondingtothreeorthogonaldirectionscalledprinci-
publishedbyvariousresearchers.Sincemostofthedatagathered pal material directions. From experimental data, normal
fromthesimulationsfellwithinthedispersionoftheinvivodata, distributions for the moduli can be gathered and correlation
the model was considered validated and to be efficient at coefficients between them can be determined to increase the
predictingstressandstrainpatternsonthesurfaceofamacaque likelihood that E 4E 4E (Peterson and Dechow, 2003; Wang
3 2 1
craniumduringorofacialfunction. andDechow,2006).
Despite such successes, some morphologists have concerns
regarding the deterministic nature of finite element models 1.3. ProbabilisticanalysesinFEA
(Grineet al.,2010).Thisistypicallynotasmuchanissueinthe
engineering world, as considerable effort is expended by engi- Here we apply the techniques of probabilistic FEA to deter-
neers to reduce variability in the performance of engineered minethevariability(orrandomness)invariousoutputvariables
products by precise control of the product geometry, material inanFEmodelofacraniumduetorandomnessorvariabilityin
property values, and environmental loading conditions (if possi- material property values. We note that in a conventional finite
ble). However, this is not true in the biological world. Material element analysis, approximate solutions are found using differ-
properties of naturally occurring biological structures, such as ential equations that are deterministic. Deterministic inputs
bone,are inherentlyvariabledue to inter-individualvariationin include boundary conditions such as loads and constraints,
sex,age,genetics,andenvironment(PetersonandDechow,2003; materialproperties,andgeometry.Deterministicoutputsinclude
WangandDechow,2006;Zapataetal.,2010;Wangetal.,2010). quantitiessuchasdeformationofthespecimenandthestressand
These variations can be manifested through varying degrees of strainfields.Inprobabilisticfiniteelementanalysis,therearetwo
anisotropyandspatialnon-homogeneity(Melchers,1987).Aniso- types of techniques used to generate variability in outputs
tropicmaterialpropertiesofbonearetypicallyaccountedforby resulting from the variability in inputs: non-statistical methods
orthotropic material models, requiring nine independent elastic and statistical methods (Sudret and der Kiureghian, 2000; Kang,
constants and three principal material directions which are 2005). Non-statistical methods involve perturbation techniques
244 M.A.Berthaumeetal./JournalofTheoreticalBiology300(2012)242–253
by introducing analytical functional approximations of the regions and a significant effect on moderate-to-high strains for
responsesaboutthemeaninput.Themethodiscomputationally homogeneous isotropic models. Second, it is hypothesized that
efficient but requires significant mathematical extensions to theeffectofvariability inmaterialproperty values onthestress
conventional FEA algorithms and has not been fully developed statewillincreaseasnon-homogeneityandanisotropyincreasein
to encompass the breadth of problems that can currently be thematerial.Third,itishypothesizedthatifthefirsthypothesisis
solved by conventional FEA tools. As a result, it is not found in true, then a significant contribution to variation in the in vivo
mainstreamcommercialFEAtools. shearstrainvaluesreportedbyStraitetal.(2005)andRossetal.
In the statistical approach, (which includes methods such as (2011)maybeattributedtovariationinmaterialpropertyvalues
the Monte Carlo simulation), input parameters are randomized betweenindividuals.Finally,itishypothesizedthattheassump-
accordingtoprescribedprobabilisticdistributions(Gaussian,log- tion of a constant coefficient of variation for all the material
normal,etc.)andasamplingalgorithm.Eachsetofvaluesforthe propertieswillresultinthesametrendinthevariationinstresses
randominputparametersproducesasetofresults(i.e.displace- and strains with respect to the degree of non-homogeneity and
ment, strain, and stress fields) through deterministic FEA. Post- anisotropyasfoundwhenthecoefficientsofvariationformaterial
processing across all results sets yields statistics for output propertiesarecalculatedfromempiricaldata.
variables(i.e.suchasdeflectionatspecificnodes,maximumvon Thefirsttwohypothesesarebasedontheobservationthatthe
Mises stress or strain, or von Mises stress or strain at specific macaquecraniofacialskeletonconsistspredominantlyofcortical
locations in the mesh). The accuracy of output statistics can be bone (the total cranium volume was 39,359mm3, with 86.5% of
improved with increased sampling. Statistical sampling techni- the cranium volume being cortical bone, 5.8% being teeth, and
ques,suchasLatinHypercubeSampling,areoftenemployedover 7.7% being trabecular bone) and hence is similar to a single-
direct Monte Carlo Simulations, to reduce the number of deter- material structure. In the linear elastic region of material beha-
ministic FEA analyses required, and response surface modeling vior,thestressstateforahomogeneous,isotropicsingle-material
approachescanalsobeusedtoreducecomputationalcosts. structureisindependentofYoung’smodulus,butthestrainstate
The probabilistic design portion of ANSYS is an implementa- isdependentonit(TimoshenkoandGoodier,1970).Thevolume
tion of the statistical approach (Reh et al., 2006). It takes into ofthecraniumwascalculatedusinganANSYSAPDLcommand,a
account the natural variability in input variables by applying a different method for calculating volume than was used in
predefinedstatisticaldistributiontothem.Theresultscanbeused Chamoli and Wroe (2011). Also, it should be noted that what
toansweranumberofquestions,suchasiftheinputvariablesare ChamoliandWroe(2011)defineasheterogeneous,wedefineas
likely to fall within a certain range, how large the range of the three region isotropic, which we consider to have a low level of
outputparameterswillbe,andwhichinputvariablesarethemost non-homogeneity. Anisotropy and lack of homogeneity violate
stronglycorrelatedtoeach,individualoutputvariable. the conditions for independence of the stress state on Young’s
Probabilisticfiniteelementanalysisisnotnewtoengineers.Itis modulus. Thus, with increasing anisotropy and lack of homoge-
widelyusedinreliability-basedapproachesthatpredicttheprob- neity, we expect increasing variability in bone stresses due to
abilityofafailureofasystem(Melchers,1987;DerKiureghianand stochastic material property values, particularly with respect to
Ke,1988;CesareandSues,1999;Lewis,1987).However,applica- Young’s modulus. These hypotheses were tested by comparing
tions of probabilistic FEA to biological systems are limited to the theresultsoftheprobabilisticanalysestoeachother.
biomedical field. Thacker et al. (2000) used probabilistic FEA to The third hypothesis could be falsified if variations in strain
estimatethecontributiontoinjuryprobabilityof29materialsand due to variation in material property values are relatively small
geometry-related parameters of a finite element model of the comparedtothelevelofvariationinstrainobservedinvivo.This
lowercervicalspine.Nicolellaetal.(2001)alsousedprobabilistic would imply that variability due to material property values
FEA ofa false hip joint todetermine the probabilityofshear and between individuals is secondary to variations in morphological
fatigue failure given certain random variables, such as muscle structure and masticatory loads between individuals. This
loading patterns and material properties of the cortical and hypothesis was tested by comparing the results of the probabil-
trabecularbone.Oneprobabilisticanalysishasbeenusedtostudy istic analyses to the in vivo data. The fourth hypothesis will be
themacaquezygomaticarch(Kupcziketal.,2007). testedbycomparingtheresultsof3modelsthathaveaconstant
Because of the limited use of probabilistic finite element coefficientofvariationof0.2appliedtothem(models1,3,and5,
analyses on non-human biological systems, one of the primary seeTable1)withtheresultsof3modelsthathavecoefficientsof
purposes of this paper is to compare the results of probabilistic variation applied to them which were calculated from empirical
analyses to the in vivo data used to validate these models. This data(models2,4,and6,seeTable1).
comparison will help to determine whether or not variability in
materialpropertiescontributessignificantlytostrainvariability.
2. Materialsandmethods
1.4. Thehypotheses 2.1. MacaqueFEmodel
Four hypotheses are tested in this paper. First, it is hypothe- TheFEmodelofthemacaquecraniumusedinthisstudywas
sized that variability in bone material property values between constructed by Strait et al. (2005). The macaque model, which
individuals will have little effect on moderate-to-high stresses consistsof379,388elementsand131,293nodes,wasdeveloped
Table1
Descriptionofthesixmodelsconstructedforthesimulations.
Physicaldescription CV¼0.2forinput CVforinputvariablesbased Levelofnon- Levelof Numberof Numberof
ofmodel variables onempiricaldata homogeneity anisotropy randominputvariables simulationsrun
3Region,isotropic Model1 Model2 Low Low 6 30
Multipleregion,isotropic Model3 Model4 High Low 38 39
Multipleregion,orthotropic Model5 Model6 High High 143 144
M.A.Berthaumeetal./JournalofTheoreticalBiology300(2012)242–253 245
Fig.1. Locationstestedonthecraniumaredepictedinnavyblue.The34locationsarenumbered1to34.20ofthe34locationswerechosenbecausetheywerethesameas
the20sitesanalyzedinRossetal.(2011)andtheother14siteswerechosentoextendthetotalnumberoflocationsbeingtestedtotheentireanteriorportionofthe
cranium.(Forinterpretationofthereferencestocolorinthisfigurelegend,thereaderisreferredtothewebversionofthisarticle.)
in Algor, a commercial FEA program. However, the model was laterally (i.e. similar to the ‘‘wishboning’’ seen in mandibles
importedinto ANSYS APDL13.0 (seewww.ansys.com) usingthe (HylanderandJohnson,1994)).Thenodesontheocclusalsurface
NASTRANfileformatsothatANSYS’sprobabilisticandsensitivity of the left M1 were also constrained in the superior–inferior
analysis could be utilized. Subsequently, a number of modifica- directiontosimulatebitingdownonahardfooditem.
tions were made to the model, but the muscle forces remained
unchanged. For more information about the construction of the
2.2. ApplyingprobabilisticanalysistotheFEM
model,pleaseseeStraitetal.(2005).
FortheisotropicmaterialmodelstestedbyStraitetal.(2005),
Atotalofsixmodelsofthemacaquecraniumwerecreatedand
shearmodulus,G,Young’smodulus,E,andPoisson’sratio,n,were
analyzed.Eachwasassignedauniquesetof(a)materialproper-
all specified based on empirical data obtained by (Wang and
ties and/or (b) CVs to the material properties. The analyses for
Dechow, 2006). However, isotropy imposes a relationship
eachmodelconsistedofanywherebetween30and144separate
betweenthreeproperties(E,G,andn)givenbytheequation:
deterministic FE simulations using the Probabilistic Design por-
G¼E=ð2ð1þnÞÞ ð1Þ tionofANSYSAPDL13.0(seeTable1).
The Probabilistic Design module requires three inputs: a
Thus,foratrueisotropicmaterialmodelonlytwoofthethree
macrocodethatdrivesANSYS,alistofinputvariables,andalist
properties may be specified independently. In the ANSYS model
of output variables. In Probabilistic Design, the macro code is
for isotropic materials, the values of Young’s modulus and
executed N number of times, resulting in N separate finite
Poisson’s ratio published by Strait et al. (2005) were utilized
element analyses. Each deterministic simulation involves a
and the program calculated the shear modulus based on the
unique set of input parameters that have been designated as
isotropicrelationshipofEq.(1).Fororthotropicmaterialbehavior,
random input variables. Statistics for each designated output
nineelasticconstantsarerequired(E ,E ,E ,G ,G ,G ,n ,n ,
1 2 3 12 13 23 12 13 variable are calculated using the set of N results obtained
andn )whichoperateinthreeorthogonaldirections(directions
23 correspondingtotheNdeterministicsimulations.
1, 2, and 3) and on three planes (planes 12, 23, and 13), which
In our analyses the random input variables were the set of
define the principal material directions and planes within each
material properties assigned to each model. Primary output
cranialregion.
variables were the minimum and maximum principal stresses
When the model was constructed by Strait et al. (2005) the
andstrainsandthevonMisesstressesandstrains,exportedfrom
teethweremodeledaspartofthecorticalbone.Forpurposesof
thesetofnodalsolutionsat34sitesonthecranium(alldepicted
thisexperiment,theteethwereseparatedfromthecorticalbone
inFig.1).Theotheroutputvariablesweree , e , and g normal
along the gum line. Because the model was meshed before the x y xy
andshearstrainsateachsite,wherethesestrainsarebasedona
teethwereseparated,itwasnotpossibletoseparateouttheroots
local xyz coordinate system at each of the 34 sites. Each local
of the teeth or the periodontal ligament (PDL) from the
coordinatesystemhaditsoriginlocatedatthenodebeingtested
corticalbone.
andthelocalxyplanewastangenttothesurfaceofthecranium.
Tosimulatethemechanicsofbiting,anaxisofrotationwasset
Forhalfofthesimulations,materialpropertyvalues(theinput
upatthetemporomandibularjoint(TMJ)byconstrainingasingle
variables) were randomized using Gaussian distributions with
nodeatoneTMJinallthreecoordinatedirectionsandconstrain-
mean values obtained from empirical data (Wang and Dechow,
ingasinglenodeattheotherTMJagainstanterior–posteriorand
superior–inferiormotion.However,thenodewasnotconstrained
against medial lateral motion, allowing the cranium to deform 1Engineeringshearstrainistwicethetensorshearstrain.
246 M.A.Berthaumeetal./JournalofTheoreticalBiology300(2012)242–253
2006) and CVs of 0.20. In the other half of the simulations, materialproperties(Table2).Inmodel1theinputvariableswere
material property values were randomized using Gaussian dis- randomizedwithaCVof0.2andaGaussiandistribution.Inmodel
tributionswithmeanvaluesobtainedfromempiricaldata(Wang 2theinputvariableswererandomizedwiththeirexperimentally
and Dechow, 2006) and standard deviations based on empirical determined standard deviations and a Gaussian distribution
data (Wang and Dechow, 2006). A truncated Gaussian distribu- (WangandDechow,2006).
tionwasassignedtoPoisson’sratiosoitneverexceeded0.4999.
(Values larger than 0.5 violate fundamental thermodynamic
principlesandarenotallowedinthetheoryofelasticity.)
2.4.2. Models3and4:regionalisotropicmaterialproperties
AMonteCarlosimulationwithLatinhypercubesamplingwas
In models 3 and 4 the cranium was divided into 35 distinct
used as a method of choosing the values for the random inputs
regions,increasingnon-homogeneityfrommodels1and2.Each
(OlssonandSandberg,2002).WhenLatinhypercubesamplingis
region was assigned its own distinct set of isotropic material
used in conjunction with a Monte Carlo simulation, an accurate
propertieswhicharedefinedinTable3.Similartomodels1and2,
statistical distribution can be obtained while running fewer
material properties were randomized with a CV of 20% and a
simulations. (For more information about Monte Carlo simula-
Gaussian distribution in model 3 and all the material properties
tions with Latin hypercube sampling and how it is used in FEA,
wererandomizedwiththeirexperimentallydeterminedstandard
please refer to Olsson and Sandberg (2002)) The minimum
deviationsandaGaussiandistributioninmodel4.
number of simulations that needed to be run, N, were N¼nþ1
wherenisthenumberofinputvariables.Ncannotbelessthan30
inordertopermitanormaldistributionfortheoutputvariables.
2.4.3. Models5and6:regionalorthotropicmaterialproperties
Because of this, 30 simulations were used to randomize input
Formodels5and6,themacaquecraniumwasdividedintothe
parameters for models 1 and 2, 39 simulations were used in
same regions as in Models 3 and 4. However, each region was
models3and4,and144simulationswereusedinmodels5and6
assigned its own set of orthotropic material properties, thus
(see Table 1). The minimumnumber of simulationswas used in
increasing anisotropy. Exceptions to this were the teeth, trabe-
eachMonteCarlosimulationbecauserunningatotalofover400
cular bone, nasal septum, and neural and basicranial regions,
FEAsimulationsbacktobackcanbecostlyintermsofcomputer
which were modeled isotropically as in models 3 and 4. The
time and memory. After the analyses were completed, postpro-
orthotropic material properties are given in Table 4. Similar to
cessinginANSYSprovidedthevaluesoftheoutputparametersfor
models3and4,materialpropertieswererandomizedwithaCV
each simulation and statistical analyses were conducted to
of20%andaGaussiandistributioninmodel5andallthematerial
calculatethemaximumshearstrainonthesurfaceofthecranium
properties were randomized with their experimentally deter-
ateachsiteforcomparisontoinvivodata(Rossetal.,2011).
mined standard deviations and a Gaussian distribution in
model6.
2.3. Comparingtheresultstoinvivodata
When the material properties were experimentally deter-
mined (Wang and Dechow, 2006), the nomenclature was stan-
Simulationresultswerecomparedtotheinvivodatagathered
dardized such that E 4E 4E , G 4G 4G , and
on macaquesand compiledin Straitet al.(2005) and Ross et al. 3 2 1 13 23 12
v 4v 4v . In general, the material properties tend to be
(2011).Graphswereconstructedwiththepublishedinvivodata 12 23 13
positivelycorrelatedtooneanother,soalargerE tendstoresult
fromnineteendifferentregionsofthemacaquecranium(Fig.2). 3
inalargerE andE .Thereexistsnoexactmathematicalfunction
The mean of experimental averages from in vivo shear strain is 2 1
thatdefinestherelationshipsbetweeneachmaterialproperty,so
depicted by the dotted black line and the dispersion (maximum
correlation coefficients relating the material properties to each
range defined by the meansþtwo standard deviations for each
other were calculated from empirical data (see Table 5). These
experiment)isdepictedbythegrayareas.Maximumshearstrains
correlationcoefficientswerespecifiedintheANSYSProbabilistic
werecalculatedfromthefiniteelementmodelateachsiteusing
Designmodule.
strain data extracted in terms of 34 site—specific local xyz
coordinate systems. Each local coordinate system had its origin
atthenodelocatedatthecenterofeachstraingagelocationand
itsxyplanelaidtangenttothesurfaceofthecranium.Maximum 2.4.4. Allmodels
shearstrainatthesitewascalculatedusingtheformula: Comparisons of results between models 1 and 3 and models
2 and 4 provide two assessments of the effect of increasing non-
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
g ¼ ðE (cid:2)E Þ2þg2 ð2Þ homogeneity in the skull that has on the various biomechanical
max x y xy
responses. Comparisons of results between models 3 and 5 and
where x, y, and z are coordinates of the site-specific coordinate models4and6providetwoassessmentsoftheeffectthatincreasing
systemandex,ey,andgxyarenormalstraininthex-direction,normal anisotropy in the skull has on various biomechanical responses.
straininthey-direction,andtheengineeringshearstrain1inthexy Finally,comparisonsbetweenmodels1and 2,3and 4,and5and
plane, respectively. This enabled the maximum shear strainin the 6 provide an assessment of the effect of assuming constant
plane of the surface to be calculated and any contributions to coefficientofvariationforallmaterialpropertyvalues,versususing
maximumshearstrainduetoout-of-planestrainstobediscarded. coefficientsofvariationsgleanedfromexperimentaldata.
This was necessary since strains measured by strain gages in the For models 1, 3, and 5, the coefficients of variation were
in vivo experiments do not account for any strain components calculatedatsamplinglocationswithmodest-to-highmagnitude
involvingthedirectionnormaltotheplaneofthestraingage. principal compressive, principal tensile, and von Mises stresses
and strains (Table 6 and Fig. 3). (‘‘Modest-to-high magnitude’’
2.4. Themodels meansstrainswithamagnitudegreaterthan250meandstresses
of greater than 3MPa.) The coefficients of variation depict the
2.4.1. Models1and2:simpleisotropicmaterialproperties magnitudesbywhichtheaveragesandstandarddeviationsofthe
In models 1 and 2 (see Table 1 for model descriptions) the stressesandstrainsatthesevariouslocationsacrossthecranium
macaque cranium was divided into three distinct materials – were affected relative to each other when a 20% variation was
trabecularbone,corticalbone,andteeth–consistingofisotropic appliedtothematerialproperties.
M.A.Berthaumeetal./JournalofTheoreticalBiology300(2012)242–253 247
Fig.2. Resultsfromtheprobabilisticanalysescomparedtotheinvivodata.Thegrayareadepictsthedispersionoftheinvivodata,andtheerrorbarsdepictthe95%
confidenceinterval(þ/(cid:2)2standarddeviations)oftheshearstrainatthevariouslocationsacrossthecranium(a)comparisonbetweenmodels1and2(3region,isotropic
materialproperties),(b)comparisonbetweenmodels3and4(multipleregion,isotropicmaterialproperties),and(c)comparisonbetweenmodels5and6(multipleregion,
orthotropic).
3. Results Inmodel1,theprincipalandvonMisesstrainsaresignificantly
more affected by variation in material properties than are the
3.1. Models1,3,and5 principalandvonMisesstresses,asisseenbythelargercoefficients
ofvariationforthestrains(seeTable6andFig.3).Asthematerial
Models 1, 3, and 5 provide a unique opportunity to see how propertiesvarywithacoefficientofvariationof0.2,modest-to-high
variabilityinstressesandstrainsareaffectedbynon-homogene- magnitude strains vary consistently with an average coefficient of
ityandanisotropywhenstandardvariabilityinmaterialproper- variation of 0.210,whereas the modest-to-highmagnitude stresses
tiesisapplied. varyconsistentlywithacoefficientofvariationofonly0.016.
248 M.A.Berthaumeetal./JournalofTheoreticalBiology300(2012)242–253
Table2 Asnon-homogeneityincreasesfrommodel1to3(seeFig.4),
Meanmaterialpropertyvaluesusedformodels1and2(threeregionisotropic). variation in the strains decreases from 0.210 to 0.180, but the
The standard deviations in the table were used to randomize the material variationinstrainsislessconsistent,asisdenotedbytheincrease
propertiesinmodel2;aCVof0.2wasusedtorandomizethemeanvaluesof
in the standard deviation of the coefficients of variation from
thematerialpropertiesformodel1.
0.007 to 0.047. Variation in the stresses, however, increases
Region Modulusofelasticity(MPa) Poisson’sratio greatlyfrom0.016to0.106,andlikethestrains,thevariationin
stressesislessconsistentinmodel3thaninmodel1.
Avg. Standarddeviation Avg. Standarddeviation
Finally,asanisotropyincreasesfrommodel3to5(seeFig.5),
Corticalbone 17,060 2760 0.28 0.10 variationinboththestressesandstrainsdecreases.
Teeth 70,000 14,000 0.30 0.06
Trabecularbone 637 127 0.28 0.056 3.2. Models2,4,and6
In model 2 the strains vary consistently with a coefficient of
variationof0.165whilethestressesvaryconsistentlyatalower
levelof0.025. Asnon-homogeneity increasesfrommodel2to 4
Table3 (seeFig.4),thelevelofvariationinstrainsincreasesgreatlyfrom
Meanmaterialpropertyvaluesusedformodels3and4(multipleregionisotropic).
0.165 to 0.560 and the strains vary much less consistently. The
The standard deviations in the table were used to randomize the material
propertiesinmodel4;aCVof0.2wasusedtorandomizethemeanvaluesof level of variation in the stresses also increases and the stresses
thematerialpropertiesformodel3. vary less consistently, but not to the level of variation in the
strains(seeTable6).
Region Modulusofelasticity(MPa) Poisson’sratio
Asanisotropyincreases(seeFig.5),theaveragecoefficientof
Avg. Standard Avg. Standard variation in the strains decreases to 0.369 and the strains vary
deviation deviation moreconsistentlythantheydidinmodel2.However,thelevels
ofvariationarestillveryinconsistent(meaningthecoefficientsof
Premaxilla 18,500 1700 0.18 0.07
variation have a large standard deviation), indicated by the
P3–M1alveolus 16,700 3000 0.24 0.13
averagecoefficientsofvariation,whichhaveastandarddeviation
M2–M3alveolus 20,600 4400 0.24 0.07
Anteriorpalate 15,300 4200 0.36 0.11 of 0.183. The variation in stresses again increases, this time to
Posteriorpalate 18,800 7900 0.26 0.17 0.150,butthestressesvarymoreconsistently.
Dorsalrostrum 19,900 3100 0.21 0.13
Lateralrostrum 18,100 5400 0.24 0.09
3.3. Comparingthemodelstoinvivodata
Rootofzygoma 17,900 3000 0.30 0.06
Anteriorzygomaticarch 20,800 8300 0.28 0.12
Posteriorzygomaticarch 12,500 900 0.27 0.20 Results from the probabilistic analyses are compared to the
Medialorbitalwall 14,600 8900 0.40 0.11 in vivo data in Fig. 4. When validating a craniofacial FE model,
Postorbitalbar 19,800 8200 0.22 0.14
validation is typically achieved by comparing strains on the
Frontaltorus 13,100 4200 0.24 0.18
model to in vivo or ex vivo strain data. Ideally, the in vivo or
Glabella 14,400 2900 0.14 0.10
Frontalsquama 14,900 7800 0.27 0.14 exvivodatawouldhavecomefromthesameindividualwhowas
Neuro-andbasicrania 17,300 5000 0.28 0.14 used to create the FE model, and the material properties and
Septum 13,194 3800 0.28 0.14 muscle forces applied to the FE model would have been taken
Teeth 70,000 20,300 0.30 0.15
fromthat same individualas well.Validationof themodels was
Trabbone 637 185 0.28 0.14
considered to be good at most locations on the cranium (where
Table4
Meanmaterialpropertyvaluesusedformodels5and6(multipleregionorthotropic).Thestandarddeviationsinthetablewereusedtorandomizethematerialproperties
inmodel6;aCVof0.2wasusedtorandomizethemeanvaluesofthematerialpropertiesformodel5.Theteeth,trabecularbone,septum,andneuro-andbasicranium
wereallgiventheisotropicmaterialpropertiesassignedinmodels3and4forbothmodels5and6.
Region Modulusofelasticity(MPa) Shearmodulus(MPa) Poisson’sratio
E1 E2 E3 G12 G23 G13 n12 n23 n13
Avg. Std. Avg. Std. Avg. Std. Avg. Std. Avg. Std. Avg. Std. Avg. Std. Avg. Std. Avg. Std.
dev. dev. dev. dev. dev. dev. dev. dev. dev.
Premaxilla 10,000 1400 13,900 3900 18,500 1700 4400 1100 5200 1000 7300 1900 0.29 0.16 0.18 0.07 0.15 0.15
P3M1alveolus 9900 3900 12,100 4000 16,700 3000 4300 1700 5800 2400 7400 2100 0.33 0.21 0.24 0.13 0.17 0.11
M2M3alveolus 12,600 3900 15,400 4000 20,600 4400 4900 1500 6400 2400 7900 2600 0.35 0.13 0.24 0.07 0.22 0.16
Anteriorpalate 7500 1400 8800 1500 15,300 4200 2600 400 2800 400 3600 900 0.41 0.06 0.36 0.11 0.26 0.11
Posteriorpalate 6400 1000 7500 1500 18,800 7900 2200 400 2500 400 3300 900 0.48 0.08 0.26 0.17 0.23 0.13
Dorsalrostrum 12,200 2100 14,000 2100 19,900 3100 5000 900 6900 1800 8900 2900 0.32 0.14 0.21 0.13 0.14 0.06
Lateralrostrum 11,500 5300 14,400 5400 18,100 5400 4700 2500 5300 2100 7300 3100 0.37 0.13 0.24 0.09 0.15 0.13
Rootofzyg.arch 8900 2900 10,900 4400 17,900 3000 3700 1700 5300 1200 8600 3300 0.53 0.15 0.3 0.06 0.18 0.12
Anteriorzyg.arch 8600 4700 12,400 5700 20,800 8300 4200 1900 4600 1500 8600 2800 0.39 0.19 0.28 0.12 0.22 0.17
Posteriorzyg. 8200 2700 10,000 2900 12,500 900 3100 1100 3800 1100 4900 1600 0.34 0.10 0.27 0.20 0.24 0.19
arch
Medorbitalwall 7100 4400 11,500 6500 14,600 8900 3600 2000 4200 2400 9000 5900 0.46 0.12 0.40 0.11 0.23 0.16
Postorbitalbar 11,300 3600 13,100 3300 19,800 8200 4400 1500 6400 2200 8000 2500 0.44 0.07 0.22 0.14 0.15 0.13
Frontaltorus 10,200 2900 11,200 3100 13,100 4200 4300 1400 5100 1900 6000 1800 0.32 0.20 0.24 0.18 0.19 0.08
Glabella 9200 2300 9700 2100 14,400 2900 3300 900 4800 1400 5100 1200 0.46 0.08 0.14 0.10 0.21 0.09
Frontalsquama 7900 3800 11,000 5000 14,900 7800 3400 1500 4300 2300 7100 3500 0.49 0.12 0.27 0.14 0.18 0.18
M.A.Berthaumeetal./JournalofTheoreticalBiology300(2012)242–253 249
Table5
Coefficientsofcorrelation(R)fororthotropicmaterialproperties,calculatedusingdatatakenfromWangandDechow(2006).E1/E2representsthecoefficientofcorrelation
betweenE1andE2,G12/G31thecoefficientofcorrelationbetweenG12andG31,etc.
Region E1/E2 E2/E3 E1/E3 G12/G31 G31/G23 G12/G23 m12/m13 m13/m23 m12/m13
Vault 0.900 0.910 0.830 0.829 0.579 0.529 0.273 0.423 0.440
Circumorbital 0.916 0.903 0.796 0.816 0.685 0.807 0.412 0.530 0.465
Zygomaticarch 0.904 0.881 0.801 0.747 0.545 0.687 0.642 0.512 0.279
Muzzle 0.933 0.925 0.946 0.681 0.319 0.312 0.600 0.383 0.708
Alveolus 0.564 0.724 0.579 0.798 0.564 0.402 0.177 0.283 0.306
Intraorbital 0.915 0.909 0.883 0.918 0.563 0.608 0.146 0.665 0.094
Palatine 0.704 0.585 0.159 0.889 0.439 0.623 0.842 0.208 0.328
Table6
ThecoefficientsofvariationfortheminimumandmaximumprincipalandvonMisesstressesandstrainswerecalculatedformodest-to-highmagnitudestresses(greater
than3MPa)andmicrostrains(greaterthan250me)ateachofthe34sitesacrossthecranium.Thecoefficientsofvariationwereaveragedandthestandarddeviationswere
calculatedtoconveythemagnitudetowhichthestressesandstrainsvariedatthevarioussitesacrossthecraniumandhowpredictablethislevelofvariationis.
Model Averagecoefficientsofvariation
Strains Stresses
InputvariableswithCV¼0.2 (1)3Region,isotropic 0.21070.007 0.01670.009
(3)Multipleregions,isotropic 0.16970.028 0.10470.043
(5)Multipleregions,orthotropic 0.14470.042 0.08170.025
Inputvariableswithstandarddeviationsbased (2)3Region,isotropic 0.16570.011 0.02570.017
onempiricaldata (4)Multipleregions,isotropic 0.56070.334 0.12470.043
(6)Multipleregions,orthotropic 0.36970.183 0.15070.038
the FEA results fell within the experimental range) but poor in materials properties are independent could be false, i.e. muscles
some locations (i.e. the left and right orbital roofs). As non- forces, morphology, and material properties could be positively
homogeneity and anisotropy increased from models 1 and 2 to correlated.
models 3 and 4 to models 5 and 6 (see Figs. 4 and 5), the FEA
results þ/(cid:2)2 standard deviations include or nearly include the
averagefromtheinvivodata.Whenthecraniumismodeledwith
the highest level of simplicity (model 1, where anisotropy and 4. Discussion
non-homogeneity are low and no empirical data is used to
calculatethestandarddeviations)theaverageinvivoshearstrain 4.1. Comparingmodels1,3,and5andmodels2,4,and6
iswithintherangeornearlywithintherangeoftheFEAresults
foronly10outofthe20locations.However,whenthemacaque Inallcaseswhenacoefficientofvariationof0.2isappliedto
craniumismodeledwiththehighestlevelofaccuracy(model6), the material properties (models 1, 3, and 5), there is less
the average in vivo shear strain is within or nearly within the variabilityintheshearstrainthanwhenacoefficientofvariation
rangeoftheFEAresultsfor18outofthe20locations. calculated from empirical data is applied. Therefore, applying a
Thisimpliesthatmodel6isamorerealisticmodelthanmodel coefficient of variation of 0.2 to the material properties under-
1,andthatmodelingwithahigherlevelofaccuracywillgivethe estimatesthevariabilityintheshearstrainsacrossthecraniumin
modelerabetterchanceofhavinginsilicoresultsthatwillmatch theFEmodel.Thisconclusionisfurthersupportedbytheresults
theinvivoresults.Thisalsoimpliesthatwhenmodelingsimplis- inTable 6 and Fig. 3, wherethe variability in theminimumand
tically,amodelthatappearstobeinvalidcouldactuallybevalid maximumprincipalandvonMisesstressesandstrainsisalways
whenuncertaintyinmaterialpropertiesisaccountedfor. lower when a coefficient of variation of 0.2 is applied to the
WhenempiricaldataareusedtocalculatetheCVsformaterial material properties compared to when a coefficient of variation
properties and the cranium is modeled simplistically (model 2), calculated from empirical data is applied. When the cranium is
thevariabilityintheshearstrainsisminimal.Alargenumberof modeled most simply, i.e. consisting of three isotropic homoge-
the average shear strains also fall within the dispersion of the neousmaterials(models1and2),thevariabilityinthemoderate-
in vivo data. As non-homogeneity increases (model 4), the to-highminimumandmaximumprincipalandvonMisesstresses
magnitude of the shear strains remains consistent but the andstrainsarelessthanwhenthecraniumismodeledwiththe
variability in the shear strains increases greatly. Finally, as highest level of accuracy (see models 5 and 6 in Table 6 and
anisotropyincreases(model6),themagnitudeoftheshearstrains Fig.3).Therefore,itshouldbeunderstoodthatwhenthecranium
changes at multiple locations across the cranium and the varia- ismodeledasconsistingofthreehomogeneousisotropicmateri-
bilityintheshearstrainsincreases. als,thelevelofvariabilityduetouncertaintyinmaterialproperty
When comparing the variability in the in vivo data to the valuesinthestressesandstrainsisminimized.Whenthecranium
variabilityintheFEAresults,thevariabilityintheFEAresultsputs ismodeledwithgreatermaterialcomplexitybyreducingthelevel
a lower bound on the total variability of in vivo data, assuming ofhomogeneity(model4),theaccuracyofthemodelremainsthe
variability in the in vivo data due to muscle forces, morphology, same but the level of variability in stresses and strains due to
and material properties are all independent. Considering that uncertainty in material property values increases (Table 6 and
therearesomelocationswherethevariabilityintheFEAresults Fig.3).
exceedsthevariabilityintheinvivodata,thiscouldimplythereis Models 2, 4, and 6 show that when the variation in material
an error with the data set, an error with the model, or that the properties is known, variability in both stresses and strains
assumption that variability in muscle forces, morphology, and increasesgreatlyasnon-homogeneityandanisotropyincrease.
250 M.A.Berthaumeetal./JournalofTheoreticalBiology300(2012)242–253
Fig.3. Coefficientsofvariationplottedagainsttheaveragemicrostrainsandstressesformodels1,3,and5.Theaveragecoefficientofvariationforthemodest-to-high
magnitudemicrostrainsandstresses(microstrainsgreaterthan250andstressesgreaterthan3MPa)canbefoundinTable6.Thebluediamondsrepresenttensilestrains/
stresses,theredsquaresrepresentthecompressivestrains/stresses,andthegreentrianglesrepresentvonMisesstrains/stresses.(Forinterpretationofthereferencesto
colorinthisfigurelegend,thereaderisreferredtothewebversionofthisarticle.)
4.2. Testingthehypotheses individuals, some of the variation in the data could be due to
variation in material properties and/or geometry of the cortical
The first hypothesis, that uncertainty in material property bonebetweenindividuals.
values will have little effect on high stresses and a large effect When the material properties were varied with standard
in high strains for homogeneous isotropic models, is fully sup- deviations calculated from empirical data, the level of variation
portedbymodels1and2,asseenbytheaveragecoefficientsof in modest-to-high magnitude stresses and strains differs greatly
variation reported in Table 5. For model 1 in which isotropic fromwhena 20%variationisappliedto thematerialproperties.
material properties were randomized with CV¼0.2, the average Models1,3,and5suggestthat,ifallmaterialpropertyvaluesofa
CV for all stress measures at all sampling sites is only 0.016, system have about the same level of randomness, variability in
whereas the average coefficient of variation for modest to high materialpropertyvaluescausesvariabilityinprincipalstressesto
strains measures 0.210. In model 3, which used CVs of isotropic increase and variability in principal strains to decrease as non-
materialpropertyvaluesbasedonempiricaldata,theaverageCV homogeneityandanisotropyincrease.However,thevariabilityin
forstresswas0.025andtheaverageCVforstrainwas0.16. materialpropertyvaluesisnotconstantinmacaquecrania.Since
The second hypothesis is supported by the probabilistic the trends in variability in the stresses and strains differ from
analysis results since, as non-homogeneity and anisotropy models1,3,and5tomodels2,4,and6(seeTable6),thisleadsus
increase,theCVformodest-to-highmagnitudestressesincreases to conclude that the assumption of constant material property
(seeTable6).Thethirdhypothesismaybesupportedthroughthe variability produces inaccurate trends in variability in stresses
results of this experiment. The error bars in Fig. 2 depict a 95% andstrainswithrespecttohowtheskullismodeledintermsof
confidenceintervalcalculatedfromtheprobabilisticanalyses.All homogeneityandanisotropy.
modelshavealarge95%confidenceinterval,whichsupportsthe Therefore, if a probabilistic analysis is to be carried out to
conclusionthatvariationinmaterialpropertiescausessignificant determine how variability in the material properties of bone
variationinmaximumshearstrainonthesurfaceofthecranium. affect stress and strain, a knowledge of the actual biological
Therefore, since the in vivo data were gathered from multiple variability in material property values is required. This nullifies
M.A.Berthaumeetal./JournalofTheoreticalBiology300(2012)242–253 251
Fig.4. Invivodataplottedagainstsimulationresultsdepictinghowshearstrainvariesasnon-homogeneityincreasesfrommodel1to3andfrommodel2to4.
thefourthhypothesis,andimpliesthatifsuchknowledgeofthe 2003). To date, there is no consensus as to what macroscopic
variability in material property values is not known, trends in mechanics quantities are correlated with these deformations. If
variability in stresses and strains cannot be trusted as being microstructuraldeformationsandboneremodelingarecorrelated
correct. to strain, then one would expect bone remodeling to be highly
sensitivetovariationinmaterialproperties,sincesuchvariability
causesalargelevelofvariationinmoderate-to-highmicrostrains.
5. Conclusion If cellular deformations are correlated to stress, then one would
expect bone remodeling to be fairly insensitive to variations in
Ideally, an FE model of an individual’s cranium would be material properties, since we found moderate-to-high stress to
constructed in which information pertaining to the material haverelativelylowvariationduetovariationinmaterialproperty
properties and muscle forces of that individual’s cranium were values.
previouslygathered,andinwhichinvivoorexvivodataonthat Of course, stress and strain are continuum mechanics con-
individualthatcouldbeusedtovalidatethemodel.However,in cepts, and the current investigation is focused on macro-level
practice this standard is difficult to achieve, and therefore an mechanics,whiletheactualmechanicalsignaturetowhichbone
understandingofhowmuchvariationinFEmodelresultsisdue responds is at the hierarchical level of the bone matrix that
tolackofaccurateindividual-specificdataisneeded.Ourresults directly affects cellular response (You et al., 2001; Cowin, 1999;
indicatethatlargevariationsinmodest-to-highstrainsandlower Cowinetal.,1995;Pavalkoetal.,2003).However,littleisknown
variations in modest-to-high stresses occur due to variation in aboutvariationinmicroscalematerialpropertiesinbonetissues
material property values. More work is needed to quantify how in three dimensions and how this variation may affect cellular
variation in morphology between individuals and variation in response.
loads (due to both intra and inter-individual variation) affect Finally, we note that our data suggest that it is incorrect to
moderate-to-highstressesandstrainsincraniofacialFEmodels. assume that variability in muscle forces, morphology, materials
Also, the high sensitivity of modest-to-high strain versus the properties, and other ‘‘cranial variables’’ are independent (the
lowsensitivityof modest-to-highstressesdue to randomnessin error bars for the in silico results in Fig. 2 cover a much larger
material property values across individuals may have implica- spread than the error bars for the in vivo results). The fact that
tions in terms of bone remodeling. It is generally accepted that positive correlations between these quantities would result in
boneremodelingoccursasaresultsofdeformationsthatoccurat less variability in stress and strain than if these quantities were
a microstructural level (Burra et al., 2010; Bonivtch et al., 2007; independent has evolutionary implications. Evolution could be
Youetal.,2001;Cowin,1999;Cowinetal.,1995;Pavalkoetal., producingarobustsystem—onewhosemechanicalperformance