Introduction

Harris lines (HL; also known as “growth arrest lines” or “transverse radiopaque lines”) are horizontal sclerotic lines formed in the metaphyseal or diaphyseal part of long bones, usually visualized using X-ray images (Nowak and Piontek 2002a), although also visible in computed tomography (Chauveau et al. 2016; Primeau et al. 2016), magnetic resonance imaging (Laor and Jaramillo 2009), or in histological slides (Miszkiewicz 2015). Most often they form on the tibia and therefore, it is the most frequently used bone for the study of HLs (Papageorgopoulou et al. 2011); nevertheless, they can form on any endochondral bone (Scott and Hoppa 2015).

HLs are a trace of a temporary bone growth arrest. Undisturbed secondary endochondral ossification has been presented and described in Fig. 1. Disruption of the process of chondrocyte proliferation or calcification means that the growth plate remains impenetrable for osteoblasts, which in turn cause mineralization along the horizontal chondrocyte layer at the end of the growth plate, forming the “primary stratum” perpendicularly to the long axis of the medullar cavity. Prolonged periods of growth arrest result in the thickening of the primary stratum, which leads to the deposition of the Harris line (Scott and Hoppa 2015). However, according to a different model, a HL forms when growth is re-established; increased cartilaginous proliferation and osteoblasts activity contribute to the thickening and forming of the transverse line, suggesting HLs are growth restart rather than growth arrest lines (Sajko et al. 2011).

Fig. 1
figure 1

Endochondral bone growth and the Harris lines deposition theory. Drawing on the left side presents undisturbed bone growth. (1) Chondrocytes (cartilage cells) from the resting cartilage zone move towards the bone diaphysis. (2) Proliferating chondrocytes form columns, proceeding towards the diaphysis. (3) Chondrocytes differentiate into hypertrophic chondrocytes which produce alkaline phosphatase, osteopontin, BSP, Osx, Bglap (osteocalcin), and Runx2 (Park et al. 2015), which induce calcification. (4) Hypertophic cells undergo apoptosis and calcify. (5) Osteoprogenitor cells differentiate into osteoblasts, causing the initiation of osteosis on the basis of calcified cartilage (Mescher 2013). The disruption of chondrocytes proliferation or differentiation (presented on the right side) makes the growth plate impenetrable for osteoblasts, which in turn begin to create a mineralized layer along the horizontal layer of chondrocytes at the end of the epiphyseal plate; they form a primary stratum perpendicular to the long axis of the medullary cavity. Prolonged periods of growth arrest result in the thickening of the primary stratum, which leads to the deposition of a Harris line (Scott and Hoppa 2015)

Among the factors that may lead to a temporary arrest of bone growth (and thus—to HL deposition), the most commonly mentioned are nutritional disorders (malnutrition, protein, vitamin, and mineral deficiencies), a history of smallpox, pneumonia or other diseases, food poisoning, or alcohol abuse (Nowak 1996). Among newer concepts, HLs can be linked to child abuse (Ross and Juarez 2016), systemic-onset juvenile idiopathic arthritis (Sifuentes Giraldo et al. 2016), tortures (Traczek 2017), psychosocial short stature (Kaspar Hauser syndrome; Khadilkar et al., 1998) osteopetrosis, hyper/hypoparathyroidism, sclerosing, spondylosis, radiation exposure, Cushing’s Syndrome, rickets, avascular necrosis, osteoporosis, congenital syphilis, Paget’s Disease, leukemia, scurvy, bone fracture (Sajko et al. 2011), and probably many more diseases and detrimental conditions. The wide spectrum of possible etiologies causes HLs to be considered as an indicator of an unspecific physiological stress or difficult living conditions rather than a specific diagnostic marker (Papageorgopoulou et al. 2011).

However, using HLs even in such a wide context is often questioned. In an in vivo experiment on rabbits, the authors noted that more HLs were formed during periods of rapid growth than in the period of nutritional stress (Alfonso-Durruty 2011). There was no statistically significant correlation between bone length and morphology and the occurrence of HL (Mays 1985; Nowak and Piontek 2002b). Moreover, numerous researchers have failed to find correlations between the age of enamel hypoplasia (EH) and HL formation, suggesting that HLs are a physiological rather than a pathological phenomenon (Piontek et al. 2001; Beom et al. 2014; Geber 2014; Zapala et al. 2016; Krenz-Niedbała 2017). Due to the similarities between the somatotropin (growth hormone) secretion curves, the long bone growth curves, and the distribution of HLs, it was suggested that HLs are associated with physiological periods of faster and slower growth (Papageorgopoulou et al. 2011).

