Biomechanical Properties of Bone and Mucosa for Design and Application of Dental Implants

Dental implants’ success comprises their proper stability and adherence to different oral tissues (integration). The implant is exposed to different mechanical stresses from swallowing, mastication and parafunctions for a normal tooth, leading to the simultaneous mechanical movement and deformation of the whole structure. The knowledge of the mechanical properties of the bone and gingival tissues in normal and pathological conditions is very important for the successful conception of dental implants and for clinical practice to access and prevent potential failures and complications originating from incorrect mechanical factors’ combinations. The challenge is that many reported biomechanical properties of these tissues are substantially scattered. This study carries out a critical analysis of known data on mechanical properties of bone and oral soft tissues, suggests more convenient computation methods incorporating invariant parameters and non-linearity with tissues anisotropy, and applies a consistent use of these properties for in silico design and the application of dental implants. Results show the advantages of this approach in analysis and visualization of stress and strain components with potential translation to dental implantology.


Introduction
The placement of a dental implant system for the restoration of a missing tooth or teeth relies on the osseointegration of the implant (screw and abutment) with the maxillary or mandibular bone, whereas abutment part success depends on its adherence and integration with soft mucosal (gingival) tissue [1][2][3]. The dental implant is exposed to different mechanical stresses from swallowing mastication and parafunctions (grinding, clenching) via the application of the forces to the crown for a normal tooth, leading to simultaneous mechanical movement and deformation of the whole structure. Depending on the degree of the osseo-and mucosal integration achieved to the time point, this mechanical force might be or might not be fully transferred to the bone and mucosa.
Non-optimal displacements can have an undesirable effect on the integrity of biomaterialtissue interface, which serves as a barrier that prevents bacterial penetration [4][5][6], and a subsequent peri-implant disease or even implant loss. Therefore, the formation of an early and long-standing biomaterial tissue interface is of paramount importance. Natural tooth and the implant regulate and stimulate respective, although not identical, biological mediator to heal by promoting the tissue remodeling process [7,8].
Hence, the knowledge of mechanical properties of the bone and gingival tissue in normal and pathological conditions is very important for the successful conception of dental implants to access and prevent potential failures and complications originating from incorrect mechanical factors combinations.
Numerous studies were carried out and clinical data have been collected regarding the shape, design, surface and geometry of abutments and dental screws and their suitability for different tissues' quality and location [1][2][3][9][10][11][12][13]]. Significant differences were Table 1. Some reported data for elastic modulus of cortical and cancellous bones.
More complex correlations of specific bone stiffness tensor components vs. bone volume fraction, mean intercept length of the microstructure and applied strain rate can be found in literature [21], but such detailed analysis is usually local and not possible to be made in vivo. Here, a primary look is on such data on bone anisotropy, which could be deployed in clinical practice with reasonable tolerances [31].
Variations and differences between the data obtained with different test methods for bone and other tissues are well known. This makes it difficult to decide the values to be used for the design and placement of the dental implants. From a biomechanology point of view [20], the most relevant data are those obtained in conditions corresponding closely to the host physiological ranges. The values obtained in more exotic conditions (high frequency ultrasound, nanoindentation, etc.) are correct for those selected methods and environments, but their practical use is inferior to numbers measured at the situations serving clinical application purposes (e.g., at 1 Hz with 30-100 µm amplitude [22]). Carrying out such "right" tests, however, is not always easy. Even in simple setups like uniaxial compression or simple bending, orientation anisotropy of the bone specimen might be difficult to account properly.

