## Abstract

Lichens fix carbon dioxide from the air to build biomass. Crustose and foliose lichens grow as nearly flat, circular disks. Smaller individuals grow slowly, but with small, steady increases in radial growth rate over time. Larger individuals grow more quickly and with a roughly constant radial velocity maintained over the lifetime of the lichen. We translate the coffee drop effect to model lichen growth and demonstrate that growth patterns follow directly from the diffusion of carbon dioxide in the air around a lichen. When a lichen is small, carbon dioxide is fixed across its surface, and the entire thallus contributes to radial growth, but when a lichen is larger carbon dioxide is disproportionately fixed at the edges of an individual, which are the primary drivers of growth. Tests of the model against data suggest it provides an accurate, robust, and universal framework for understanding the growth dynamics of both large and small lichens in nature.

## 1. Introduction

Lichens are symbioses of fungi and photosynthetic algae or bacteria, and are ubiquitous on our planet [1,2], found at the poles, in boreal, temperate and tropical forests, and in deserts and other biomes. They grow on rocks, bark or leaves, soil, and other substrates, weathering the rock, stabilizing soil and providing animals with food, shelter and camouflage. In the environment, lichens are important carbon and nitrogen sinks [3,4], dominating approximately 8% of Earth's terrestrial ecosystems [1].

Lichens grow slowly, and may become very old. Experimental data on the growth rates of crustose and foliose lichens, which are close to circular disks when mature [5], suggest many species grow with similar dynamics. Small lichens grow slowly, but at a steadily increasing rate; growth rates level to a constant as individuals reach a larger size. Despite the variability inherent in ecological data taken from nature [6], this pattern is generally observed and is supported by the field of lichenometry, which uses lichens to date geological events [7], for example rockfalls [8], by using thallus size and measured or extrapolated growth rates to calculate an age of the substrate. Although lichenometry is a popular technique [9], the forces shaping the growth rates of lichens remain unclear [10,11].

Various models have been developed to describe lichen growth, with the most successful emphasizing the fixation and movement of carbon within a thallus [12–14]. The models assume atmospheric carbon is fixed uniformly over the surface of a thallus, and that *internal transport* of carbon to the edge causes radial expansion. Qualitatively, these assumptions reproduce observed patterns of growth: when the lichen is small the entire structure contributes to carbon flux towards the edge, hence growth rates initially increase with the area of the lichen. But above a critical size, internal transport within the lichen cannot keep up with growth; expansion reaches a steady state where only a fixed band close to the outer edge contributes to carbon flux and biomass growth.

Tests of the models' assumptions are difficult and depend on detailed knowledge of carbon fixation, respiration and carbon flux within a thallus. A thorough discussion of what is known and not known about the carbon economy of lichens is provided by [15]. To date, models have estimated the different parameters by fitting predictions to available growth rate data [13,14].

We propose a simpler model of lichen growth, based on a previously overlooked, fundamental fluid mechanical constraint on carbon flux. The model predicts a universal limit to lichen growth as a direct consequence of the diffusion of carbon dioxide in air, with no assumptions about the specific nature of metabolic rates or carbon movement inside a lichen. The central idea is adapted from a widely known study of stains left by evaporating coffee drops [16]. Stains take a specific pattern because of the evaporation of water from the coffee drop, and in our model, absorption of carbon dioxide by the lichen takes the same role in lichen growth as water evaporation does in the coffee drop problem. As carbon is fixed by a thallus, the carbon dioxide in the air close to the lichen surface is depleted. Depletion causes an uneven diffusive flux of carbon dioxide towards the lichen, stronger near the edges versus at the centre of a lichen (figure 1). This mechanism, considered by itself, reproduces the saturation of the growth rate, even if the internal transport of carbon within a lichen is entirely neglected. The model leads to a quantitative prediction for the expansion rate of large lichens, that depends only on the diffusion constant *D* of CO_{2} in air; the density of carbon as CO_{2} in air *ρ*_{air}, the density of carbon in the biomass *ρ*_{lichen}, the height *H* of the lichen, and the fraction of time the lichen is photosynthetically active, *e*.