Nevertheless, the lack of correlation between HL and EH may result from the disappearance of early-formed HLs as a result of bone remodeling at a later age (Mays 1995). In turn, the formation of HLs in a period of rapid growth corresponds to the theory of HLs as growth recovery lines. Moreover, if one consider the formation of HLs as a trace of physiological processes only, it is necessary to explain the differences in the frequency of HL occurrence between individuals and populations found in numerous studies (Piontek et al. 2001; Beom et al. 2014; Geber 2014; Zapala et al. 2016; Krenz-Niedbała 2017). In addition, even if it will be concluded that HLs are poorly suited for determining the health of individuals (as it may be a result of non-pathological processes), they can still be useful for determining the overall health at the population level (Nowak and Piontek 2002b): assuming that the occurrence of HL unassociated with any pathologies is (more or less) constant, statistically different frequency of HL between the compared populations would result from different occurrence of HL associated with pathologies or physiological stress.

Regardless of the chosen HL formation model, it can be assumed that a HL is formed more or less in the time following the appearance of a stressor hindering growth to length. It is easy to deduce that the HL will be located where the growth plate was located when the growth was stopped/restored. The position of the HL is therefore related to the period of incidence of physiological stress inhibiting bone growth, which makes it possible to estimate the age at which the subject was exposed to it (Papageorgopoulou et al. 2011). Such information can be valuable in the study of archeological populations. The mere counting of HL allows to compare the level of physiological stress between specified groups. Calculation of the age of HL formation enables more comprehensive statistical approach. It can be concluded, which period of childhood or puberty in the study population was most abundant in adverse effects (Zítková et al. 2004; Krenz-Niedbała 2014; Boucherie et al. 2017)—and compared between specific groups (Jerszyńska and Nowak 1996). It is additional easily obtainable data that could improve the description of life of archeological populations. Therefore, various methods have been developed to determine the age of HL deposition (Allison et al. 1974; Hunt and Hatch 1981; Clarke 1982; Maat 1984; Hummert and Van Gerven 1985; Byers 1991).

In this review, six methods for calculating the age of HL origin will be presented: the Allison et al. (1974); McHenry and Schulz (1976), Hunt and Hatch (1981), Clarke (1982), Hummert and van Gerven (1985), Maat (1984) and Byers (1991) methods. To the best of the authors’ knowledge, these are all known methods. In addition, a modification which enables calculations on non-adult bones will be proposed to the last of them. All tables required for calculations were included in Supplementary File 1, and all necessary measurements for each method—in Fig. 2. Summarized comparison of key features of all methods is presented in Table 1. Method output and detailed comparison are presented in the “Comparison of all methods” section.

Fig. 2
figure 2

All measurements needed for individual methods, on the example of the tibia. For more details on each method, please refer to the main text. a Allison’s method. b Hunt and Hatch’s method. The DPOC measurement needs to be referred to the appropriate table c Clarke’s as it is made for other methods. The Dys/T ratio needs to be referred to the appropriate table. d Byers’ method. Both T and Dys or P (depending on the distal or proximal position of the HL) need to be used in the formulas presented in Table 2. e, f Maat’s method. The DPOC/DEPI ratio needs to be referred to the appropriate table. g Hummert and van Gerven’s method. DPOC/Dx needs to be referred to the appropriate table. All necessary tables can be found in Supplementary File 1

Table 1 Summarized characteristics of the methods reviewed. For more details, refer to the appropriate sections of this paper

However, before the presentation of methods, one important issue must be raised: what should be considered a Harris line?

What line should be considered a Harris line?

After death, HLs remain unchanged; they can be observed on the remains of more than 45,000 years (Nowakowski 2018). However, during a lifetime, the bone undergoes constant remodeling processes which results in a possible disappearance of the HL (Papageorgopoulou et al. 2011). This hinders the interpretation and use of HLs in archeological studies of the population, especially for the remains of older individuals.

Usually, when examining the bones of adult individuals, clear lines running along the entire length of the bones will not be found, but thin remains of distinct HLs will be visible instead. There are different approaches to what should be considered a HL. Some authors recommend taking into account lines of at least 25% of bone width (Goodman and Clark 1981; Piontek et al. 2001), others: at least one third (Alfonso et al. 2005; Chauveau et al. 2016) or half (Hummert and Van Gerven 1985; Scott and Hoppa 2015) of the bone. To the best of the authors’ knowledge, there are no unambiguous guidelines suggesting the optimal solution in this matter.

Another approach is the arbitrary categorization of HLs; Maat (1984) categorized HLs into three types: weak (visible only upon careful inspection), medium, and strong, recognizing the third type as the most reliable determinant of physiological stress. However, to the best of the authors’ knowledge, this division has not been applied in other studies.