Reported Oral Soft Tissues Properties
The experimental testing of viscoelastic properties of dog, porcine and human mucosa and other oral section (palate, gingiva, periodontal ligaments) in vivo and in vitro was carried out or summarized in many studies [17,19,26,[34][35][36][37][38][39][40]. As expected, there is a significant scatter in the data (Table 2)-more than for bones-depending on tissue location, conditions, and the deployed testing method. Many studies do not provide a well-defined rationale as to why these specific parameters have been chosen for the tests and how other parameters' selection might have affected the data obtained. Furthermore, as stated in [21], the testing of isolated soft tissues presents substantial challenges, due to preparation variations, grips and clamping effect (especially in tensile and shear modes), the inaccurate measurements of geometry and cross-sectional areas (for stresses estimation), preconditioning methods, deformation rates, etc.
In rheology (such as oscillatory shear), even more errors and artefacts are often seen in the literature data. These may comprise, i.a., a priori set linear viscoelastic (LVE) models (not suitable for biological tissues), wrong boundary conditions in real tests vs. models; ignorance of present momentum diffusion or viscoelastic waves; disturbances from secondary flows, especially at low Reynolds numbers (considered as "laminar flow"); non-controlled elastic instabilities at high Weissenberg numbers; waves propagation in dynamic tests [9,20,27]. This makes the selection and comparison of consistent data values for soft tissues much more difficult than for a bone.
It is known that soft tissues including oral ones (mucosa, gingiva, palate, PDL) exhibit complex nonlinear time-dependent behavior [17,39,41,42]. However, most of the models in literature assume material and geometric linearity with the homogeneity and isotropy of the mucosa. Due to evident limitations of such approach, more complex models such as hyperelastic ones were deployed in description of these soft tissues (a hyperelastic material modelling formulates a potential energy function (strain energy potential) per unit of reference volume [17] as a function of the strain at a typical point in the material). This function can be dependent either on strain tensors of a nonlinear deformation field, on the invariants of these strain tensors, or directly on the principal stretches [17].
However, hyperelastic material model application is limited, especially in its demand for a substantial number of parameters or fitting functions if nonlinearity and anisotropy have to be included in the calculations [22,23,43]. It is difficult to estimate to what extent these fitting values could be translated to other subjects with different heterogeneous anatomical microstructures or generalized for practical use under complex physiological responses.

Other Parameters of Clinical Importance
Among the parameters which might be treated as having significant clinical relevance, osseointegration quality improvement and the minimization of the pressure-pain threshold (PPT), interstitial fluid (hydrostatic) pressure (IFP) and pressure leading to residual ridge resorption (destruction of the supporting bony tissues) can be mentioned.
Osseointegration is a complex outcome (or an endpoint) which can be evaluated only after some time (from months to years), combining objective measures from an examination, as well as patient satisfaction with physiological functionality and aesthetic effect ("quality of life"). Other outcomes can be, in principle, calculated or estimated using experimental and modelling (in silico) data. The biomechanical models of "hard" (bone) and "soft" (mucosa) tissues aim to interpret, analyze and (where possible) predict various aspects of these tissues' responses to dental prostheses [17].
It is also known that the success of a dental implant depends on a variety of biomechanical factors, including the design and position of the implant, implant-abutment connection, cantilever length, surface roughness, bone quality and type, depth of insertion, arch configuration, the nature of bone-implant interface, and occlusal conditions [1,2,9,12,14,29]. Whereas such biomechanical factors have been simulated earlier, the results are often contradictory, due to differences in model parameters, construction, materials models selection and meshing [29].
In regard to mechanical stimuli, many reported studies are expressing acting stress in an equivalent format such as von Mises stress [17,29,31,33,44]. However, the use of von Mises stress has substantial limitations: it has been most widely applied for engineering problems as a yield criterion, stating that the yielding of a material occurs once the second deviatoric stress invariant reaches a critical (yield stress) value [17,21,31,45].
Yielding (nonelastic irreversible deformation) is clearly defined for metals, plastic and ceramic materials, but it is questionable for live dental tissues. Furthermore, since von Mises stress is a module (absolute value) of the deviatoric stress, it is always positive: one cannot see from the von Mises data whether the area or point in question is under tension or compression.
Because the tissue's behavior in tension and compression is different and often nonlinear, the use of symmetric stress characterization is a clear oversimplification. Indeed, authors [29] noted that besides single values alike von Mises stress, the magnitudes of the stress concentrations, stress distribution, and displacement components of specific points should be analyzed as they could provide a valuable information about the deformation of the model and assist in the interpretation of the results.
For example, micromovements along the bone implant interface must stay within proper tolerances, as micromotions beyond those tolerances result in connective tissue encapsulation or adherence failure [1,2,9,11]. It is thus preferable to see whether there are tensile, compressive and shearing areas to identify possible risks. This could be done by visualizing selected stress and strain tensor components.
The average data for PPT identified by different studies are usually within a 0.1~0.4 MPa range [17]. Most of these data were acquired for the cases of dentures on mucosa, and it is not clear which pressure level inside the bone tissue would be considered as one causing a noticeable pain. It is known than cortical bone can withstand much higher pressures, since bone permeability is low, and respectively its physiological perfusion requires a significant pressure gradient [16,21,30].
For hydrostatic pressure (IFP) values between 0.3 to 2 kPa in gingival mucosa have been reported according to data [17]. Authors [46] have reported a "safe" level of internal pressures for soft tissues <5~10 kPa (preferably <4 kPa) and for a short-term loading <13~20 kPa: once it exceeds the vascular pressure differential, blood flow will be reduced and may temporarily cease, potentially leading to local anoxia and localized ischaemia [16,[46][47][48]. Both PPT and IFP can be, in principle, calculated with in silico FEA post-processing and compared with the criteria above to understand what potential risk might be associated in different locations.
For residual ridge resorption, about 20 kPa of intermittent pressure and 7 kPa for continuous pressure at mucosa were reported to cause alveolar bone ridge resorption [17]. These data are compatible with the one shown above for IFP (<20 and <5~10 kPa) respectively [46], so that they might be recommended as watch criteria in results analysis.
The summary of these parameters, which would be qualified as those for "acceptable clinical risks", could be highlighted as shown in Table 3. These values could be recommended to monitor in simulations and the optimization of the personalized implants design. Table 3. Parameters of clinical relevance related to biomechanics (in kPa).

Bone Stiffness Approximation
Bone and other biological tissues are non-linear viscoelastic materials, so for timedependent, direction-and strain-rate dependent behavior, rather complex constitutive equations have to be used. Usually, the tissue behavior of an anisotropic viscoelastic material is described based on the application of the Boltzmann superposition integral [21,31], with (time-dependent) second rank stress (σ ij ) and strain (ε ij ) tensors and four-rank stiffness tensor T ijkl . This links stress (σ) and strain (ε) of the material as: where T is the stiffness tensor and S is the compliance tensor (the inverse of T). For the reduced notation of T leading to the 2-rank stiffness tensor C, there are 36 remaining independent elements (2) for the lowest symmetry case and 12 nonzero independent elements (C 11 . . . C 33 , C 44 , C 55 and C 66 ) for an orthotropic solid [21,31,49].
Resulting stress-strain equations for anisotropic materials are easily becoming rather cumbersome, requiring a special computation to be used in the practice. Nevertheless, for the more precise planning of the implant and its expected mechanical stability it is important to know the anisotropy of the combined cortical and trabecular (cancellous) bone as a whole, since the screw displacements and deformations have to be compliant to both bone tissues at the same time. As the abutment is tightly fixed to the screw, the movement of the screw and the abutment are highly correlated: they are moving together as a single assembly object.
Abutment purpose is to have sufficiently tight adhesion to the local gingival tissue laying on the bone, so both gingiva and the bone would face the same resulting displacement of the fixture, but this displacement will be felt differently by bone and gingiva, due to intrinsic differences in their properties. Hence, the same displacement will result in different strains in gingiva, cancellous and cortical bones. Due to properties' anisotropy and non-linearity, it is improper to describe all these players with single values of "elastic modulus".
Most of the finite element analysis (FEA) software uses an assumption that the bone is linearly elastic isotropic material, which cannot catch the differences in local strains and stresses faced by the tissues [33]. For anisotropic materials instead of "elastic modulus", the whole stiffness matrix (2) needs to be considered. The latter can be simplified for higher symmetry materials like an orthotropic one, when many off-diagonal components (3a) vanish [45,50]. There compliance matrix S is shown in (3a) is in the reduced notation for orthotropic symmetry, with E ii being directional elastic moduli (11 = X, 22 = Y, 33 = Z directions as in Table 1) and G ij are respective shear moduli [45]. The Poisson ratios (ν) in general are not symmetric ν ij = ν ji (3b).
In Table 1, reported values of the stiffness matrix components are those usually deployed in bone behavior modeling. The coordinate indexing adopted here is as commonly used, i.e., it takes X-(11), Y- (22) and Z-(33) orthogonal directions [45,49] and cross-components as XY (12), YX (21), YZ (23), etc., with a row normalization (3a, 3b). The allocation of global coordinates is usually as Z for infero-superior direction, Y for bucco-lingual direction, and X for mesio-distal direction [33]. Local coordinates, which describe the local instant direction, depend on the point of view (=they can be rotated), keeping the orthogonality and alignment with the tissue structure.
The rotation of the coordinates for expression of the 4-rank tensor components T ijkl can be performed with two consequent rotations [50], as shown in Figure 1. There original position of Z-coordinate is aligned with the axis of the implant in the superior direction, Y-coordinate in the buccal direction, and X-coordinate in distal direction. The first rotation (shown in Figure 1 as blue axis) is done around Z-axes by yaw angle Ψ, resulting in new X' and Y' axes. The second rotation (shown in Figure 1 as red axis) is done along this new X' axis by pitch angle Θ. This leads to a new set of local coordinate axes X'-Y"-Z", which can be directed to any point of view just by these two rotation operations (note that the order of rotation is important [50], Figure 1).
î ï The rotation of the coordinates for expression of the 4-rank tensor components Tijkl can be performed with two consequent rotations [50], as shown in Figure 1. There original position of Z-coordinate is aligned with the axis of the implant in the superior direction, Y-coordinate in the buccal direction, and X-coordinate in distal direction. The first rotation (shown in Figure 1 as blue axis) is done around Z-axes by yaw angle Ψ, resulting in new X' and Y' axes. The second rotation (shown in Figure 1 as red axis) is done along this new X' axis by pitch angle Θ. This leads to a new set of local coordinate axes X'-Y"-Z", which can be directed to any point of view just by these two rotation operations (note that the order of rotation is important [50], Figure 1).

Figure 1.
Schematic of the implant in the mandible and respective coordinate axes with rotations as shown in the text. The first rotation ("blue" with the yaw angle Ψ) around Z-axes is followed by the second ("red" with the pitch angle Θ) around new X'-axes counter-clockwise.
Finally, to calculate new rotated stiffness components, the unrotated ones are to be multiplied with the rotation matrix Q and summed over as in Equation (4). The rotated stiffness tensor and its new inverse (compliance tensor) can now be used for the calculation of anisotropic stresses and strains in the implant-bone-gingiva system. Figure 1. Schematic of the implant in the mandible and respective coordinate axes with rotations as shown in the text. The first rotation ("blue" with the yaw angle Ψ) around Z-axes is followed by the second ("red" with the pitch angle Θ) around new X'-axes counter-clockwise.
Finally, to calculate new rotated stiffness components, the unrotated ones are to be multiplied with the rotation matrix Q (4b) and summed over as in (4a). The rotated stiffness tensor and its new inverse (compliance tensor) can now be used for the calculation of anisotropic stresses and strains in the implant-bone-gingiva system.
Before using known data which consist of separate elastic and shear moduli (Tables 1 and 2), they must be converted first into the stiffness tensor C components. Here these data of elastic moduli (E i , G ij ) were converted with Equation (4) into the reduced stiffness matrix components C ij and respective Poisson ratios ν ij using "RotaStiff" software developed by the author [50]. This conversion can also be done even with a spreadsheet, but such computations are more cumbersome.
The rotation (4b) does not lead to a symmetric tensor (C ij = C ji ) also due to differences in Poisson's ratios (for example, for data [33] one has C 21 = 8.169 and C 12 = 8.162, and C 13 = 9.678 vs. C 31 = 9.725 GPa, Table 4). For practical applications, specific tissue properties of the patient are desired instead of averaged literature values, despite the latter still being able to be used as a good approximation if their scatter in values is acceptable. Table 4. Reported data for unrotated stiffness tensor components (GPa) and Poisson ratios for cortical bones [29,31,33]. Data in italics in shaded cells are calculated in this work. Most of the differences in "bone quality" measured using X-ray or CT examination are related to the bone density [33]. It was shown [33,51,52] that for cancellous bone, its apparent density ρ a in g/cm 3 (without bone marrow, which does not contribute to load bearing capacity) can be rather well calculated from CT contrast values in HU (Hounsfield units) as: ρ a = 1.011 · 10 −3 × HU Together with the experimental values reported in [53] and using the conversion of single elastic properties into the stiffness matrix, authors have obtained the following stiffness matrix (GPa) for cancellous mandibular bone as functions of ρ a : This allows a direct computation of the stiffness matrix from the HU values (5). For cortical bone, a similar approach can be used, and with data [44,54], authors have obtained the following stiffness matrix values for mandibular cortical bone (GPa) vs. bone true (not apparent) density ρ (g/cm 3 ): For both cancellous (6) and cortical (7) bones stiffness, components C ij are not symmetrical, and they depend on bone densities in different ways.

Soft Tissues Stiffness Approximation
As follows from Table 2, there is a great scatter in the soft tissue data depending on the origin, preparation, testing method and the interpretation of the data. There is also a natural difference in tension, shear and compression method outcomes and it is not very clear which of these data are actually relevant for in vivo situation. When the implant is loaded in the physiological range, the displacements are usually below 50-100 µm along the direction of applied force [22]-micromovements over 150 µm are to be avoided [55]. During these micromotions, gingival tissue around the abutment undergoes some mixed deformation mode with local tension, shear and compression at the same time. Based on reported data ( Table 2) it is not however possible to artificially derive some averaged properties by combination of tensile and compressive results, which could be directly applicable to this case. Furthermore, there are no reliable data showing the anisotropy of such soft tissues, compared to the data available for the bone (Tables 1 and 3).
In this analysis, authors have selected the publications that have numerical and graphical data of the stress-strain-time relation for at least one reported experiment with sufficient details (test method, frequency, time, preconditioning, deviations, etc.). These data were digitized with GraphClick software (Arizona Software Ltd., Neuchâtel, Switzerland) and evaluated with BEST software (Seqvera Ltd., Helsinki, Finland) to extract invariant values [24,27]. Some sources in Table 2 were rejected, where only a value or a range was shown without the explanation of the method of analysis or calculation, or which otherwise were lacking essential information. Results where 'elastic moduli' or similar property was evidently obtained with the differentiation of stress by strain and extrapolation to zero strain or strain rate, were not analyzed due to larger errors and lesser relevance for clinical applications.
The invariant values [24] were obtained with BEST method [27] and can be used for the estimation and prediction of biomaterials and tissue properties without the use of a material model, as was demonstrated in [9,10,20,56,57]. There, the dependence of the material strain ε from the applied stress σ under pseudodifferential time-convolution (accounting for loading history in this case): where α is a material memory value, E 0 is the averaged time-invariant (i.e., not timedependent) intrinsic modulus, τ 0 is invariant characteristic time, and Γ() is the gammafunction. The product <E 0 (τ 0 ) α > is the time-averaged viscostiffness (a pseudo-property) of a material [20,27,46,56,57]. Characteristic invariant time can be related to the material Deborah (De) number: the specimen reaches steady deformation behavior when the observation (measurement) time is larger than the invariant time.
The test data analyzed in this work were divided into two classes: pseudo-static (creep) and dynamic. For creep under constant applied stress, the method [27,56,57] predicts a log-time dependence of the strain with E 0 c as the averaged invariant creep modulus, τ 0 cinvariant characteristic creep time, α c -material creep memory parameter. For dynamic loading case, the method predicts a sub-dimensional dependence of the frequency for the strain amplitude [27,46,57]. When the input stress is harmonic, the strain depends on observation time t, frequency ω and constant stress amplitude σ dyn , with E 0c ω as the averaged invariant dynamic modulus (i.e., not dependent on time and frequency), τ 0 ω -invariant characteristic dynamic time, α ω -dynamic material memory parameter. This Equation (8) comprises the time-convolution of the specimen loading history without a general need of complex algebra, assumptions of a material model (Maxwell, Burger, standard linear solid, Prony series, Mooney-Rivlin, Ogden, hyperelastic, neo-Hookean, etc.) or use of local differentiation [20,23,58,59].
The values of E 0 , τ 0 and α are time-invariant, in a sense that they do not depend on the time of experiment. Equation (8) can be numerically, explicitly computed without the need for assumptions of linearity of E 0 , τ 0 and α, whether they are for a creep (static) or dynamic case. According to [20,27], this allows one to overcome common restrictions for the linearity of tissue properties in many models [19], namely a scaling property (homogeneity) and a superposition property (additivity) (addition and multiplication are two linear operations which are applicable to complex numbers-in general, to the commutative rings [58,60]). This might not generally hold for (linear) Fourier transformation, commonly used in linear viscoelasticity, as most tissues and biomaterials properties functions are frequency-or strainrate-dependent [61][62][63]. For example, if the loss tangent, commonly expressed as ratio of imaginary to real moduli in viscoelasticity, is not constant, a rigorous application of complex math transformation cannot be made correctly, leading to wrong predictions, artefacts or just improper conclusions [62]. Authors [63] have underlined that the presentation of oscillatory stress or strain response in a real Fourier series is not at all equivalent to complex Fourier transformation, especially when only the first harmonic real and imaginary moduli are calculated as this itself imposes linearity transformation requirements (in the real world, there are no complex physical quantities: it is not possible to have a mass like 3 + 2i kg or time like 10 − 3i seconds-so no complex moduli or viscosity could exist in reality, being just mathematical abstractions).
For the time-convolution Equation (6), such linearity conditions assumption is not required due to idempotent processing [60]: the values of these memory parameters are always positive (because of causality principle: no response is generated before the stimulus has been applied). When they are in the range 0 < α < 1, they represent the fading memory. The smaller the value, the greater is the effect of short-time memory (= immediate reaction to the stimulus). Zero value means that the material has only a short-time memory, i.e., ideally elastic behavior, whereas α = 1 means ideally viscous behavior [27,46,56,57,59]. For the cases 1 < α < 2, the material only "feels" a long-term memory and the immediate elastic reaction on the stimulus is damped: the material mechanical behavior could turn into inertial damped oscillations [27,57].
It is noteworthy to emphasize that invariant values of E 0 , τ 0 and α describe the whole system behavior (tissue components, its fluids, possibly with blood content, etc.) interacting with the applied mechanical stimuli, and not the single materials proprieties values. This is different from the traditional way of mechanistic description of every tissue component (like collagen fibrils, extracellular matrix, etc.) separately, followed by assembling these descriptions into some set of approximate equations.

Statistical Method and Data Visualization
The data assessment procedure for soft tissues after digitizing literature data was realized in this work as follows. First, the presence of leverage points was detected with hat matrix diagonal components, which were not falling under Stephen's rule (these points were removed). An influencer (outliers) points analysis was made by calculating Cook's distances, and those data points exceeding unity value (if any) were also removed. The consistency of obtained coefficients with quantum regression method (Seqvera Ltd., Helsinki, Finland) was independently checked by the application of the Theil-Shen estimator and the goodness of fit significance-by the Nelson-modified Anderson-Darling test [64] (data not shown: the original algorithm in [64] assumes all positive input numbers and fails if a negative normalized value appears in the data. This feature has been improved in the present analysis). The heteroscedasticity of residuals was estimated with RUNS test and possible residuals' autocorrelation-with Durbin-Watson parameter.
All the data passed these criteria were used in Equation (8) to obtain biomechanical invariant values [24], which were considered to be BLUE (BEST linear unbiased estimators). This calculation method does not impose linearity, scaling and superposition requirements for soft tissues [19,36,37]. A visualization of moduli rotation in 3D by Equations (4) was made with MathCAD 14 software (PTC Inc., Needham, MA, USA). This procedure comprises the calculation of the bone stiffness matrices as Equations (1)-(3), stepwise "blue" rotation of the tensors by 0~90 • anticlockwise around Z-axis, followed by "red" rotation by 0~90 • around X' axes ( Figure 1), recalculation of the rotated stiffness for each rotation step with Equations (4), inversion of the stiffness tensor Equation (1) to obtain the compliance matrix, the calculation of respective Poisson ratio components (3b) and the 3D plotting of the calculated values.

Clinical Case Simulation
Based on the estimated tissues properties, a simulation of the clinical case has been made with the following conditions. We assumed that a patient CT-scan has been performed, pre-assessed by the doctor and the primary tissue data (geometry, density, size, location) have been extracted. A mandibular implant is placed at the molar position into normal quality bone (along unrotated Z-axis, Figure 1). Cortical bone thickness was taken as 1 mm and its density as 1.8 g/cm 3 ; medullar (cancellous) bone thickness as 6.5 mm with its density 600 HU (the apparent density of 0.607 g/cm 3 as by Equation (5)). Gingival thickness is taken as 2 mm average and its properties are approximated with BEST method, as shown in Equation (8). The implant micromotion amplitude in Z-direction was set to 50 µm at 1 Hz, imposed on the axially symmetric surface (Figure 2) leading to traction along the XZ-and YZ-planes ( Figure 1). Due to the symmetry, the same traction displacement will also act in rotated planes X'Z and Y'Z ( Figure 1), but not in X'Y", X'Z" or Y"Z" planes, which require a rotation conversion by Equations (4)  Equation (8). The thread and shape of the screw are not included in the calculations for this demonstration (as they are implant-type specific). The base XY plane is assigned with roller conditions (no vertical Z-displacements) and symmetry planes imposed at XY and YZ planes-in this way, ¼ of the whole structure is simulated. All other surfaces are set to free displacements ( Figure 2).

Stiffness Estimation and Visualization for an Arbitrary Direction
The estimation of the stiffness (2) and compliance (3) matrices was calculated for fixed bone densities, as shown in Section 3.2 with a stepwise (40 steps) rotation by 90°, first around the Z-axis and then around the X-axis (Figure 1). For every step, the tensor components (2) were recalculated and then converted into components of elastic moduli (here viscous contribution was introduced with Equation (8) for every matrix component). Examples of moduli for XX (E11) and shear XY (G12) are shown in Figure 3 in phase space. The coordinate system reflects the unrotated (initial) state, and the point of view-respective rotation.
For example, for Figure 3a (E11), the view vector is directed from the coordinates origin towards the reader, corresponding to the rotations by yaw Ψ~30° and pitch Θ~15° (Figure 1). The resulting stiffness matrix C and respective Poisson ratios matrix for this case for cortical bone is shown matrixes (9), where C values are in GPa and Poisson ratios matrix for cortical bone, as seen from the point of view of Figure 3a (orthogonally to the picture: from the reader towards coordinates origin). It is seen how stiffness anisotropy affects moduli depending on which direction it might be seen. All components of the C matrix (9) are now non-zero as they are off principal axes. Despite being the dynamic All tissues (cancellous bone, cortical bone, gingiva) were assumed to face the same micromotion as a boundary condition at the interface (= there is no assumption on the adhesion quality of the screw and the abutment, the traction is supposed just to be transferred to the contact surface). All materials (tissues) were set to "linear elastic material" in the COMSOL software, but with the customized properties described for every matrix element as by Equations (5)-(8) introducing the non-linearity and history-dependent viscoelasticity. As initial data for bone in Table 1 did not have viscoelastic components, a loss factor of 0.03 was adopted for 1 Hz based on general knowledge [9,16,21]. For gingival tissue, non-linear viscoelastic contribution was embedded in the stiffness matrix with Equation (8). The thread and shape of the screw are not included in the calculations for this demonstration (as they are implant-type specific). The base XY plane is assigned with roller conditions (no vertical Z-displacements) and symmetry planes imposed at XY and YZ planes-in this way, 1 4 of the whole structure is simulated. All other surfaces are set to free displacements ( Figure 2).

Stiffness Estimation and Visualization for an Arbitrary Direction
The estimation of the stiffness (2) and compliance (3a) matrices was calculated for fixed bone densities, as shown in Section 3.2 with a stepwise (40 steps) rotation by 90 • , first around the Z-axis and then around the X-axis (Figure 1). For every step, the tensor components (2) were recalculated and then converted into components of elastic moduli (here viscous contribution was introduced with Equation (8) for every matrix component). Examples of moduli for XX (E 11 ) and shear XY (G 12 ) are shown in Figure 3 in phase space. The coordinate system reflects the unrotated (initial) state, and the point of viewrespective rotation.
With Equation (8), the literature data for gingival tissues were processed and the gingi- For example, for Figure 3a (E 11 ), the view vector is directed from the coordinates origin towards the reader, corresponding to the rotations by yaw Ψ~30 • and pitch Θ~15 • (Figure 1). The resulting stiffness matrix C and respective Poisson ratios matrix for this case for cortical bone is shown matrixes (9a, 9b), where C values are in GPa and Poisson ratios matrix for cortical bone, as seen from the point of view of Figure 3a (orthogonally to the picture: from the reader towards coordinates origin). It is seen how stiffness anisotropy affects moduli depending on which direction it might be seen. All components of the C matrix (9a) are now non-zero as they are off principal axes. Despite being the dynamic simulation (1 Hz), the stiffness matrix components and Poisson ratios (9b) are all real, as due to method by Equation (8), no complex math via Fourier transform was needed.
The literature data for gingival tissues were processed with Equation (8) and the gingival/mucosal invariant data have been extracted (Table 5). Analyzed creep values were consistent in the invariant time scale (13~17 s) but differ in creep modulus, possibly due to about 200 times difference in applied stress in those experiments. In dynamic experiments, the data showed the opposite trend: the invariant moduli values are compatible, but invariant times have much higher difference, Table 5. This can be explained by dissimilar experimental conditions and that in [46] gingiva was not separated but kept attached to the underlying bone and submitted to a wide strains range (up to 0.4). These last data were used in in silico simulation for the gingival part of the model, as they are the most coherent for the clinical case considered.

In Silico Simulation
The total displacement amplitude of the tissues (µm) is shown in Figure 4. As expected, there are non-uniform, anisotropic displacements in all three tissue types, being the highest near to the implant traction surface-contact area (Figure 2, right). The maximal positive pressure appears near this interface and is the highest in the cortical bone in the X-direction (mesio-distal), Figure 5. This zone, however, is rather small (~100 µm depth) and changes into negative pressure zone deeper in the tissues. Neutral zone stays in gingival tissue and in the deeper cancellous bone, which might be considered positive from the point of view of having lower PPT and IFP values (Table 3). A detailed analysis of these threshold zones is possible, but it likely makes more sense for an exact implant shape and thread geometry. [46] abutment load at 1 Hz and 37 °C τ 0 ω = 23 s support; up to 50 μm displacement

In Silico Simulation
The total displacement amplitude of the tissues (μm) is shown in Figure 4. As expected, there are non-uniform, anisotropic displacements in all three tissue types, being the highest near to the implant traction surface-contact area (Figure 2). The maximal positive pressure appears near this interface and is the highest in the cortical bone in the X-direction (mesio-distal), Figure 5. This zone, however, is rather small (~100 μm depth) and changes into negative pressure zone deeper in the tissues. Neutral zone stays in gingival tissue and in the deeper cancellous bone, which might be considered positive from the point of view of having lower PPT and IFP values (Table 3). A detailed analysis of these threshold zones is possible, but it likely makes more sense for an exact implant shape and thread geometry.  Bone is the tissue which bears most of the load, as seen from the strain energy density by 4-6 orders of magnitude ( Figure 6). There is also a noticeable gradient of strain energy density (by ~1000 times) inside the gingival tissue. The values of the strain energy density magnitude in gingiva are, however, small (as there are very little deformation, Figure 4), so they are unlikely to cause any intermittent effects. On the other hand, long-term effects might be noticeable if the strain redistribution will proceed with time.
Simulation can also show different specific stress and strain components and other variables, as well as their time variation within the loading cycle (~1 Hz) or in a longer scale (creep) for detailed analysis. It might be interesting for clinical cases to identify possible critical areas and risk points, especially for personalized shapes of the implant-for example, a small, 0.2~0.5 mm change in thread or diameter in some areas might have a significant impact on clinical success. Bone is the tissue which bears most of the load, as seen from the strain energy density by 4-6 orders of magnitude ( Figure 6). There is also a noticeable gradient of strain energy density (by~1000 times) inside the gingival tissue. The values of the strain energy density magnitude in gingiva are, however, small (as there are very little deformation, Figure 4), so they are unlikely to cause any intermittent effects. On the other hand, long-term effects might be noticeable if the strain redistribution will proceed with time. might be noticeable if the strain redistribution will proceed with time.
Simulation can also show different specific stress and strain components and other variables, as well as their time variation within the loading cycle (~1 Hz) or in a longer scale (creep) for detailed analysis. It might be interesting for clinical cases to identify possible critical areas and risk points, especially for personalized shapes of the implant-for example, a small, 0.2~0.5 mm change in thread or diameter in some areas might have a significant impact on clinical success.  Simulation can also show different specific stress and strain components and other variables, as well as their time variation within the loading cycle (~1 Hz) or in a longer scale (creep) for detailed analysis. It might be interesting for clinical cases to identify possible critical areas and risk points, especially for personalized shapes of the implantfor example, a small, 0.2~0.5 mm change in thread or diameter in some areas might have a significant impact on clinical success.

Discussion
In this work, a simulation approach has been presented comprising a combined application of non-linear anisotropic tissue descriptions with semi-analytical and numerical computation. Here, the last one formally deploys an embedded "linear elastic material" description (a standard FEA feature), having however non-linear, invariant components of the stiffness tensor. Such methods are usually referred to as "reduced order models", capable of a very quick operation for adequate description of complex systems involving multi-phase, multi-materials combination in real-time dynamics [65,66]. This has a clear advantage over sophisticated non-linear methods-for example, in [67] it was assessed that for an incompressible second-order reduced polynomial model with a second-order quasi-linear viscoelastic extension in Prony series (with four unknowns), the analysis of the test data took about 24-40 h of computational time on a standard desktop computer. For ten or fifteen unknowns and a higher mesh density (more degrees of freedom), such an approach is unfeasible even with sufficient computational power.
The implant design and selection are highly personalized and present an important part of the whole implantation procedure. Therefore, for practical reasons, it is desirable that collected data and patient conditions could be used for such selection in a fast and safe way. It is less likely that a doctor or dental lab would engage themselves in sophisticated computations, taking time and generating a substantial amount of data needed to be analyzed, filtered and converted into decision making.
Here, numerical data compression has been achieved with the time convolution method and complemented by algorithm allowing visualization of the tensor components in 3D-rotated coordinates. Whereas the stiffness component visualization is not obligatory for the implant design and selection, it helps us to see the directions where the tissues stiffness is e.g., the weakest (zones of increased risks). Such visualization assists the user in the understanding of the tissues' anisotropy and its variation vs. view angle or direction of load application in 3D. The latter might be considered, for example, as a direction in which the force is applied to give an immediate response how this force would be transformed into respective strain components and pressure (important for PPT, IFP and residual ridge resorption; Table 3).
The most significant findings of the present work are: • Transitions from the single "elastic modulus" values towards true stiffness matrix components, which are real material functions of density, strain, time, frequency, etc.
give a more precise and comprehensive description of the tissues with patient-specific data (such as bone density in loco). • Visualized rotation and generated 3D plots can be of an assistance of identification of the weaker and stronger directions in the X-Y-Z coordinates that may help the digital planning of surgery (for example, where the lowest stiffness part might coincide with thinner or lighter bone in that direction).

•
For soft tissues (gingiva, mucosa), invariant values (modulus E, characteristic time τ and material memory α) have shown being extracted from experimental data without the assumption of a material model. The same can be made for bone, if viscoelastic effects need to be considered without the need of complex numbers processing.

•
In the case shown, the data used were obtained closer to the physiological range, allowed to mimic the behavior without needs to break down the values into specific tissue components. These invariant values comprise time-convoluted data and are better predictors for materials performance comparison than traditional stiffness (stress/strain ratios). • Non-linear viscostiffness data can be placed into software packages for a detailed simulation, and there, the directional rotations are performed automatically. Here, the detailed simulation of the implant screw and abutment is not always necessary-the imprint surface transferring the load is of the greatest interest as a boundary condition. By changing the patient-specific parameters, a user can simulate a number of relevant scenarios in a short time to see whether a specific implant in that location would potentially have risky zones. This can simplify calculations for more realistic clinical cases with 3D implant placement planning and outcome estimation.
The algorithm for the deployment of this method in practice might be as follows: • Examining the patient and determination of the bone/tissues status (quality) and parameters for the implantation site. • Selection of the implant type (dictated by the general considerations and safety).

•
Entering the patient's parameters and implant geometry into the 3D model. • Setting up "normal" (a nearly perfect osseointegration) or "worst case" (e.g., absence of cohesion between the screw or the abutment) scenarios and performing numerical calculations. • Post-processing the results, leading to the identification of potential risk zones and parameters.

•
In the case of need, adjustment or changing the implant type, shape, thread geometry, etc. and repeating the computation for the optimal outcome.
The limitations of the recent study are linked to likely over-generalized case of the tissues' properties: for the bone, the data scatter is lower, but soft tissues have much more variations. Authors might recommend that all the relevant tissues biomechanical properties should be tested in the conditions as possibly close to the real situation, and not in artificially created ones (even if such conditions were considered as a standard). In the present analysis, screw thread and details of the implant surface were not considered, as they require more detailed meshing and more computations. For realistic patient cases, the use of 1 4 of the full tissue can be too simplified as the mandible (as well as maxilla and other sites) are not fully symmetric at XZ and YZ planes (Figure 1).
For the analysis of literature data, there was little possibility to validate those test conditions and the data reliability coming from very different sources. Authors [63] have warned that in many published datasets, issues related to inertia effects (high frequencies) or instrument torque limits (low frequencies) are usually not sufficiently documented. Hence, even if instrument inertia is eliminated, the sample itself will always have finite