Numerical characterization of the load distribution in ball screws

Load distribution in ball screws is a representation of the ball contact stress, and it is fundamen-tal to understanding the behavior of these machine elements. This work aims to conduct a multi-variable analysis of the load distribution in ball screws. For this purpose, a numerical tool is developed for the generation and calculation of ball screw FEM models, which has been validated against the state of the art. Many diﬀerent design variables are studied in order to obtain a general characterization of the morphology of the load distribution in ball screws. The two most characteristic features, the non-uniformity at a local and global level are identiﬁed, along with as the possible causes of their appearance and the consequences that they may cause.

for the study of load distribution in ball screws. The experimental analysis of load distribution in ball screws is a complicated task due to the poor accessibility of the balls and the difficulty in measurement signals. These experimental tests with 2D prototypes enable obtaining results to validate theoretical models. However, they are unable to correctly represent 3D ball screws. Shimoda and Izawa [12] studied ball screw load distribution with half-pitch mechanics equations, which were complemented with their experimental analysis. Mei et al. [10] presented an analytical method to build a model to analyze the static load distribution by considering the Hertzian contact theory, nut/screw axial deformation and geometry errors. Xu et al. [16] considered the contact angle variation based on Wei and Lin's elastic deformation model [15]. These models show an uneven load distribution that decreases as the distance to the point of load application increases. However, they do not consider the effect of lateral deformation that causes non-uniformity at the local level. Later, Lin and Okwudire [7] proposed a model where the effect of axial, torsional and lateral deformations were considered, along with, Hertzian contact deformations and geometric errors of the ball screw components. It was shown that the lateral deformation directly affects the load distribution even if the applied load is purely axial.
This phenomenon and its effect combined with geometry errors was studied in more detail by Zhen and An [18]. Lastly, Zhao et al. [17] added the effect of the turning torque caused by assembly errors in their model. In addition, Wei and Kao [14] presented their model focusing on high load cases. Recently, Liu et al. [8] developed a model in which they consider the positioning of the nut along the shaft. All these works propose numerical models to obtain the load distribution in ball screws and each one, presents slight improvements with respect to the previous one. However, the analysis of the load distribution is carried out for particular geometry and load cases, taking into account the influence of the individual factors under study. The nature of the load distribution, with its general characteristic features, has not been studied in any case.
The FEM is a commonly used tool for modeling and optimizing ball-to-raceway contact based machine elements such as bearings and ball screws, as Lostado et al. [9] show in their double-row tapered roller bearing model. In this paper, a strategy to generate finite element models of ball screws is developed where an analytical model representing ball-to-raceway contact, is integrated in a commercial FEM software. Due to the reduction of the computational cost and the versatility of this method, these models can also be used as components of larger assemblies for specific applications. Then, once the proposed model has been validated with the contact ball bearing and also for ball screw results from literature, an analysis of the load distribution is carried out. Moreover, the influence of different ball screw variables, such as pitch, contact angle, number of start threads, ball size, slenderness and load arrangement and mounting options, are highlighted. Finally, a general characterization of the morphology of the load distribution, which could be used for optimization purposes, is developed.

Connector-based equivalent Contact Model (CCM)
Equivalent contact models have been developed for machine elements based on contact between raceways and balls by using non-linear elastic connectors. Dadalau et al. [4] carried out a linear guide system equivalent FEM model, where the balls were substituted with non-linear springs. Oyanguren et al. [11] adapted this model to ball screw thermo-mechanical models.
Connector-based equivalent contact model (CCM) is a ball screw FEM model that combines the precision and versatility of a high-order model with the efficiency of a low-order model by replacing the contact interactions with elastic joints based on an analytical study. With the aim of developing the model, the following assumptions are made: 1. Contact deformations between balls and raceways do not generate plastic strain. 2. Dynamic effects, such as centrifugal force and gyroscopic moment, are neglected. This assumption is valid when the rotational speed is low [13].  3. The load distribution analysis is made under static conditions. 4. Frictional forces in the contact are not taken into account. 5. While calculating the contact rigidity, the ball body's own stiffness is neglected because it is far larger than that of the contact itself.
This model allows calculating the load distribution and analyzing the influence that the different variables have on it. The basis of the model is described next and the study of ball screw behavior and variable analysis is carried out.