What is more, the classification and counting of the Harris line is also hindered by possible discrepancies between observers (MacChiarelli et al. 1994; Grolleau-Raoux et al. 1997). A possible solution to this problem is the proposed semi-automatic Harris line detection system: HL-tool (Suter et al. 2008; Papageorgopoulou et al. 2011). However, although undoubtedly useful, reliable, and giving repeatable results, this program works quite slowly, especially on older computers. The program, tested by the authors of this review on a computer with the Windows Vista™ operating system, an Intel®Core™ 2 Quad CPU 2.5 GHz processor, and 4 GB RAM, needed 4 to 8 min for one bone, excluding the preparation of X-ray images for analysis.

Hence, children’s bones are definitely worth the attention in studying HLs. Because the subjects are young, the HL is more likely not to be resorbed. Clear, pronounced HLs are also less probable to give different results between different observers. Numerous authors reported increased number or visibility of HL on juvenile bones (Piontek et al. 2001). More surprising is the small number of methods which enable the determination of the chronology of HLs on juvenile bones. Moreover, current methods allow to estimate the age of HL formation in the distal part of the juvenile tibia only. Therefore, it was necessary to develop a new method which would have wider application—and which is presented in the “Byers’ modified formulas for juvenile bones” section.

The Scott and Hoppa study (Scott and Hoppa 2015) is also worth mentioning in the context of HL detection. The authors compared the visibility of HLs on bone X-ray images taken in the anterior-posterior view (A-P) and medial-lateral (M-L) view. Although A-P images are usually used for observation, the authors achieved better visibility with the M-L plane.

The Allison et al. (1974) and McHenry and Schulz (1976) methods

These are the two oldest of the methods discussed; the former was used to determine the chronology of formation of HLs on the tibia, the latter—on the femur. Due to their high similarity and identical assumptions, they can be treated as one method.

Basic assumptions

The Allison method assumes that (a) the average bone length at birth is 90 mm, (b) the bone grows by the same length every year, and (c) 40% of growth occurs on the distal part of the bone, and 60%—on the proximal part.

Necessary measurements

To calculate the age of HL deposition using the Allison method, one needs to: (a) subtract 90 mm from the total bone length (T), (b) multiply the result by 60% to calculate the increase in the proximal part (or by 40% for a HL found on the distal part), and (c) divide the obtained length into 16 equal pieces, each correspondent to 1 year. Next, the position of the Harris Line should be compared to the obtained “ruler.”

The McHenry method introduces one modification; it focuses on HLs on the femur, in which a 71% increase in length occurs on the distal part and a 29% one on the proximal part. Consequently, the percentages by which the bone length is multiplied change accordingly. However, McHenry does not specify in his paper if the bone length is also shortened by 90 mm or any other value.

Summary of the method

The method is depicted in Fig. 2a. The authors themselves note that the assumptions are not very precise; the length of the bones after birth can be different than 90 mm (Scheuer and Black 2000), and bone growth varies with age (Hoppa 1992; Ruff 2003). In addition, the total bone length is measured, ignoring the fact that HLs form at the epiphyseal plate, not at the end of the epiphyses (Ameen et al. 2005). Nevertheless, despite the obvious lack of precision, the method enables a rough categorization of HLs into those that were created in the early or late period of growth. In addition, it is relatively simple and fast to use and does not require a comparison of the results with tables or growth curves.

Hunt and Hatch’s (1981) method

This method is based on a mathematical model of bone growth. Although the calculations themselves may seem difficult to understand, they are actually quite simple to use, especially when using the spreadsheets in Supplementary File 1. In short, the method requires the following steps: (1) creating growth curves for the study population, (2) creating curves showing the dependence of distance between the HL and the primary ossification center (POC) and the age of HL deposition, (3) determining the POC of the examined bone, and (4) measuring the distance from the HL to the POC and comparing the result with the curve defined in point 2.

Creating growth curves for the study population

The formulas were derived from long-term studies on the growth of long bones in the population of Denver, Colorado, collected by McCammon (1970). Hunt and Hatch derived the growth curves described by the formula:

$$ \overline{g}(t)=\frac{a_1}{1+{e}^{-{b}_1\left(t-{c}_1\right)}}+\frac{a_2}{1+{e}^{-{b}_2\left(t-{c}_2\right)}} $$

The \( \overline{g}(t) \) function describes mean diaphyseal bone length in time t. The first part of the equation (variables with 1 in the subscript) describes the growth in the prepubertal term, the second part (variables with 2 in the subscript)—in the adolescent term. Variables a are the upper limit of bone growth, b—the maximum slope at the point of inflection of the sigmoid left logistic function, c—age in the period of the fastest growth, and t—age in years (from 1 to 18).

Understanding how to derive this formula may be quite challenging, but it is not required for its proper use, so it will be omitted from this review; for more details, please refer to the original publication by Hunt and Hatch.

This method allows to create bone growth curves for the humerus, radius, ulna, fibula, femur, and tibia, and allows to estimate the age of HL formation in the last two types of bones. An exemplary growth curve for the tibia is shown in Fig. 3a. The remaining curves (with all necessary calculations) are included in Supplementary File 1.