We present data collected for a population of *Xanthoparmelia* lichens growing in Petersham, MA, USA, and revisit published data of lichens in other genera. Data provide empirical support for our model, moreover, the data of all species collapse onto a single growth curve, suggesting the growth dynamics of small lichens are also governed by universal principles. We extend our model assuming that, even for small lichens, carbon dioxide flux is the limiting growth factor. Height, a rarely studied feature of the growth dynamics of circular lichens, emerges as the parameter critical to reproducing experimental data. Moreover, the model implies newly established lichens will grow as spherical balls, a result consistent with our observations. The model suggests new experiments targeting both the shape and ultimate height of young lichens: if our hypothesis is correct, the height and shape of a developing lichen must follow the specific pattern predicted by our model. Although height and shape are measurable parameters, and appear critical to the growth dynamics of lichens, they have been previously overlooked in the study of lichen growth.

## 2. Diffusion of CO_{2} above the lichen enforces an upper limit for growth speed

Because the air is nearly still within the boundary layer close to the lichen, advection can be neglected and at steady state, the concentration of carbon dioxide, *c*, obeys the Laplace equation
2.1The density of photobionts is greatest close to the surface of most lichens, and at a steady state CO_{2} will be completely depleted very near the surface of a thallus. We thus impose absorption boundary conditions, i.e. that *c*_{lichen} = 0 at the surface of the lichen. Far away from the lichen, at the edge of the boundary layer, *c*(∞) = *c*_{CO2}, where *c*_{CO2} is the average concentration of CO_{2} in free stream air. The condition of a boundary layer of nearly still air above a lichen, with a thickness larger than the other length scales in the problem, is satisfied in most environments. Although lichens may live outside of canopies where wind may be strong, natural terrain is never perfectly flat. Typical values of roughness are larger than the height of a lichen even for level terrain with no vegetation [18], thus the mean wind close to a lichen is likely negligible in most cases. Fluctuations of wind velocity would affect carbon dioxide transport by modifying the diffusion constant.

The mass flux *q* of carbon per unit area into the lichen is *q* = − *M*_{c}*D*∇*c*, where *M*_{c} is the molar mass of carbon. We assume that the flux into the substrate the lichen is growing on, for example a rock, is zero. The Laplace equation with this type of boundary condition is a classic problem in physics, and its solutions are well understood. Originally encountered in electrodynamics, where *c* would be an electric potential and *q* the field strength on the surface of an ideal conductor, the problem is also encountered in fluid dynamics, for example in the evaporation of droplets, a problem nearly identical to the present problem except that the direction of the flux is reversed [16,19]. The pattern of diffusion over an absorber depends on its shape: in the context of lichens, we are interested in two asymptotic limits with respect to the two length scales of the thallus, height *H* and circular radius *R*. A schematic cross section of a typical, growing lichen, taken from an authoritative source [17], illustrates the two limits (figure 1). Large lichens (*R* > *H*) resemble circular disks while small lichens (*R* < *H*) are more rounded. Numerical solutions of equation (2.1) obtained with finite element simulations using COMSOL Multiphysics (see figure 1 caption) demonstrate the two distinctly different patterns of carbon dioxide flux toward a small versus a large lichen (figure 1*b*,*c*).

## 3. The growth rate for large lichens

When the lichen is large, its morphology resembles a flat, circular disk (figure 1). For a disk the flux is non-uniform in the radial direction *r* and very high at the edge
3.1where *ρ*_{air} = *M*_{c}*c*_{CO2}. Equation (3.1) corresponds to the electric field close to a circular conducting disc, see [20] §3.12. Although real lichens are not perfectly smooth, the corrugations and structures on a lichen's surface are irrelevant as long as the lichen remains close to a disc (compare figure 1*c* and inset). The use of equation (3.1) is thus justified and the total flux of carbon toward the lichen causes an equal increase in lichen mass, *m*:
3.2plugging equation (3.1) on the r.h.s. of equation (3.2) and using *m* = *πR*^{2}*Hρ*_{lichen} with constant *H*, the radial growth rate of the lichen is
3.3which is independent of lichen radius *R*. Equation (3.3) provides a direct formula for the expansion rate of maturely shaped lichens in terms of measurable parameters. Since smaller individuals are more rounded and grow slower, we refer to equation (3.3) as the maximum growth rate. With *D* = 16 mm^{2} s^{−1}, an average density of carbon in the lichen *ρ*_{lichen} = 0.63 × 10^{3} kg m^{−3} ([21] and references therein), an average density of carbon in air *ρ*_{air} = *M*_{c}*c*_{CO2} = 2.1 × 10^{−4} kg m^{−3}, and taking for example a reference lichen height of *H* = 4 mm, we can estimate
3.4Lichens will rarely grow at this maximum growth rate because thalli are only photosynthetically active when habitats are the right combination of temperature, moisture, etc. and optimal environments may only occur sporadically. Photosynthesis will also be depressed when thalli are very wet because the diffusion of CO_{2} will slow within the water layer covering the thallus (e.g. [22,23]). To compare our model predictions with real growth speeds, we introduce an additional parameter *e* that represents the mean fraction of time a lichen is photosynthetically active (0 < *e* < 1):
3.5Of the parameters in this model, *D* and *ρ*_{air} are independent of the species and growth conditions and we assume *ρ*_{lichen} is also largely invariant. Hence the model predicts that maximum growth rates of different lichens will vary only when lichens differ in average height *H* or photosynthetic activity *e*.