Analytical description of the interaction between raceway and ball surfaces
3D high-order ball screw finite element models, which include ball contact as is, involve very high computational requirements. The excessive mesh-fineness required to correctly represent the contact and the typical non-linearity associated with this type of problem makes it computationally very cost-expensive to analyze such high-order models. Therefore, a method capable of solving this limitation is proposed in the present work.
Each ball is replaced by a non-linear elastic connector that represents the rigidity of the contacts between balls and raceways (both the contact with the shaft and the nut in a single rigidity). As explained previously, the rigidity of the ball body deformation is neglected since it is much larger than the contact stiffness itself and does not affect the serial spring configuration. The contact rigidity is determined by the Hertzian elliptical contact theory provides the relationship between force (P ) and displacement (δ) in the contact between two curved bodies [6]: Where C = C s + C n represents the Hertzian contact stiffness for both contacts (screw-ball and nut-ball respectively) and it depends on the curvatures in contact and elastic material characteristics of both bodies. For each contact, C parameter is then calculated as: Where δ * is a dimensionless quantity obtained by solving the first elliptic integral [6]. E and ν are the Young's modulus and Poisson's ratio of each body. I and II represent ball and nut or ball screw. Σρ is the principal curvature sum of two contact bodies.
For the screw-ball contact: For the nut-ball contact: Where d w and d m are the ball and pitch circle diameters, respectively; f s and f n are the conformity ratio of the screw and nut raceways; α is the contact angle and λ is the nominal helix angle (without any manufacturing error).
Solving these equations, the relationship between contact load and displacement to be used as the connector element elastic behavior is obtained.

Implementation of CCM in a FEM model
Once the load-displacement relationship has been obtained, it is implemented in a finite element model using Abaqus as the FEM solver, where the ball's elastic behavior is represented. For this purpose, axial type spring connector elements (CONN3D2 in Abaqus) are employed and the non-linear elastic behavior is set according to analytical calculations.
Before introducing the connector elements, the high-order model is generated. As a high-order model, nut and screw components are precisely meshed with 3D hexahedron elements (C3D8). This allows for appropriately generating ball screw models with a wide range of particular specifications. Based on the parametric geometric definition, the mesh of both components is generated.
The axial type (1D) connector does not constrain any relative motion. The relative motion, u 1 , acts along the action line connecting the two nodes, measures the change in distance separating the two nodes and is defined as: where l is the actual distance between both nodes and l 0 is the initial one. The initial position of the constitutive nodes and the way they are attached to the raceways are relevant for obtaining accurate results, which implies having to contact kinematics under loaded conditions (Fig. 1). Understanding contact kinematics is relevant at this point.
Before any load is applied, the center of the ball O b and the center of the raceways O s and O n are in the same action line following the initial contact angle α (Fig. 1). It can be assumed that the component where the load is applied is the displaced one (in this case the screw), while the other stays fixed (in this case the nut) [3,18,17]. The three mentioned centers remain in a single line, with a common pressure angle for both contacts. This is acceptable at low rotation speeds where the centrifugal force and gyroscopic moment acting on the ball are negligible [13]. However, for high-speed operations, a model where the centrifugal force effect on the ball center position is considered would be recommended.
The positioning of the nodes constituting the connector in the centers of the raceways allows the variation of the contact angle to be correctly represented. Each node is attached to its corresponding raceway. The connector is therefore positioned in such a way that it works in tension (that is the extension of the connector a compressive ball contact force) and the contact kinematics is properly represented, as shown in Fig. 1   Once the constitutive nodes of the connector have been defined, the way the reference nodes are joined to the raceways needs to be studied, which is done using rigid coupling-type joints in which all degrees of freedom are restricted. The main issue is that the displacement of the raceways is redundantly taken into account. First, the contact deformation is included in the elastic behavior of the connector. Second, since the raceways are not rigid surfaces, the transmission of force in the joint causes an additional deformation in the raceway.
The solution to this issue is to apply a correction to the rigidity of the connector. The stiffness of the raceway is calculated in a parallel model for a load applied to the reference node, and the stiffness of the connector is corrected accordingly.
Lastly, the join area of the raceway that is attached to each connector should be determined. The previously explained correction applied to the rigidity of the connector makes this choice irrelevant. However, after several tests, it is observed that the most consistent results are obtained when the selected raceway join area is the largest possible. The structured build of the mesh implies that for a small portion, the union occurs in very few elements, causing a less stable calculation and an increase in boundary errors.
The implementation of the connector method in a finite element case is shown in Fig. 2. The connector element is defined between of the two reference nodes that are linked by a rigid joint to their corresponding raceway join area.
The whole process, from the ball screw definition to the results, is carried out by the interaction of Matlab and Python codes using Abaqus as the FEM solver. The flowchart followed to carry out this interaction and the model generation process is given in Fig. 3.