Fig. 3
figure 3

Example growth curves derived from Hunt and Hatch’s study. a Mean tibia growth. b Curve for the estimation of the age of HL deposition

Next, the bone growth curves should be fitted to the study population. To do so, it is necessary to determine the mean diaphyseal length of distinct bones of adults in this population. Then: (a) measure the total length of the bones, (b) subtract the average epiphysis thickness (determined by the authors—see the table in Supplementary File 1) from it to receive the variable \( \overline{y} \). The next step is to modify the \( \overline{g}(t) \) function using the average diaphyseal length obtained from the population under study. The equation looks as follows:

$$ \overline{z}(t)=\frac{\overline{y}\ }{a_1+{a}_2}\times \overline{g}(t) $$

Now, simply substitute numbers from 1 to 18 for t, and use the results to create a specific bone growth curve. This curve can be used to estimate the age of a particular individual.

Curves representing the interdependence of the distance between the HL and the primary ossification center (POC) and the age of HL deposition

Based on the research of Anderson et al. (1963), the authors of this method assumed that 29% of femoral growth occurs at the proximal end and 71% at the distal end. For the tibia, the proportion of proximal and distal growth is 57 and 43%, respectively. In order to create curves that will determine the position of the HL with respect to the POC depending on the age of birth, the following formulas should be used:

$$ {\displaystyle \begin{array}{c}{0.29}^{\ast }\ \overline{z}(t)\hbox{--} \mathrm{for}\ \mathrm{the}\ \mathrm{proximal}\ \mathrm{part}\ \mathrm{of}\ \mathrm{the}\ \mathrm{femur}\\ {}{0.71}^{\ast }\ \overline{z}(t)\hbox{--} \mathrm{for}\ \mathrm{the}\ \mathrm{distal}\ \mathrm{part}\ \mathrm{of}\ \mathrm{the}\ \mathrm{femur}\\ {}{0.57}^{\ast }\ \overline{z}(t)\hbox{--} \mathrm{for}\ \mathrm{the}\ \mathrm{proximal}\ \mathrm{part}\ \mathrm{of}\ \mathrm{the}\ \mathrm{tibia}\\ {}{0.43}^{\ast }\ \overline{z}(t)\hbox{--} \mathrm{for}\ \mathrm{the}\ \mathrm{distal}\ \mathrm{part}\ \mathrm{of}\ \mathrm{the}\ \mathrm{tibia}\end{array}} $$

Numbers from 1 to 18 should be substituted for t. The results can be used to create appropriate curves (Fig. 3b).

Determining the primary ossification center (POC)

As described above, 29% of femur growth occurs in the proximal end and 71%—in the distal end. Accordingly, the POC is located at 71% of the length of the bone if measured from the distal end or at 29% of the length if measured from the proximal end. For the tibia, the proportion of proximal and distal growth is 57 and 43%, respectively, so the POC will be at 43% of the length of the bone, measured from the distal end.

The last stage of the procedure is to measure the distance from the HL to the POC and match the obtained measurement to the growth curve for a given bone and its proximal or distal part.

Summary of the method

The Hunt and Hatch method allows to estimate the age of HL formation in the bones of adult individuals, both for the distal and proximal parts of the tibia and femur. It is based on a fairly accurate model that can be adapted to different populations. The accuracy of this method and the wide range of applications are undoubtedly its greatest advantages.

Nevertheless, the method is quite sensitive to individual deviations; for particularly high or low individuals, the method will provide an underestimated or overestimated age of HL formation, respectively. A significant disadvantage of this method is its laboriousness; both understanding and using the Hunt and Hatch method takes a long time, and calculations can discourage less mathematically skilled anthropologists. However, supplementary materials attached to this publication should significantly facilitate the use and understanding of this method.

Clarke’s method (1982)

Basic assumptions

Clarke’s method was an attempt to improve the methods of Allison and McHenry. It does not assume even bone growth throughout the period of puberty and adolescence. Instead, the author created a table based on the Gindhardt study (1973), which reflects the variable bone growth; it reflects the ratio of bone length in an immature individual of a specific sex and age to the length of an adult bone (Supplementary File 1).

The reference table was adapted to calculations for HLs on the distal part of the tibia; according to Anderson et al. (1963) study, it is assumed that 43% of the tibia’s growth at the distal end and 57% at the proximal end. Therefore, each of the results in the table was multiplied by 43%.

Necessary measurements

Application of the method requires the following steps:

  1. a)

    Measuring the bone length (T), determining the sex of the subject

  2. b)

    Measuring the distance from the distal end of the bone to the HL (Dys)

  3. c)

    Calculation of Dys/T × 100% equation and comparison of the result with the appropriate table