## 4. Tests of the universal growth saturation

To quantify the saturation of the growth rate, we took advantage of our own direct measures of the growth of 53 individuals (figure 2*a*). Originally, a group of 55 foliose lichens growing on the French tombstone of the North Cemetery, Petersham, MA, USA (42°31′50.20′′ N; 72°11′22.19′′ W) was used to measure growth rates in nature. Target lichens are morphologically uniform and belong to the genus *Xanthoparmelia*. The taxonomy of species within *Xanthoparmelia* is controversial [24,25], and without genetic data we cannot assign a species epithet to the population. Inscribed letters and numbers were used to make a map of the entire population and identify individual thalli from year to year. Each thallus was measured each autumn for 7 years, starting in 2005. A transparent piece of plastic was placed over the thallus and the diameter of the thallus traced with permanent marker. Tracings were digitized and the area, *A*_{n}, calculated from digitized images for each year *n*. Fifty-three new individuals born in 2006, 2007, 2008 and 2009 were added to the survey, and during the survey 33 individuals died. By the end of the survey, in 2011, data had been collected for a total of 75 individuals. Because our model concerns single, isolated, entire lichens, we discarded data for all lichens that fragmented (six individuals), merged with others (nine individuals) or where we noted possible disease or other kinds of damage during the period of observation (seven individuals). We retained a dataset of 53 lichens. For each individual, the growth rate d*R*/d*t* was determined as (*R*_{n+1} − *R*_{n})/1 year, where is the calculated radius of the lichen in year *n*. Raw data for *R* versus *t* as well as d*R*/d*t* versus *R* are shown in figure 2*a* inset and main panel. While data were collected with care and rigor, scatter is an inevitable result of individual variation and environmental heterogeneity (i.e. the relative amounts of shade and sunlight at the top versus the bottom of a tall, columnar tombstone), and is an inherent feature of ecological data. We also identified nine published datasets on the growth rates of seven additional species; all growth rates were similarly measured in nature [11–13,26–29].

All data show saturation of the growth rate, as visualized by the placement of red, flat lines in figure 2. Data also provide evidence for a striking pattern of linear growth rate d*R*/d*t* ∼ *R* among smaller lichens. We use a simple functional form to interpolate between d*R*/d*t* ∼ *R* at small sizes and growth saturation d*R*/d*t* ∼ d*R*/d*t*_{max} at long timescales:
4.1where *r*_{0} marks the transition between linear growth and growth saturation. We then fit equation (4.1) to the data (grey curves in figure 2*a*,*b*), which provides estimates for the two unknown parameters d*R*/d*t*_{max} and *r*_{0}. If carbon dioxide flux is the factor limiting growth over the entire lifetime of a lichen, then we expect the transition to occur when the lichen becomes flat, i.e. when *R* grows larger than *H* (figure 1). Based on this prediction, we identify *r*_{0} ∼ *H* and using equation (3.5) our fitting parameters thus provide estimates for *e* and *H*. To test the consistency of the model, we next searched for published data which independently and experimentally measure the parameters *e* and *H* and found information for three of the seven species considered [11,27,30]. Our calculated values of the fitting parameters *H* and *e* are consistent with their natural range of variation, demonstrating consistency of the model and suggesting the model captures at least the essential elements of growth dynamics (figure 2*c*).