Validation process of the equivalent FEM model
The premise of CCM is that it should provide the same results as a high-order model. It is therefore necessary to validate the method by benchmarking it against a high-order model. Due to the amount of balls involved, high-order models of ball screws are computationally very cost-expensiv. The validation is therefore carried out using ball bearing models, which are equivalent in terms of contact definition and follow a very similar modeling process as ball screws. A reference bearing is chosen, and its rigidity is analyzed. Both models are defined and tested against axial and radial loads. The high-order model is generated with a complete mesh, which through iterations, presents a very fine mesh in the contact areas and a thicker one in the remote areas, as shown in Fig. 4.
In addition to this model validation, results from the CCM on ball screws are compared with cases presented in existing literature. Three load distribution study cases are chosen for this purpose: Okwudire's Case study 1 [7], Wei's heavy load ball screw [14] and Zhao's numerical calculation [17].

Ball Bearing model validation
The validation of the CCM against a high-order model is done using a commercial bearing as reference. The 7306-B-XL-2RS-TVP is used as the reference bearing and is calculated against purely axial load and combined axial and radial load. The results to be compared are the rigidity of the bearing in both directions. The applied load is up to 40 kN, which slightly exceeds the basic dynamic load rating of the bearing (35.5 kN). Figure 5 (a) shows axial load-displacement relationships.
Both ball high-order FEM model and CCM present very similar curves, with the maximum relative error lower than 0.5%. Radial rigidity comparison (b) presents analogous results, with a maximum relative error between both models lower than 0.7%.

Ball screw model validation
As a complement to any validation, a comparison is made against cases solved by other models of ball screws found in the literature. The analysis of this distribution enables knowing the load state of the ball screw, so that it is the result shown in all the case studies of ball screws.
It is known that the load distribution in ball screws generally presents a non-uniformity, wherein some balls support greater loads than others. The non-uniformity can be classified at local and global levels in  order to know and give value to this non-uniformity and its cause. The dimensionless parameter ratio r is set as an indicator of global non-uniformity of load distribution. This ratio represents the relationship between the maximum load and the average load. The second parameter, s, refers to the sinusoidal level in the local area. This parameter indicates the deviation of the load distribution curve from a second-order polynomial fitting curve acting as a baseline. Figure 6 shows a typical ball screw load distribution, with the mentioned baseline and the mean value represented. The two ratios r and s are obtained as: Both ratios will always be higher than 1 (except for an ideal case with uniform load on all balls, where r = s = 1), with their value increasing as the non-uniformity to which they refer increases. It should be noted that, according to the definition of the parameters, the s parameter should always be equal to or higher than the r parameter.
Following are the individual comparison cases. Okwudire's Case Study 1 [7] consists of a single nut ball screw with reduced dimensions (2 turns) to which an axial load is applied on the screw. The shaft is allowed to bend and, with it, lateral deformations can occur. The results of both models are shown in Fig. 7, where it can be seen that both curves maintain a certain similarity. Both the mean and the maximum values are similar and overall distribution coincides. However, there is a notable deviation in the center of the nut, with a peak in the CCM and a local minimum in Okwudire's model. The main reason for this deviation is probably the possible differences in the conditions applied in each case. According to the parameters, both r and s are small and very similar, which indicates that the non-uniformity is mostly local and is due to lateral deformations. The small dimensions of this ball screw limit the non-uniformity at a global level.
Wei's heavy load ball screw [14] analyzes a high load case using a nut composed of two independent reels with a total of 8 turns. The purely axial load is applied to the nut and the screw is fixed. Figure 8 shows the comparison of both models in this case. The results give load distribution curves that follow a  very similar trend. The maximum, minimum and average values as well as the position of peaks and valleys coincide significantly. The most remarkable deviations appear at the point close to the separation between reels. Again, these deviations are probably due to slight differences in the conditions of the case. In this case, the r parameter is higher than previous case, as shown by the downward form of the curve. The s parameter and, therefore, the non-uniformity at the local level, even though it is also high, loses relevance as it is much lower than the non-uniformity at the global level Lastly, Zhao's numerical calculation [17] studies the effects of turning torque on a medium nut (3 turns). Both the axial load and the turning torque are applied to the screw and the nut is fixed. The results in Fig. 9 show the comparison of both models. The curves follow a very similar distribution, except for the final balls, where small differences appear. The influence of the turning torque that causes the increase of the sinusoidal amplitude can be clearly appreciated. Proof of this is that the r and s values are very similar and high. This means that the non-uniformity of this case is mainly of local level, caused by the lateral deformations and the torque. The comparisons show that, in all cases, the results are similar in terms of load distribution values and trends; therefore, the CCM is very similar to those found in the literature.