Necessary measurements are presented in Fig. 2c.

Summary of the method

The method is quite fast as it requires only two simple measurements. However, it has some serious imperfections: it takes into account the total length of the bone, omitting the fact that the line is formed within the growth plate, not at the end of the epiphysis, which generates a constant error. In addition, the method can only be used for the distal part of the tibia.

Maat’s method (1984)

This method makes it possible to determine the age of HL formation at the distal end of the tibia. The author created a table (based on the research of Maresh (1955)) describing the length of the tibia diaphysis and its distal part depending on the age and sex, as well as the ratio of the distal length of the diaphysis at a given age to its average length in adults (Supplementary File 1).

Necessary measurements

Maat’s method requires the determination of: (1) the primary ossification center (again, 43% of tibia length measured from the distal articular surface), (2) the POC’s distance to the epiphyseal plate, (3) the distance between the primary ossification center and the HL. Then, divide measurement 3 by measurement 2 and refer the result to the table prepared by the author.

Summary of the method

This method is relatively easy to use, it requires few measurements and simple calculations, and it cannot be accused of inaccuracies; it includes the thickness of epiphyses; the comparable tables are based on reliable growth curves and it uses ratios instead of absolute measurements. However, this method is limited only to the distal part of the tibia. It can also be used to examine children’s bones—if the age of the child was determined, the table can estimate the length of the bone diaphysis of a “would-be” adult and use this hypothetical length in the calculations.

The Hummert and van Gerven method (1985)

This method is designed mainly to assess the age of HLs in the bones of children and adolescents. It was dedicated to the assessment of HLs in the medieval population of Kulubnarti (Sudan). Data from this population was used to create a growth chart employed further in this method, assuming that the bone growth in this population may differ from the European or American populations, on which the growth curves of the previous methods were based.

Basic assumptions

Calculations on non-adult bones require the determination of the approximate age of the examined individual. Based on the previously collected data (Hummert 1983), the authors created a table that defines the ratio of the mean distal part of the tibia (again, 43%) in the subsequent developmental years to the length of the distal part of the tibia in the evaluated individual of known age (Supplementary File 1).

Necessary measurements

The Hummert and van Gerven method is another one that requires determining the POC according to the same assumptions as in the Hunt and Hatch (1981) and Maat (1984) methods, 43% length of the bone, measured from the distal end. The second measure is the distance between the HL and the distal end of the bone. Then, the second measurement should be divided by the first one and referred to the table created by the authors.

Summary of the method

The authors recommend not to rely on age estimation based solely on the bone length but also on other age determinants. They estimated the developmental age basing on sequences of dental formation and eruption (Uberlaker 1978). Among newer methods, age can be estimated based on the degree of epiphyseal union (Scheuer and Black 2000), the frontal sinus development (Moore and Ross 2017), and measurements of the ilium (Corron et al. 2017) or girdle bones (Cardoso et al. 2017). One quite obvious disadvantage of the Hummert and van Gerven method is the restriction to the distal part of the tibia. Moreover, the table obtained by the authors refers to the medieval population of Kulubnarti, so it is necessary to create a new table for a population with different dynamics of bone growth.

Byers’ method (1991)

The last method is both easy to use and accurate. It is inventive, based on logical assumptions, and it allows to determine the chronology of the HL on four types of bones (tibia, femur, radius, and humerus). Since the authors of this review propose a modification to these formulas, it is necessary to introduce the concept of this method and show its derivation. However, as with the Hummert and van Gerven models, understanding the derivation of the formulas is not required to use this method effectively.

Basic assumptions

Byers based his considerations and calculations on observations already included in the presentation of the previous methods: (a) Harris lines form in the place of the epiphyseal plate, (b) bone growth dynamics vary with age, but, based on literature data, tables/growth curves describing the growth of individual bones over time can be designed, and (c) bone growth differs at the distal and proximal parts of the bone.

Based on the works of Maresh (1955) (in which bone measurements of children aged 10–12 included the measurements with the epiphyses and without them), Byers calculated that the epiphyses constitute 8% (± 0.68%) of the total length of the humerus and the radius, 11% (± 0.29%) of the femur, and 13% (± 0.33%) of the tibia for both sexes. The growth tables were based on studies on children from Boston and Yellow Spring (Anderson and Green 1948; Maresh 1955; Anderson et al. 1963; Gindhart 1973). In addition to the aforementioned proportions of growth for the distal and proximal parts of the femur (29%/71%) and the tibia (43%/57%), he took into account the proportions for the humerus (81% on the proximal end) and the radius (75% on the distal end).

Formulas derivation

All Byers’ formulas boil down to determining the ratio of bone length at the time of HL formation (Td) to the total length of the subject’s bone (T), which will then be referred to the growth table for individual bones. It can be defined with the Td/T × 100% formula.