Rescaling time and radius with and we can collapse all data onto the curve (figure 3): 4.2Data collapse is especially tight suggesting that the mechanism underlying growth kinetics is robust and applies to a vast number of crustose and foliose lichen species.

## 5. The growth rate for small lichens

We next explore the consequences of our hypothesis that the entire growth curve, including growth at early stages, is limited by carbon dioxide diffusion in the air. We first calculate the diffusive flux from equation (2.1) toward a small lichen, assuming its shape is approximately a spherical cap with contact angle *θ* and radius *R*. When the contact angle *θ* is large the lichen is round, and when *θ* is small the lichen is flat (see side view sketch in figure 4). This problem has been solved analytically in [19], with the result
5.1where . Equation (5.1) implies that the total flux of carbon dioxide toward the lichen translates to an increase in mass d*m*/d*t* of the spherical cap. Note that equation (5.1) asymptotes to equation (3.2) in the limit for small contact angle, i.e. for a disk. Using *m* = *πρ*_{lichen}*R*^{3}*g*(*θ*) where *g*(*θ*) = (cos^{3}*θ* − 3 cos*θ* + 2)/(3 sin^{3}*θ*) and after simple algebra we obtain the change through time of the contact angle *θ*:
5.2where *F* = *f*(*θ*)/*g*′(*θ*) and *G* = 3*g*(*θ*)/*g*′(*θ*). The negative term on the r.h.s. of equation (5.2) tends to flatten the lichen, whereas the positive term tends to round it up, making it closer to a sphere. Over a long period of time, growth saturates and the lichen's radius increases linearly in time according to equation (3.5): the prefactor of the positive term thus decreases ∼1/*t*^{2} faster than that of the negative term decreasing ∼1/*t*. Eventually the negative term takes over and the lichen flattens, confirming that growth saturation will occur when the lichen is flat. Before flattening, the shape of the lichen undergoes a transition the details of which depend on the exact value of parameters found in equation (5.2), the growth rates of individuals and the initial conditions.

To probe this transition we compute lichen shape over time by solving equation (5.2) using our own experimental data (figure 2*a*). Given *R*(*t*) from data, equation (5.2) predicts the dynamics of the contact angle, and using known values for *ρ*_{lichen}, *ρ*_{air} and *D* (see caption of figure 4), we are left with a single non-dimensional parameter *e*. Figure 4*a* illustrates the shape dynamics of a *Xanthoparmelia* lichen starting from a hemisphere with radius *R*_{0} = 1.5 mm. We used three different values of *e* spanning one order of magnitude, chosen from within the known range of natural variation (figure 2*c*). Large values of *e* correspond to intense photosynthetic activity and consistently yield thicker lichens. Note that the lichen initially rounds up, quickly reaches a maximum contact angle, and then slowly flattens out (figure 4*b*). To verify the generality of this observation, we repeat the analysis with all 10 datasets (figure 2*b*). We use the value of *e* obtained from the fitting procedure described above and start with a hemispherical lichen of radius 1.5 mm as above. The results confirm that the lichens round up quickly before slowly flattening out, as observed for our own dataset. While the maximum contact angle depends on the initial conditions as well as the value of *e*, the transition is qualitatively robust. Interestingly, newly established lichens are in fact very round (A Pringle 2005–2007, personal observation). Figure 5 shows a photograph of a newly established lichen. However, exact and quantitative measurements of lichen shape over time are needed to fully corroborate our prediction.

## 6. Discussion

A simple physical model provides a universal growth limit relevant to tens of thousands of circular, crustose or foliose lichen species. While earlier models focused on internal carbon dynamics [12–14], ours is based only on the flux of carbon dioxide in the air above a lichen. Fundamental constraints on carbon dioxide diffusion are enough to predict observed patterns of growth. Our model clarifies the universal nature of growth limits: lichens cannot manipulate the diffusion of carbon dioxide in air, and differences among species in the efficiency of carbon transport within a thallus are irrelevant to a limit that is solely determined by large scale morphological features.

Other mechanisms will also influence the growth of lichens, and may be sources of error in our model. Growth may be modulated by aspects of algal or bacterial quality, supplies of scarce resources including nutrients scavenged from rainwater or substrates, competition with adjacent individuals, or disease. However, aspects of the model suggest the limiting factor for the growth of circular lichens is, in fact, the diffusive flux of carbon dioxide into a thallus. The growth limit is extremely fundamental and relies on just three basic assumptions: a boundary layer with nearly still air around a lichen, a thallus with a roughly disk-like morphology, and photosynthetic organisms capable of depleting carbon dioxide close to the surface of the thallus. We do not assume that photobionts at the lichen edge deplete carbon dioxide more quickly than photobionts at the centre. When these three assumptions are met, the limit we describe in equation (3.5) has to hold, independently of any other mechanisms or differences among species. Using a typical thallus height of 4 mm, the carbon density of lichens, and the material properties of surrounding air as variables, our prediction of the maximum possible growth speed is about 26 mm yr^{−1}, and is surprisingly close in magnitude to the fastest reported growth speeds, e.g. 13 mm yr^{−1} for a species of *Parmelia* [11,33]. Moreover, the model generates several general and easily verifiable (or falsifiable) predictions: (i) the transition from d*R*/d*t* ∼ *R* and d*R*/d*t* ∼ constant will depend purely on geometry, and will happen after a lichen's radius becomes greater than its height (ii) values for average height *H* and photosynthetic activity *e* can be found by fitting growth data (as measured by changes in radius) to equation (4.1) (iii) for large lichens radial growth is caused primarily by photosynthesis at the lichen's edge. Note that (iii) is not an assumption, but a prediction of the model, as it follows directly from the idea that lichen growth is limited by carbon dioxide uptake. Note also that we do not consider reproduction in our model. We speculate that any photosynthates generated at the centres of larger lichens are disproportionately used to grow reproductive structures, which are often concentrated towards the centre of a thallus and in turn may affect the height and roughness of an individual.