Geometrical characterization of ball load distribution
Load distribution allows to better know the real load capacity of the ball screw. A study of the load distribution morphology and a variable analysis is carried out to determine the influence that some of the most used design parameters have on the distribution.
A reference ball screw is defined to perform the analysis. The selected ball screw has common dimensions for high-load and intermediate-speed applications. The duty cycle of these kind of ball screws include loading and unloading repetitive motions. The parameters of this reference screw are presented in Table 1. The parameters of the reference ball screw are maintained while the variables are analyzed individually. These variables are as follows: number of start threads, pitch, contact angle, ball size, slenderness and load arrangement.

Number of start threads
The number of start threads refers to the number of independent ball tracks sharing the same axial position. Multi-start ball screws are commonly used to increase the lead while maintaining the number of effective turns and, therefore, the length of the nut.
Single, double and less common triple-start ball screws are evaluated. The strategy is keeping the length of the nut, the total number of turns, the total number of balls and the pitch all constant. The lead increases proportionally with the number of start threads. The results are given in Fig. 10.
The curves show a very similar load distribution for all three layouts in general, where the similar r ratio values confirm this. The small differences do not follow any particular pattern and are caused by small deviations in the number of balls. However, the local level parameter, the so called s, provides more information. These differences among them mean that the sinusoidal level at the local scale decreases significantly as the number of start threads increases. This is because an increase in the number of start threads leads to an increase in the contact-load symmetry and, therefore, in the lateral stability of the ball screw. Consequently, the lateral deformations that cause this local non-uniformity are reduced.
This effect is concealed because the overall non-uniformity, which is reflected by the ratio, is much greater. This is due to the large dimensions and form factor of the ball screw.

Pitch
In single-start ball screws, the pitch refers to the axial distance when the helix completes one revolution. It relates the rotation speed of the ball screw to the linear feed rate of the nut and so is a main design parameter. The results in Fig. 11 show the influence of this variable.
The curves show that global non-uniformity increases with pitch, and the r values confirm this trend. However, for the range of values studied, which is quite wide for a ball screw with these characteristics, the difference in ratio is relatively small. The increase in pitch leads to an increase in the ratio, but with a low influence. The reduction in speed for the same feed rate would probably compensate for the increase in maximum loads on an overall dynamic analysis of the ball screw. Increasing the length of the nut and the distance between the balls leads to an increase in non-uniformity on a global level.

Contact angle
Contact angle α is the pressure angle in the ball to raceway contact. This angle reflects the ratio between the axial and contact load of the ball screw. The influence of the contact angle can be seen in Fig. 12 The most obvious conclusion is that the average contact load decreases as the contact angle increases. This is because the load is purely axial. The axial component of the contact load remains fixed and, therefore, the modulus of the contact load is a function of the angle. The non-uniformity of the load distribution remains the same in all three cases, as shown by the ratio values, so it can be said that the contact angle has a negligible influence on this point.

Ball size
The ball size is one of the design parameters and refers to the diameter of the balls inside the ball screw. Figure 13 shows the influence of this ball size on the behavior of the ball screw.
Changing the size of the ball while keeping the other geometrical parameters constant leads to a variation in the number of balls. The two figures show two ways to represent the same results.
The first one, Fig. 13 a, shows the load distribution in reference to the balls. The values indicate the specific force of each contact. In this case, the larger the ball size, the fewer the number of contacts and the greater the load per contact. However, if the value of the contact load is relative to the number of contacts and thus to the body of the nut, the contact force is equalized in all cases, as shown in Fig. 13 b.
Regarding the non-uniformity, the ratio values show that the size of the ball does not have a direct influence. Although in the first figure, the shape and the maximum values suggest that the bigger the size, the greater the non-uniformity, the ratios show that the non-uniformity hardly varies, mainly because of the increase in concordance of the average load.