While T is a trivial measurement, Td is harder to determine; the location of the HL tells only where the epiphyseal plate was located at the time of deposition. Therefore, it was necessary to develop a formula allowing to estimate Td based on the distance of the distal (Dys) or proximal (P) HL to the bone end. In this example, it is assumed that there is a HL on the proximal part of the tibia (Fig. 4):

  1. a)

    The distance from the HL to the proximal end (P) is known. Because of a 43% tibia growth on the distal part, Dys will be 43/57 × P or 0.75P. The length of the diaphysis (bone shaft) at the time of deposition (Sd) will be T − (Dys + P) or T − 1.75P.

  2. b)

    Since Sd was calculated, Td (total length of bone during HL deposition) needs to be calculated. Since the length of the epiphyses is 13% of the whole tibia, it can be assumed that: Td = Sd + 0.13Td. Then, Sd = 0.87Td and further: Td = 1.15Sd. Therefore, Td can be represented as Td = 1.15(T − 1.75P).

Fig. 4
figure 4

a Explanation of Byers’ method on the example of the tibia with one HL on the proximal part. The aim of this method is to determine the Td/T ratio. T can simply be measured, while Td needs some calculations. Subtract T − (Dys + P) to get Sd, which in turn needs to be multiplied by 1.15 to receive Td. Since Dys is unknown, it is necessary to estimate this length by multiplying P by 0.75. b Explanation of the modification of Byers’ method for calculations on non-adult bones. The lack of epiphyses prevents obtaining T and P measurements; therefore, this method focuses on the diaphyses, not the entire bone length. The aim of this method is to determine the Sd/S ratio and it uses similar measurements as the original Byers’ method

The formulas shown in Table 2 can be derived analogously.

Table 2 Original Byers equations and theirs modification

Necessary measurements

Only two measurements are required: (1) the total length of the bone (T) and (2) the distance from the HL to the nearest end of the bone (Dys or P). Both values can be used in appropriate formulas and the result can be compared with the values in the tables (Supplementary File 1).

Byers’ modified formulas for juvenile bones

Byers’ formulas seem to have only one disadvantage: they cannot be used to calculate the age of HL deposition in juvenile bones. However, the authors propose here a simple modification to Byers’ method which will allow to make such calculations.

In children and adolescents, it is impossible to obtain the total length of the bones (due to lack of epiphyses), but only the length of diaphysis. The modification to the formulas is based on the assumption that \( \frac{T_{\mathrm{d}}}{T}=\frac{S_{\mathrm{d}}}{S} \), where Sd is the length of the diaphysis at the time of formation of the HL, S—length of the diaphysis of the evaluated bone; de facto the total available bone length. Therefore, the main aim of the formulas is changed to Sd/S × 100%. Again, S is a trivial measurement, and Sd needs to be calculated.

Because of the lack of epiphyses, it is not possible to measure P or Dys (namely distance between the HL and the end of bone). It is possible, however, to measure P′ or Dys, which will be the distance between the HL and the end of the diaphysis. Therefore, Sd = S − (P′ + Dys). It is assumed that the P′/Dys proportions are the same as P/Dys for all bones; therefore, all P and Dys in the original Byers’ formulas can be replaced with P′ and Dys’.

There is also a need to modify the following elements of Byers’ formula: replace T with S and then remove the part of the formula “adding” the assumed length of the epiphyses to the calculated length of the bone shaft at the time of HL formation. The modified formulas are presented in Table 2.

Finally, it is required to compare the result with the modified Byers’ table (Supplementary File 1) and read from it the approximate age of HL formation. Similarly to the Hummert and van Gerven methods, in order to be able to use it, the approximate age of the subjects should be determined beforehand.

Comparison of all methods

Models

To compare the outcomes of different methods, two models were created. The first consists of one enormously HLs-rich male tibia, which was scaled to the following lengths, 44, 39.2, and 34.5 cm (Fig. 5). Then, the age of formation of every single HL was calculated with all applicable methods. For the purpose of the Hunt and Hatch method, 39.2 cm is assumed as mean male tibia length.

Fig. 5
figure 5

Bone models for the comparison of methods. For the purpose of comparison of methods, one tibia with 27 HLs was created and then scaled to 44, 39.2, and 34.5 cm. For the purpose of the Hunt and Hatch method, it is assumed that the mean tibia length for males is 39.2 cm. The outcomes of a different method are presented in Table 3

The second model includes the distal part of juvenile tibia, 26.5 cm long bone of an 11-year-old, 20.5 cm bone of a 7-year-old, and 16.5 cm long bone of a 5-year-old (Fig. 6). Calculations for the proximal part of the tibia are not possible with the Maat or Hummert and van Gerven method; therefore, this part was not included in this model.

Fig. 6
figure 6