Lichens are not easy to grow in the laboratory and most data are taken in nature. But available field data are inherently noisy; variability likely stems from intraspecific variability, differences in the external environments around individual lichens, and the challenges of accurately recording growth rates of very small or slow growing thalli. The mechanistic principles we offer may also emerge as a useful guide to data fitting; although the general trend is for the growth rate of large lichens to saturate, often several curves can be used to fit the same data [6]. And while tests of the available hypotheses used to explain lichen growth remain challenging, our model suggests various routes forward.

The model's predictions can be tested with experiments, as can hypotheses associated with our analysis of small lichens. Any data that measure radial growth speed would test predictions associated with the growth limit equation (3.5), and independent measures of *H* and *e* would validate the fit of growth rate data to equation (4.1); if the entire three-dimensional shape of the growing lichen was also measured, including the contact angle *θ*, data would verify (or reject) the entirety of our theory. Aspects of geometry (*R*, *H* and *θ*) might be measured directly, while average photosynthetic activity can be measured with tools targeting isotopes or gas exchange. Testing whether the edges of a lichen are highly photosynthetically active, compared to its centre, might involve chlorophyll fluorescence imaging to generate a highly resolved spatial map of photosynthesis. The variations of *H* among individuals of a species can be large, because of the intrinsic variability of environments and individuals, and in part because of reproductive dynamics. Available data for *H* and *e* are snapshots, often from a single individual or single point in time. Rigorous tests of our model would require multiple measures of both *H* and *e* for many individuals in multiple habitats.

The predicted growth pattern depends on a small boundary layer of static air above the lichen surface, and diffusion within the boundary layer. While these assumptions will generally be met in the field, conditions might be modified in the laboratory, if entire lichens were moved indoors. Forcing constant, highly unsteady flow over a lichen surface would make advection dominate over diffusion and in these environments the predicted growth law would disappear.