Slenderness
Another aspect to be analyzed is the slenderness or the form factor of the ball screw. It has been seen in other results that the length of the ball screw is a factor that influences the load distribution; specifically, (a) Ball size influence. Load distribution for each ball. As the ball size increases, the number of balls decreases, so the load per ball has to increase (b) Ball size influence. Load distribution for the nut body. Even if there is different a number of balls for each case, the nut length remains equal. Fig. 13: Ball size influence the length of the nut is related to the non-uniformity at the global level. However, to better analyze this aspect, the slenderness of the ball screw is studied instead of focusing only on the length. In this manner, the overall ball screw dimensions are left out and the focus is on the form factor.
The slenderness, in this case slr, is defined as the ratio between the length of the nut l 0 and the diameter of the pitch D pw . The geometrical data is kept constant and the number of turns is increased from 2 to 16. The reference ball screw of the previous sections corresponds to the 10 turns one (slr = 3.086), with the difference that, in this case, it has only one reel while in the previous cases, it had two. Figure 14 shows the load distribution according to the slenderness. The most obvious conclusion is the non-uniformity increases as slenderness increases. The r ratio increases as slr increases, while s remains almost constant. The area closest to the point of application of the load suffers a very high force and the force decreases as it moves away from that point. Therefore, it is shown that as the slenderness increases, the non-uniformity also increases.

Load arrangement and mounting options
An aspect to consider when analyzing load distribution is the load arrangement and mounting options of the ball screw. Four different cases are analysed. In the first two cases, Fig. 15 a and b, the nut is compressed while the shaft is in both tension and compression. In the two second cases, the nut is in tension and the shaft is in tension and compression. This way, all the possibilities of load arrangement in the axial direction are analyzed. The results of the load distribution for the four cases are shown in Fig. 16. The curves show two different behaviors. The first one is given for the cases C-t and T-c and is characterized by a very severe non-uniformity (r close to 2). This is due to the fact that the load on both components is applied on the same side, with the other side being free. This means that the area near the application side suffers a much greater strain than the opposite side. The second is given for the cases C-c and T-t, and the non-uniformity is considerably reduced (r close to 1.5). In these cases, the load is applied Fig. 16: Load disposition influence. C means that the component is being subjected to compression loads while T means tension loads. Capital letters refer to the nut and lower-case letters to the screw. The applied load in each case is axial and to the right.
on both sides (on the shaft on one side and on the nut on the opposite side). Thus, the ends have higher loads than the central area, with a much smaller difference than in the previous cases. The sinusoidal ratio s remains very similar in all cases Therefore, it can be stated that the arrangement of the boundary conditions plays a very important role in the load distribution.

Discussion
The results shown indicate that the load distribution in ball screws can be very variable and depends on many different factors. However, it is possible to identify a generic trend being followed to a greater or lesser extent by all the analyzed cases. As mentioned in this paper, the load distribution per ball is characterized as a variable distribution.
The r and s ratios, in fact, are two ways to analyze the non-uniformity of load distribution in ball screws. Analyzing the results, it can be seen that the different variables or specific characteristics of the case studies can affect these ratios. It is possible to deduce the reasons for the non-uniformities and, therefore, to obtain a characterization of the load distribution that allows for a better understanding of the behavior of these machine elements.