Juvenile bones models for the comparison of methods. For the purpose of comparison of methods, three juvenile tibiae with six HLs were created. The Maat and Hummert and van Gerven methods can be used to evaluate the distal part of the tibia only; therefore, the presented model is based solely on this part of the bone. The outcomes of a different method are presented in Table 4

Bland-Altman plot

The outcomes of different methods were compared with the Bland-Altman test (Bland and Altman 1986). It is intended to clearly indicate the differences between the outcomes of two distinct methods. The x-axis presents the mean value of each pair of measurements, the y-axis—the difference between the results of the methods for particular measurements. The more consistent the methods, the more densely distributed the dots and the lower mean difference obtained. The plots were created in the R program (R Core Team 2018) with the BlandAltmanLeh package (Lehnert 2015).

The results are presented in Table 3. Selected Bland-Altman plots are presented in Figs. 7 and 8. Bland-Altman plots for all possible pairs of methods are presented in Supplementary File 2.

Table 3 Results of different methods
Fig. 7
figure 7

Selected Bland-Altman plots comparing the methods used for adult bones. Bland-Altman plots clearly show concordance between the two compared methods. The middle dashed line shows the mean difference between the results of those two methods. The lower and higher dashed lines show ± 1.96 standard deviations. The x-axis presents the mean results calculated for each HL and y-axis—the difference between the first and second of the compared methods for each measurement. All results are given in years. a Results of the Allison method for two bones of different lengths. Varying length significantly affects the outcomes of this method. b The results of the Hummert and van Gerven method for long and short tibiae. The difference is even greater than in the Allison method. c The Hummert and van Gerven method for mean bone compared with the Byers method. Both methods show great concordance, which implies that the Hummer and van Gerven method is accurate for mean bones. The Byers method is also compared with the Allison (d), Clarke (e), and Maat (f) methods

Fig. 8
figure 8

Bland-Altman plots comparing the methods used for juvenile bones. While the Maat and Byers methods give very similar results, the Hummert and van Gerven method produces different outcomes. The cause is rooted in the different source of the comparison table used in the last method

Results and discussion for the adult bone model

The Clarke, Byers, and Maat methods give results independent from the bone length because these three methods are based on ratios rather than absolute measurements. On the contrary, the Allison and Hunt and Hatch methods are very sensitive to different bone lengths. The Allison is based on ratios; however, the assumed 90 mm of initial bone length (regardless of its total length) strongly influences the following calculations in bones of different lengths and, subsequently, the outcomes. The mean difference between the estimated age of HL formation for 44 and 34.5-cm bones is about 0.75 years and ranges from 0 to 2.5 years (Fig. 7a).

Even worse difference between long and short tibiae is obtained with the Hunt and Hatch method (mean difference 3.1 years, range 1–7 year., Fig. 7b). When compared with other methods, age for long tibiae is overestimated and for short—underestimated. It is noteworthy that the outcomes of the Hunt and Hatch method for 39.2-cm tibia give results comparable with the accurate Byers or Maat method (Fig. 7c). It can therefore be assumed that the Hunt and Hatch method may easily be improved by plotting distinct growth curve for each examined bone, based on its length. In such proceedings, the Hunt and Hatch method would become de facto ratio-based. However, without such modification, this method is inaccurate in the estimation of HLs formation in particularly long or short bones.

The Clarke method tends to underestimate the age of HLs formation (besides the earliest ones) when compared with the Byers or Maat method (Fig. 7d), as the thickness of epiphyses was not taken into consideration by Clarke in his calculations.

Byers and Maat methods should be considered the most accurate—they are based on correct premises and give almost identical results (Fig. 7e). However, as the Byers method is faster and allows to perform calculations on a broader range of bones, the authors of the current review strongly recommend choosing this method.

Results and discussion for the juvenile bone model

In this model, the results of the Maat, Hummert and van Gerven, and modified Byers methods were compared (Table 4). While the Maat and modified Byers methods give very similar results (differences ranging from − 0.5 to 0.5 years, Fig. 8b), the Hummer and van Gerven method gives quite different results (range 0–1 years, mean difference 0.4 years, Fig. 8a, c). It does not mean that the Hummert and van Gerven method is unreliable—the cause is rooted in the comparison tables.

Table 4 Results of different methods for juvenile bones

The comparison tables for Maat and Byers are based on the same dataset and therefore the results produce very similar outcomes, Hummert and van Gerven—on the medieval Nubian population. The growth curves for medieval Nubian and twentieth century American population differs—which is reflected on the Bland-Altman plots.

This comparison shows that the assumption of the modified Byers method is correct. Therefore, it can be applied to the proximal part of the tibia, femur, radius, and humerus—which was not possible with the previous methods.

General remarks

In this review, the Byers method is recommended for calculations on adult bones and the modified Byers method—on non-adult bones. However, there are some issues which need to be raised, regardless of the method used.