The model may also explain other features of lichen growth which have already been observed in nature and during experiments. For example, the centre of a thallus often falls out when a lichen is sufficiently large, although the growth rate of the thallus at its edges remains unchanged [11]. Lichens without centres grow as quickly as control lichens with intact centres [34], and shading everything but the small band at the edge of a thallus also has no effect on growth [11]. Our model predicts that as a lichen gets older and larger, the flux of carbon dioxide into the centre of a thallus decreases and eventually becomes negligible: the lack of carbon may cause the centre of a lichen to fall apart or die. Moreover, because most carbon dioxide intake is concentrated at lichen edges, the centre of the lichen will have no direct impact on peripheral growth, explaining why the growth rates of edges remain constant even when centres are missing. The model and available data may also explain why experiments designed to track the movement of carbon within lichens are unsuccessful [11]. Carbon influx is greatest at lichen edges, and this influx drives growth; there may be no movement of carbon within the thallus. In fact carbon dioxide flux in the air above a lichen is at the origin of the growth pattern in our model, and so our thinking is profoundly different from the assumptions of previous models, which focus on the movement of carbon within a thallus and assume transport will shape growth dynamics.

Finally, the mechanism we identify as underpinning lichen growth may be relevant to the growth of other microorganisms with a similar geometry. The canonical example is the growth of bacterial colonies on a Petri dish, epitomized by *Bacillus subtilus* [35]. The geometry of a bacterial colony is similar to the geometry of a circular lichen; a bacterial colony grows out from a centre as a disk over a Petri dish. As with lichens, as the colony grows, there are changes at the centre of the colony [36,37]: in the bacterial colony a bistable switch causes cells to transition from being motile and possessing flagella to expressing extracellular matrix [38]. The first cells that express matrix are within the centre of the colony whereas the motile cells remain on the outside. As the colony develops further there is another transition at the centre of the colony, to sporulation. Intriguingly, the mechanism for the growth of a *Bacillus* colony is similar to the mechanism we describe here, but instead of the diffusion of carbon dioxide, the growth of the bacterial colony is limited by the diffusion of nutrients in the agar [39]. As with the lichen, but from the standpoint of the nutrient diffusion problem in the agar, the colony is a perfect absorber. In other words the mathematical framework that we describe translates directly to bacterial growth. Equation (3.3) predicts the relationship between the colony growth rate and the diffusivity of nutrients, the height of the colony and the relative concentrations of nutrients and bacteria. Using parameters for *Bacillus subtilis* biofilms (see [40] and references therein), diffusivity of glycerol *D* = 5 × 10^{−4} mm^{2} s^{−1}; density of carbon in the medium *ρ*_{agar} = 1.8 × 10^{−4} g cm^{−3}, density of carbon in the cell *ρ*_{cell} = 0.18 g cm^{−3}; biofilm thickness *H* = 0.15 mm and volume fractions from 0.2 to 0.6, we obtain a maximum growth rate –75 μm h^{−1}, comparable with typically observed growth rates of 140 μm h^{−1}. As the colony grows, nutrient concentrations drop in the centre of the colony [40]. Laplace equation (2.1) describes how when the colony grows the nutrient flux drops in the centre of the colony; this nutrient depletion has been demonstrated to cause the bistable switch. A quantitative comparison of these features with bacterial growth dynamics is beyond the scope of our current manuscript but remains an intriguing direction for future research.

Our model highlights the subtle roles physics and fluid dynamics can take to shape the growth and morphologies of organisms. Constraints on carbon dioxide diffusion may also limit the growth of lichens with other body plans, including shrubby, fruticose lichens. But whether constraints do control the growth of lichens with more complex morphologies remains an open, and fascinating, question.

## Data accessibility

Data collected for this study (figure 2*a*) are accessible as supplementary material; all other data (figure 2*b*) have been obtained from the referenced literature.

## Authors' contributions

A.S., J.F. and M.P.B. performed theoretical analysis; A.P. collected experimental data; A.S., J.F., M.P.B. and A.P. designed research, analysed data and wrote paper.

## Competing interests

The authors declare no conflict of interest.

## Funding

This research was begun at the Radcliffe Institute, and has been supported by the Human Frontiers Science Program, and the National Science Foundation Division of Mathematical Sciences DMS 1411694 and DMS 1715477. M.P.B. is an investigator of the Simons Foundation, the French government, through the UCA-JEDI Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-15-IDEX-01 and CNRS PICS 2FORECAST.

## Acknowledgements

Thanks to Jennifer Johnson, David Mauzerall, and the Pringle Laboratory, for reading and commenting on earlier versions of the manuscript. Thanks also to the Harvard Forest, for enabling research on lichens at the North Cemetery.

## Footnotes

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.4108325.

- Received January 23, 2018.
- Accepted May 15, 2018.

- © 2018 The Authors.

Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.