Global level non-uniformity
Non-uniformity at a global level refers to the deviation between the ball loads depending on the axial position they are in along the nut. This non-uniformity means that there are areas of the nut that suffer higher-than-average loads and are therefore more prone to wear or surface fatigue.
Non-uniformity at a global level is found to a greater or lesser extent in the load distribution of all the cases analyzed. It has been seen that there are variables and factors that affect this non-uniformity in greater measure than others.
Load arrangement and slenderness are shown as determining factors in this aspect. Slenderness is a parameter that indicates the relative length of the ball screw (slr = l 0 /D pw ). Since this non-uniformity is subject to the axial positioning of each ball along the nut, it is logical that the relative length of the nut has a direct effect. The mounting arrangement determines the position of both the active and reactive axial forces. The load level of the balls increases according to the proximity to these points of application of the loads because of strain concentrations.
The influence of these two factors suggests that the main reason for this overall non-uniformity is the axial deformation of the components. The balls that are closer to the point of application of the load receive the highest contact load intensity. As the position of the ball advances along the nut, the axial deformation of the components reduces the intensity of the contact loads. To demonstrate this, the load distribution is analyzed as a function of the stiffness of the components. The previous ball screw is taken as a reference and Young's elastic modulus is taken as a variable. The reference ball screw has the mechanical properties of the steel. The rest have elastic modulus 10, 100 and 1000 times higher. the results are given in Fig. 17.
Through the r ratio values, it can be seen that as the stiffness of the materials increases, the overall non-uniformity decreases. The materials in this analysis, with the exception of reference steel, have no real physical significance, however, they serve the purpose of studying the influence of stiffness on load distribution. It is therefore proven that the reason for the appearance of an overall non-uniformity is mainly due to the axial deformations suffered by the components of the ball screw.

Local level non-uniformity
Non-uniformity at a local level refers to the deviation between the ball loads depending on the radial position they have in the helix. This non-uniformity means that there are periodical load peaks at each pitch of the helix.
All the tested configurations have, local non-uniformity to a certain level. In most of the cases, the local non-uniformity, quantified by the s ratio, does not practically depend on the parameter studied. Only the number of start threads and the rigidity of the material, analyzed in the previous point, seem to have a direct influence on the non-uniformity at the local level.
This non-uniformity follows a cyclical pattern, with a period equal to the pitch of the helix. The value of the load therefore depends on the radial position of the ball in the helix. These factors allow us to determine that the reason for the appearance of non-uniformity at the local level is the lateral instability caused by the helix itself. However, as a demonstration, an analysis is carried out in which radial loads are introduced at different levels. In addition to the reference case, study cases with radial loads of 5% and 10% of the applied axial load are analyzed in Fig. 18.
Graphically, it can be seen that local non-uniformity increases as the radial load applied increases. Where data of the s ratio shows the same trend, local non-uniformity is a consequence of a relative radial displacement between the nut and axis. Under the influence of a radial load, it is the bending force itself that causes this displacement, while under the only influence of axial load, it is the helix's own ball arrangement that causes this lateral instability.

Conclusions
This paper analyzes the capacity of the CCM model to generate, simulate and analyse finite element models of ball screws. The CCM is used to study the load distribution on ball screws in order to develop a general characterization of its morphology that allows a better knowledge of the behavior of these machine elements.
The importance of the characterization of the load distribution in ball screws lies in the fact that it offers great insight into the behavior of the screw. The load distribution enables knowing the load level of each ball and, what is more interesting, the load state of the nut. This load state, directly affects properties such as rigidity, precision and the ball screw's service life.
The principal conclusions drawn from this study are the following: 1. The CCM connector-based model is equivalent to high-order FEM models, with a highly reduced computational cost. The great versatility and integration capacity in assemblies of larger FEM models gives it an advantage over lower-order numerical models. 2. Ball screw load distribution depends mainly on its design parameters and the specific conditions of the load case. In any case, the load distribution is far from being uniform. This non-uniformity appears in two different ways: non-uniformity at a global level and non-uniformity at a local level. 3. The global non-uniformity is due to the axial elastic deformation of the ball screw components and depends mainly on the slenderness and load arrangement of the ball screw. 4. The local non-uniformity is due to the lateral deformation caused by the inherent lack of symmetry of the ball screw helix. 5. Both non-uniformities lead to the presence of areas with higher-than-average loads that directly affect the service life of the ball screws.

Acknowledgments
The technical and financial support of Shuton, S.A. is greatly acknowledged.

Ethical Approval
Not applicable

Consent to Participate
Not applicable

Consent to Publish
Not applicable

Authors Contributions
-LS performed the conception and design of the study and drafted the manuscript.
-AO contributed in the conception and design of the study and revised the manuscript critically for important intellectual content. -JL contributed in the conception and design of the study.
-MI contributed un the acquisition of data.
-AA contributed in the analysis and interpretation of the data.
-IU contributed in the conception and design of the study and revised the manuscript critically for important intellectual content.

Funding
Not applicable

Competing Interests
The authors declare that they have no known competing interests.

Availability of data and materials
Data supporting the results reported in the article can not be shared publicly. However, such data is available and can be obtained after consulting with the corresponding author.