Every reviewed method is linked with one false although indispensable assumption: “course of bone growth is identical for everyone; it differs only between the sexes.” It is obviously untrue; the beginning of the adolescent growth spurt and its duration varies between individuals (Anderson et al. 1963). The estimation of HLs formed in this period is always biased, regardless of the method used. To the best of the authors’ knowledge, there is no possible way to eliminate this error; however, a study on a large population should overcome the severity of this problem in statistical inference.

The majority of methods are based on some bone growth databases. However, these databases should not be considered universal. Besides the obvious difference in the mean adult bone length between different populations, the course and dynamics of bone growth may also vary significantly; these differences are well depicted by Bland-Altman plots comparing the modified Byers and Hummert and van Gerven methods—which are based on datasets originating from quite different populations.

If this issue is not taken into consideration, it may disturb the results and their interpretation. For example, Piontek et al. (2001) used the Byers method for adult and the Hunt and Hatch method for non-adult bones. They calculated that the mean age at formation of HLs in juvenile bones is 9–11 years of age, and in adult male bones, 11–12 years of age. However, the Hummert and van Gerven method underestimates the age at formation of HLs when compared with the Byers method (as it is presented in the Bland-Altman plot). If they would use the same datasets, most probably they would find similarities rather than differences in this aspect of their study.

Moreover, the creation of new growth curves or comparison tables dedicated to the study population should always be considered, as it was done by Alfonso et al. (2005). However, such approach is rarely used and a vast number of studies are based on the comparison tables provided by authors of distinct methods.

As shown in the “Comparison of all methods” section, most probably every study based on the Allison, McHenry, Clarke, or Hunt and Hatch method is biased and could be reevaluated. The Bland-Altman plots provided in Supplementary File 2 should facilitate the reinterpretation of results in the publications in which the aforementioned methods were used.

Summary, perspectives, and challenges

The interest in Harris lines seems to have slightly increased over the last 18 years, basing on Google Scholar database records (Fig. 9). However, in the last few years, authors have been focusing on the presence and number of HLs rather than the age at their formation. According to publications indexed in the Google Scholar database, there was only one study using the Byers method published in the last 5 years (Krenz-Niedbała 2014) and, similarly, there was only one using the Hummert and van Gerven method (Boucherie et al. 2017).

Fig. 9
figure 9

Articles on HLs indexed by Google Scholar. The graph shows the number of papers published in the years 2000–2018 (by 21.11.2018) containing the phrase “Harris line”

The age at formation of HLs seems to be avoided recently. It is understandable, due to the controversies mentioned in the introduction, such as lack of correlation with other physiological stress indicators (Alfonso et al. 2005), disappearance of HLs resulting from ongoing bone remodeling (Mays 1995), or appearance of HLs unassociated with any physiological stress (Alfonso-Durruty 2011). Moreover, the seemingly difficult and time-consuming calculations may also discourage from estimating the age at formation of HLs.

The present review is aimed at reversing this trend. The authors hope that this study, together with the attached figures and supplementation materials, will greatly simplify the understanding and calculation of the chronology of Harris lines. The authors also encourage to perform the study primarily on juvenile bones; the HLs on them are less likely to be resorbed. Moreover, these HLs are more pronounced and clear, which could increase the intra- and interobserver agreement. New methods enable the estimation of the age at formation of HLs on juvenile femurs and the proximal part of tibiae—and thus increase the number of obtainable results, which will improve the statistical proceedings on archeological study populations.

The controversies connected with HL etiology should be considered as a challenge rather than an obstacle. Many authors who failed to find statistically significant correlation between the HL and other indicators admit that bone remodeling and the following HL disappearance may be responsible for the lack of such correlation (McHenry and Schulz 1976; Mays 1985, 1995, Nowak and Piontek 2002a, b). A similar study conducted on juvenile bones (which are less likely to resorb HL) could give statistically significant results. Moreover, the studies which used the Allison, McHenry, Clarke, or Hunt and Hatch method could also be reevaluated and repeated, as most probably they are biased.

Also, the histological model of formation of HLs could be improved—there is no unanimity whether the HLs are rather “growth arrest” (Scott and Hoppa 2015) or “growth recovery lines”(Sajko et al. 2011). Current biochemical models assume the involvement of GH and IGF-1 (Alfonso-Durruty 2011); however, bone growth and mineralization may be connected also with other proteins and pathways, such as TAZ/Runx2 and PI3K/Akt (Alkagiet and Tziomalos 2017; Zhang et al. 2017), or could involve miRNA (Han et al. 2017) and long non-coding RNA (Xi et al. 2018). Better understanding of the biochemical base of HL formation could improve the interpretation of results involving these lines.

In conclusion, there is much to discover in the field of Harris lines.