An anisotropic interaction model for simulating fingerprints

Evidence suggests that both the interaction of so-called Merkel cells and the epidermal stress distribution play an important role in the formation of fingerprint patterns during pregnancy. To model the formation of fingerprint patterns in a biologically meaningful way these patterns have to become stationary. For the creation of synthetic fingerprints it is also very desirable that rescaling the model parameters leads to rescaled distances between the stationary fingerprint ridges. Based on these observations, as well as the model introduced by Kücken and Champod we propose a new model for the formation of fingerprint patterns during pregnancy. In this anisotropic interaction model the interaction forces not only depend on the distance vector between the cells and the model parameters, but additionally on an underlying tensor field, representing a stress field. This dependence on the tensor field leads to complex, anisotropic patterns. We study the resulting stationary patterns both analytically and numerically. In particular, we show that fingerprint patterns can be modeled as stationary solutions by choosing the underlying tensor field appropriately.


Introduction
Large databases are required for developing, validating and comparing the performance of fingerprint indexing and identification algorithms. The goal of these algorithms is to search and find a fingerprint in a database (or providing the search result that the query fingerprint is not stored in that database). The database sizes for fingerprint identification can vary between several thousand fingerprints e.g. watchlists in border crossing scenarios or hundreds of millions of fingerprints in case of the national biometric ID programme of India.
Extended author information available on the last page of the article Clearly, fingerprint identification is of great importance in forensic science and is increasingly used in biometric applications. Unfortunately, collecting databases of real fingerprints for research purposes is usually very cost-intensive, requires time and effort, and in many countries, it is constrained by laws addressing important aspects such as data protection and privacy. Therefore, it is very desirable to avoid all these disadvantages by simulating large fingerprint databases on a computer.
The creation of synthetic fingerprint images is of great interest to the community of biometric and forensic researchers, as well as practitioners. The SFinGe method (Cappelli et al. 2000) has been proposed to this end by Cappelli et al. (2000). This method can produce fingerprint images which look realistic enough to deceive attendees of a pattern recognition conference, however, systematic differences between real fingerprints and synthetic images by SFinGE regarding the minutiae pattern have been found which allow to distinguish between the two (Gottschlich and Huckemann 2014). Recently, the realistic fingerprint creator (RFC) (Imdahl et al. 2015) has been suggested to overcome the issue of unrealistic minutiae distributions. SFinGe and RFC are both based on Gabor filters (Gottschlich 2012) for image creation. A different approach to fingerprint creation has been introduced by Kücken and Champod (2013). They strive to directly model the process of fingerprint pattern formation as it occurs in nature and their approach is inspired by existing knowledge from biology, anatomy and dermatology. Two commonalities of Gabor filters based and biology-inspired approaches are that both start with random initial conditions and both perform changes in an iterative fashion. Kücken and Champod suggest a model describing the formation of fingerprint patterns over time based on the interaction of certain cells and mechanical stress in the epidermis (Irmak 2010).
In principle, a nature-inspired model nourishes the hope of producing more realistic fingerprints and potentially also to gain insights into the process of natural fingerprint pattern formation. An extensive literature (Champod et al. 2016;Dell and Munger 1986;Irmak 2010;Kücken and Champod 2013;Moore and Munger 1989;Morohunfola et al. 1992;Wertheim 2011) in the biological community suggests that fingerprint patterns are formed due to the interaction of mechanical stress, trophic factors from incoming nerves and interactions between so-called Merkel cells. Merkel cells are epidermal cells that appear in the volar skin at about the 7th week of pregnancy. From that time onward they start to multiply and organise themselves in lines exactly where the primary ridges arise (Kücken and Champod 2013).
The development of fingerprints can be described by three phases (Kücken and Champod 2013). In the first phase, growth forces in the epidermis and shrinkage of volar pad create compressive mechanical stress, modeled by Newell (2004, 2005). The rearrangement of Merkel cells from a random configuration into parallel ridges along the lines of smallest compressive stress forms the second phase. This phase can be regarded as the actual pattern forming process, was first modeled by Kücken and Champod (2013), and is studied in this paper. In the third phase, the primary ridges are induced by the Merkel cells.
The fingerprint development based on the rearrangement of Merkel cells was first modeled by Kücken and Champod (2013). They propose that Merkel cells are the missing link between the stress distribution in the epidermis and the developing pattern due to their mechanosensing ability. For their mathematical description they use an agent-based model to describe the pattern formation process in the second phase of the fingerprint development where the underlying stress field from the first phase Newell 2004, 2005) is considered as an input. Due to the lack of specific information not all details of their model can be confirmed by experimental observations. Hence, they aim to propose a model as simple as possible that captures the essence of the interaction between Merkel cells and stress distribution. For instance, the sensitivity of their model to initial conditions is consistent with the long standing belief that the pattern arrangement is unique and even for identical twins the fingerprints are different. However, the resulting patterns in the model proposed by Kücken and Champod (2013) do not seem to be stationary which is desirable for describing the formation of fingerprints accurately.
Note that a large range of models exist in literature for describing biological pattern formation, including reaction-diffusion models (Kondo and Miura 2010;Turk 1991;Witkin and Kass 1991) and the elastic instability mechanism, see Ball (2009), Koch and Meinhardt (1994), Meinhardt (1982) for good summaries on this topic. A generic partial differential equation, well-known for its pattern-forming behavior, is the Swift-Hohenberg (SH) equation (Swift and Hohenberg 1977). It produces patterns which are locally stripe-like, and upon inspection of simulations (e.g. Stoop et al. 2015), it seems that SH equations can, in principle, produce any patterns occurring in fingerprints, including defects such as triradii and loops in the fingerprint vernacular, and minutiae ends. To the best knowledge of the authors, however, SH equations have never been studied for actual fingerprint simulations. Besides, the well-known existence of an underlying stress field Newell 2004, 2005) is not included in these pattern formation models.
To describe the central phase of the fingerprint development process, i.e. the rearrangement of Merkel cells in the second of the three phases, as accurate as possible the underlying stress field, created in the first phase of the fingerprint development process, has to be considered as an input of our class of models. Motivated by the approach by Kücken and Champod we propose a general class of evolutionary particle models with anisotropic, biology-inspired interaction forces in two space dimensions. In contrast to the Kücken-Champod model, our forces are bio-inspired and we are able to show that fingerprint patterns can be obtained as stationary solutions to our model, an important feature of a biologically meaningful fingerprint development model (Galton 1892;Maltoni et al. 2009;Yoon and Jain 2015). Indeed, the stability of line patterns was the focus of most studies analyzing effects of growth on fingerprints. Sir Francis Galton was among the first to demonstrate scientifically the permanence of the configuration of individual ridges and furrows (Galton 1892). These findings were subsequently confirmed in intensive pediatric research such as Babler (1991).
In our model, we consider a tensor field, modeling the underlying (inhomogeneous) stress field, as one of the inputs of our interaction forces. Besides, the interaction force between two Merkel cells depends on the distance vector between these two cells. We model the coefficient functions of the interaction forces as damped harmonic oscillators, a well-established modeling assumption in cell biology. Besides, this choice reflects the exponential decay of the interaction over larger distances, implying that interactions over very large distances can be neglected, and reinforces an interplay between repulsive and attractive forces as the distance between two cells increases. This choice of the interaction forces is consistent with the general modeling assumption that interaction forces should be short-range repulsive to avoid collisions between cells, and attractive over larger distances to obtain cell accumulations. Note that a similar model is proposed in Burger et al. (2018) and its stationary states are studied both analytically and numerically in the spatially homogeneous case.
Our class of models can be regarded as an biology-inspired adaptation of the Kücken-Champod model (Kücken and Champod 2013) and we describe our modeling assumptions in detail, resulting in a reproducible pattern formation for fingerprints. We show that the resulting stationary patterns depend strongly on the underlying tensor field and the given initial conditions. Perturbations in the initial configuration of the Merkel cells result in perturbed stationary patterns. This situation is analogous to the fingerprints in identical twins who have very similar fingerprints in terms of direction of the ridges and qualitative features of fingerprint lines, but the exact location of ridges and minutiae differs (Jain et al. 2002;Srihari et al. 2006;Tao et al. 2012). Since environmental (within the mother's womb) and genetic conditions are almost identical for twins the differences in defect location are solely due to small perturbations such as the initial configuration of the Merkel cells and the stress field in the epidermis (Kücken and Champod 2013), implying that the fingerprint patterns of underlying identical tensor fields are different but similar. More varied fingerprints can be obtained by changing the underlying tensor field in the model.
In this work, we consider N interacting particles on a domain Ω ⊂ R 2 whose positions equipped with initial data x j (0) = x in j , j = 1, . . . , N . The term F(x j − x k , T (x j )) in (1) denotes the force which a particle at position x k exerts on a particle at position x j . This force depends on an underlying stress tensor field T (x j ) at location x j . The existence of such a tensor field T (x j ) is based on the experimental results in Kim and Holbrook (1995) where an alignment of the particles along the local stress lines is observed. We define the tensor field T (x j ) by the directions of smallest stress at location x j by a unit vector field s = s(x) ∈ R 2 and introduce a corresponding orthonormal vector field l = l(x) ∈ R 2 , representing the directions of largest stress. Then the force is given by for coefficient functions f s and f l .
In previous works on the Kücken-Champod model (Kücken and Champod 2013) and its generalization (Burger et al. 2018) a dynamical system of ordinary differential equations of the form (1) was considered where the force that particle k exerts on particle j is given by i.e. the sum of repulsion and attraction forces, F R and F A , respectively. Here, the attraction force depends on the underlying tensor field T (x j ) at x j , modeling the local stress field. The matrix T (x j ) encodes the direction of the fingerprint lines at x j , defined by and orthonormal vector fields s = s(x), l = l(x) ∈ R 2 . For studying the pattern formation with an underlying spatially homogeneous tensor field T producing straight parallel ridges, e.g.
is considered (Burger et al. 2018). The repulsion and attraction forces in the Kücken-Champod model (Kücken and Champod 2013) and its generalization in Burger et al. (2018) are of the form and respectively. Note that the direction of the attraction force F A and hence also the direction of the total force F are regulated by the parameter χ in the definition of the tensor field T . The parameter χ introduces an anisotropy to the equation leading to complex, anisotropic patterns. For χ = 1 the model (1) with interaction forces of the form (3) for repulsion and attraction force (5) and (6) reduces to a gradient flow and F(d) = −∇W (d) for a radially symmetric interaction potential W . The continuum equation associated with the isotropic particle model (7) is given by x) is the macroscopic velocity field. This continuum model, referred to as the aggregation equation has been studied extensively recently, mainly in terms of its gradient flow structure, the blow-up dynamics for fully attractive potentials and the rich variety of steady states, see Ambrosio et al. (2005), Balagué et al. (2013aBalagué et al. ( , b, 2014, Bernoff and Topaz (2011), Bertozzi et al. (2009Bertozzi et al. ( , 2012, von Brecht and Uminsky (2012), von , Cañizo et al. (2015), Carrillo et al. (2003Carrillo et al. ( , 2006Carrillo et al. ( , 2011Carrillo et al. ( , 2012aCarrillo et al. ( , b, 2016a, Raoul (2010, 2011), Li and Toscani (2004), Raoul (2012), Villani (2003) and the references therein. There has been a trend recently to connect the microscopic and the macroscopic descriptions via kinetic modeling, see for instance (Bellomo and Soler 2012;Carrillo et al. 2010;Ha and Tadmor 2008) for different kinetic models in swarming, Fornasier et al. (2011), Ha and Liu (2009) for the particle to hydrodynamics passage and Karper et al. (2015) for the hydrodynamic limit of a kinetic model. It seems that not many results are currently available in the field of anisotropies. In Evers et al. (2015Evers et al. ( , 2017 anisotropy is modeled by adding weights to the interaction terms. One can show that the model in Evers et al. (2015Evers et al. ( , 2017 is related to our model if a tensor field T is introduced as the velocity direction. Fingerprint simulation results are shown for certain model parameters in Kücken and Champod (2013) where the underlying tensor field is constructed based on fingerprint images using the NBIS package from the National Institute of Standards and Technology. However, Kücken and Champod (2013) is purely descriptive, the choice of parameters is not discussed and the model (1) was not studied mathematically. The model (1) was studied analytically and numerically for the first time in Burger et al. (2018). Here, the authors justify why the particles align along the vector field lines s provided the parameter χ is chosen sufficiently small so that the total force is purely repulsive along s. Besides, the authors investigate the stationary states to the particle model (1) for a spatially homogeneous underlying tensor field where the chosen model parameters are consistent with the work of Kücken and Champod (2013). For the simulation of fingerprints, however, non-homogeneous tensor fields have to be considered, making the analysis of the model significantly more difficult. No analytical results of the long-time behavior of (1) for non-homogeneous tensor fields are currently available. Besides, numerical results for the given model parameters and different non-homogeneous tensor fields are shown over time in Burger et al. (2018) and one can clearly see that the resulting patterns are not stationary. The simulation results for realistic tensor fields for the simulation of fingerprints in Kücken and Champod (2013) seem to be far away from being stationary too. This is illustrated in Figure  9 in Kücken and Champod (2013) where snapshots of the solution are shown for a spatially homogeneous tensor field which should have been parallel lines for steady states. In the biological community, however, it is well-known that fingerprint patterns with their ridge lines and minutiae configuration are determined during pregnancy and remain the same during lifetime (as long as no fingerprint alterations occur). Hence, we are particularly interested in stationary solutions of the system (1).
The goal of this work is to develop an efficient algorithm for creating synthetic fingerprint patterns as stationary solutions of an evolutionary dynamical system of the form (1) as illustrated in Fig. 1d for the underlying tensor field in Fig. 1c.
As a first step we study the existence of stationary solutions to (1) for spatially homogeneous underlying tensor fields and extend these results to certain spatially inhomogeneous tensor fields. Based on these analytical results as well as the stability analysis of line patterns in Carrillo et al. (2018) we can expect stable stationary patterns along the vector field s. Since the solutions to the particle model (1) with the parameters suggested by Kücken and Champod do not seem to be stationary, we investigate Fig. 1 Original fingerprint image and lines of smallest stress s = s(x) for the reconstructed tensor field T = T (x) with an overlying mask of the original fingerprint image in black, as well as stationary solution to the interaction model (1) for interaction forces of the form (2) and randomly uniformly distributed initial data the impact of the interaction forces on the resulting pattern formation numerically. In particular the size of the total attraction force plays a crucial role in the pattern formation. We adjust the model parameters accordingly and simulate fingerprints which seem to be close to being stationary, resulting in an extension of the numerical results in Burger et al. (2018) for inhomogeneous tensor fields. Based on real fingerprint images as in Fig. 1a we determine the underlying tensor field T with lines of smallest stress s by extrapolating the direction field outside of the fingerprint image based on Gottschlich et al. (2009). In Fig. 1b we overlay a mask of the original fingerprint image on the estimated tensor field with direction field s and in Fig. 1c only the direction field s is shown. Besides, we include a novel method for the generation of the underlying tensor fields in our numerical simulations which is based on quadratic differentials as a global model for orientation fields of fingerprints (Huckemann et al. 2008).
In the fingerprint community major features of a fingerprint, called minutiae, are of great interest. Examples include ridge bifurcation, i.e. a single ridge dividing into two ridges. We study how they evolve over time, both heuristically and numerically. Finally, we propose a new bio-inspired model for the creation of synthetic fingerprint patterns which not only allows us to simulate fingerprint patterns as stationary solution of the particle model (1) but also adjust the distances between the fingerprint lines by rescaling the model parameters. This is the first step towards modeling fingerprint patterns with specific features in the future.
Studying the model (1) and in particular its pattern formation result in a better understanding of the fingerprint pattern formation process. Due to the generality of the formulation of the anisotropic interaction model (1) this can be regarded as an important step towards understanding the formation of fingerprints and may be applicable to other anisotropic interactions in nature.
This work is organized as follows. In Sect. 2 the Kücken-Champod model (Kücken and Champod 2013) is introduced and we propose a new bio-inspired modeling approach. Section 3 deals with the existence of steady states to (1) in the form of parallel, equidistant lines for spatially homogeneous tensor fields and its extension to locally spatially homogeneous tensor fields, implying that measurable quantities, such as the almost constant distance between the stationary line patterns, can be predicted with the model. In Sect. 4 we adapt the parameters in the force coefficients (10) and (11) of the Kücken-Champod model in such a way that fingerprint patterns can be obtained as stationary solutions to the particle model (1). Based on these results, we propose the bio-inspired model, described in Sect. 2, to simulate fingerprints with variable distances between the fingerprint lines. For the creation of realistic fingerprints we consider a novel methods for obtaining the underlying tensor field based on quadratic differentials as well as images of real fingerprint data.

Description of the model
In the sequel, we consider particle models of the form (1) where the force F is of the form (2) or (3) where the repulsion and attraction forces are given by (5) and (6), respectively. Note that a model of the form (1) can be rewritten as and can be derived from Newton's second law for particles of identical mass m. Here, λ denotes the coefficient of friction and F j is the total force acting on particle j. Rescaling in time τ = m λ t for small > 0, setting results in the rescaled second order model for small > 0. Starting from (9) the first order model (8) was justified and formally derived in Bodnar and Velazquez (2005) and similar to the rigorous limit from the isotropic second order model (9) to the isotropic first order model (8) as → 0 in Fetecau and Sun (2015) one can proof the rigorous limit of the anisotropic model (9).
Note that setting = 0 in (9) leads to (8), corresponding to instantaneous changes in velocities.

Kücken-Champod model
In the papers (Burger et al. 2018;Kücken and Champod 2013) systems of evolutionary differential equations of the form (1) are considered where the total force, the attraction and the repulsion forces are of the forms (3), (5) and (6), respectively, and the underlying tensor field T is defined as (4). The coefficient functions f R and f A of the repulsion force F R (5) and the attraction force (6) in the Kücken-Champod model are given by and for nonnegative constants α, β, γ , e A and e R , and, again, To be consistent with the work of Kücken and Champod (2013) we assume that the total force (3) exhibits short-range repulsion and long-range attraction along l and we choose the parameters in an initial study as: These parameters are chosen in such a way that the resulting plots of the force coefficients are as close as possible to the ones shown by Kücken and Champod (2013).
Here, the parameter χ ∈ [0, 1] determines the direction of the interaction. For χ = 1 the attraction force between two particles is aligned along their distance vector, while for χ = 0 the attraction between two particles is oriented exactly along the lines of largest stress (Burger et al. 2018). In Fig. 2a the coefficient functions (10) and (11) for the repulsion and attraction forces (5) and (6) in the Kücken-Champod model (1) are plotted for the parameters in (12).
The sums of the coefficients of the forces f R + f A and f R + χ f A for χ = 0.2 are illustrated in Fig. 2b. Note that f R + f A and f R +χ f A are the force coefficients along l and s, respectively. For the choice of parameters in (12) repulsion dominates for short distances along the lines of largest stress to prevent the collision of particles and the force is long-range attractive along the lines of largest stress leading to accumulations of the particles. The absolute value of the attractive force decreases with the distance between particles. Along the lines of smallest stress the particles are always repulsive for χ = 0.2, independent of the distance, though the repulsion force gets weaker for longer distances.  (10) and f A in (11) of repulsion force (5) and attraction force (6), respectively, as well as total force coefficients along the lines of largest and smallest stress for χ = 0.2 (i.e. f A + f R and 0.2 f A + f R , respectively) for parameter values in (12)

Bio-inspired model
We propose a system of ordinary differential equations of the form (1) where the forces are of the form (2). Note that plugging the repulsion and attraction forces (5) and (6) as well as the definition (4) of the tensor field T into the force term (3) results in forces of the form (2). Hence, to generalise the Kücken-Champod model we require for the coefficient functions f s and f l : We model the force coefficients f s and f l in (2) as solutions to a damped harmonic oscillator. Like for the coefficient functions (10), (11) in the Kücken-Champod model we consider exponentially decaying forces describing that short-range interactions between the particles are much stronger than long-range interactions. Besides, the repulsion and attraction forces suggested in the Kücken-Champod model dominate on different regimes. For a more unified modeling approach one may regard this interplay of repulsion and attraction forces as oscillations. This motivates to model the force coefficients f s and f l in (2) as solutions to a damped harmonic oscillator which is also a common modeling approach in cell biology. Hence, we consider the following ansatz functions for the force coefficients f s and f l : for real parameters c, c s , c l , e s 1 , e s 2 , e l 1 , e l 2 , a s , a l . The constants e s 1 , e s 2 , e l 1 , e l 2 control the decay rates of the force coefficients. Since the force coefficients f s and f l both vanish over large distances, this implies that the constants e s 1 , e s 2 , e l 1 , e l 2 are  (14) all negative. Note that c, c s , c l are scaling parameters for the size of the interaction forces. Since f s has to be an exponentially decaying, repulsive force coefficient (i.e. f s ≥ 0) with (possibly) small adaptations, we require that the term c exp(e s 1 |d|) decays exponentially fast and dominates in the definition of f s . Hence, we assume that c is a nonnegative constant with |c| > |c s |. The force coefficient f l is assumed to be short-range repulsive, long-range attractive. Since the cosine function can be regarded as a short-range repulsive, long-range attractive function, this implies that c is nonnegative, consistent with the assumptions before, and |c| > |c l |. Besides, we control the frequency of the oscillations along s and l by positive constants a s , a l , respectively. A possible parameter choice satisfying the above assumptions is given by and we will see that for this parameter choice it is possible to obtain stationary fingerprint patterns and that rescaling of the coefficient functions f s and f l leads to stationary patterns with scaled line distances. The force coefficients f s and f l for the parameters in (14) are shown in Fig. 3. In comparison with the force coefficients F A + f R and 0.2 f A + f R along l and s, respectively, the force f s along s is also purely repulsive, while the force f l is less attractive which is necessary for obtaining stationary patterns as discussed in Sect. 4.2.

General setting
In the sequel, we consider the particle model (1) with force terms of the form F(x j − x k , T (x j )), such as (2) and (3). As in Burger et al. (2018) we consider the domain Ω = T 2 where T 2 is the 2-dimensional unit torus that can be identified with the unit square [0, 1) × [0, 1) ⊂ R 2 with periodic boundary conditions. These boundary conditions have proven to be very useful to simulate interactions on microscopic scales where the simulation domain is large compared to the size of the interacting particles. Besides, periodic boundary conditions are the natural choice in terms of the mathematical analysis and the derivation of the associated macroscopic model. Note that the particles on the domain Ω are separated by a distance of at most 0.5 due to the periodic boundary conditions. Motivated by this we require for j ∈ {1, 2} and all where e j denotes the standard basis for the Euclidean plane. The forces satisfy this assumption if a spherical cutoff radius of length 0.5 is introduced for the forces in (2) or (3), respectively. This assumptions guarantees that the size of the domain is large enough compared to the range of the total force. In particular, non-physical artifacts due to periodic boundary conditions are prevented. A cutoff radius is also very useful to make numerical simulations more efficient. Since our model describes the second phase of the fingerprint development (Kücken and Champod 2013), i.e. the rearrangement of Merkel cells from a random configuration into parallel ridges, we consider randomly uniformly distributed initial data on the torus T 2 in the numerical simulations.

Mathematical analysis of steady states
To use the particle system (1) for the simulation of fingerprints it is of great interest to have a better understanding about the form of the steady states. The steady states are formed by a number of lines which are referred to as ridges. As discussed in Sect. 2 we consider purely repulsive forces along s. In this section, we study the existence of steady states for the particle model (1) for spatially homogeneous and certain inhomogeneous tensor fields T analytically. The stability of these line patterns is further investigated in Carrillo et al. (2018). In particular, the authors show that line patterns for purely repulsive forces along s can only be stable if the patterns are aligned in direction of the vector field s.

Spatially homogeneous tensor field
For spatially homogeneous tensor fields T it is sufficient to restrict ourselves to the tensor field given by s = (0, 1) and l = (1, 0) since stationary solutions to the Kücken-Champod model for any other tensor field can be obtained by coordinate transform (Burger et al. 2018). Further note that steady states are translation invariant, i.e. if x 1 , . . . , x N is a steady state, so is x 1 + z, . . . , x N + z for any z ∈ R 2 . Hence it is sufficient to consider one specific constellation of particles for analysing the steady states of (1). Because of the stability analysis in Carrillo et al. (2018) we restrict ourselves to line patterns along s = (0, 1), i.e. we consider patterns of vertical lines. Note that two-dimensional vertical stripe pattern of width Δ for any Δ > 0 do not satisfy the steady state condition by the analysis in Section 3.3.2 in Burger et al. (2018), i.e. stable line patterns are one-dimensional structures.
Proposition 1 Given |d| ∈ (0, 1] such that n := 1 |d| ∈ N and let N ∈ N be given such that N n ∈ N. Then n parallel equidistant vertical lines of distance |d| of N n uniformly distributed particles along each line are a steady state to the particle model (1) for forces of the form (2) or (3) where the repulsion and attraction forces are of the form (5) and (6), respectively.
Note that the choice of the distance |d| of the parallel vertical lines is consistent with the periodic boundary conditions. Proof Because of the translational invariance of steady states it is sufficient to consider any n equidistant parallel vertical lines of N n particles distributed uniformly along each line. Without loss of generality we assume that the positions of the particles are given bȳ Because of the periodic boundary conditions of the domain as well as the fact that the particles are uniformly distributed along parallel lines, it is sufficient to require that for steady states of the particle model (1). Note that for forces of the form (2) or (3) where the repulsion and attraction forces are of the form (5) and (6), respectively, we have As a first step we show that Note thatx k ∈ {0} × [0, 1] for k = 1, . . . , N n andx N = (0, 0) by the periodic boundary conditions, i.e. we consider all the particles of the vertical line with x 1coordinate x 1 = 0. If N n is odd, then (18) is satisfied by the balance of forces (17). For even N n we have for k = 1, . . . , N 2n − 1. Besides, | = 0.5 and the assumption of the finite range of the forces in (15), implying that (18) is satisfied. If there is an odd number n of parallel equidistant vertical lines, then the condition for steady states (16) is satisfied by (17). For n even, the forces due to particles on the vertical lines at x 1 = k|d| balances the interaction forces due to particles on the vertical lines at x 1 = (n − k)|d| for k = 1, . . . , n 2 − 1 by (17), so it suffices to consider the particles on the vertical line at x 1 = n 2 |d|, i.e. the particles at positionsx k for k = N 2 , . . . , N 2 + N n − 1. Note that . . , N 2 + N n − 1 and the assumption of the finite range of the forces in (15). This implies that the condition for steady states (16) is satisfied. Hence,x 1 , . . . ,x N form a steady state of the microscopic model (1).
Corollary 1 Given d ∈ (0, 1] such that n := 1 d ∈ N and let N ∈ N be given such that N n ∈ N. Then n parallel, but not equidistant, vertical lines of N n uniformly distributed particles along each line are not a steady state to the particle model (1) for forces of the form (2) or (3) where the repulsion and attraction forces are of the form (5) and (6), respectively.
Remark 1 Even though parallel, equidistant lines form a steady state for any distance |d| the line patterns in Proposition 1 are not stable for every |d| ∈ (0, 1]. The maximum distance between parallel equidistant lines is given by the cutoff radius R c of the force coefficient f l or, equivalently, by the distance R c such that f l (|d|) vanishes for all |d| ≥ R c . In particular, a steady state of parallel, equidistant lines of distance R c is also stable under perturbations. This implies that a steady state to (1) of parallel, equidistant vertical lines for a given choice of force coefficients f s and f l can be transformed into a steady state of parallel, equidistant vertical lines with a different line distance by rescaling the force coefficients appropriately.

Non-constant tensor fields
Many non-constant tensor fields can locally be regarded as spatially homogeneous tensor fields. Note that by the assumptions in Sect. 2.3 we consider forces of finite range. In particular, we have local forces for the forces (2) with coefficients (13) and parameters (14) as well as for forces of the form (3) with force coefficients (10), (11) and parameters (12). Applying the results from Proposition 1 and Corollary 1 to a locally spatially homogeneous tensor field implies that the resulting steady states are locally parallel, equidistant line patterns where the distance of the line patterns crucially depends on the range of the interaction forces. In particular, this suggests that the steady states to (1) are given by roughly parallel, equidistant lines whose distance is almost constant. By rescaling the force coefficients the (almost constant) distance between parallel lines can be adapted. This shows that the almost constant distance between (stationary) ridges can be predicted with the model.

Simulation of fingerprint patterns
In this section we investigate how to simulate fingerprint patterns by extending the theoretical and the numerical results in Burger et al. (2018). In particular, we consider more realistic tensor fields for the formation of fingerprint patterns and study the dependence of the parameter values in the Kücken-Champod model on the resulting fingerprints.

Local fields in a fingerprint image
In order to use the particle model (1) to simulate fingerprint patterns a realistic tensor field is needed. It is well known that fingerprints are composed of two key directional features known as cores and deltas. Hence, we consider the construction of the tensor fields for these two features first. In Huckemann et al. (2008) Huckemann et al. propose to use quadratic differentials for generating the global fields in a fingerprint image. The local field is then generated by the singular points of the field. A core is the endpoint of a single line (cp. Fig. 4b) and a delta occurs at the junction of three lines (cp. Fig. 4a).
For simplicity we consider the origin (ζ = 0) as the only singular point, but the idea can be extended to arbitrary singular points ζ ∈ C. As outlined in Huckemann et al. (2008) one can model the field near the origin ζ = 0 by considering the initial value problem for a smooth, positive function φ = φ(r ), r ∈ R, and z 0 ∈ C. For φ = 2 3 the solution to the differential equation is given by Note that the shape of the solution curves does not change for reparametrizations, provided φ > 0. By varying z 0 ∈ C the associated solution curves form a delta at the origin (ζ = 0) as illustrated in Fig. 4a. Hence, we require z dz 2 > 0 for a delta at the origin. Note that z = |z| exp(iarg(z)) where arg(z) denotes the principal argument of the complex number z ∈ C. Further note that dz can be regarded as the direction of the smallest stress at z ∈ C if R 2 is identified with C. As outlined in Sect. 2 the direction of smallest stress is denoted by the unit vector s = s(z) for z ∈ R 2 implying dz = ± exp(−iarg(z)/2). Thus, the lines of smallest stress on a domain Ω ⊂ C can be obtained by evaluating dz for all z ∈ Ω. Note that exp(−iarg(z)/2) and − exp(−iarg(z)/2) result in the same lines of the stress field.  (20) and (22) to the initial value problems (19) and (21), respectively, generating fields of quadratic differentials for a delta and a core Similarly, the initial value problem generates a field with a core at the origin. Up to reparameterization the solution is given by and the solution curves are illustrated for different initial conditions z 0 ∈ C in Fig. 4b. The condition dz 2 z > 0 has to be satisfied for a core at the origin, implying dz = ± exp(iarg(z)/2), and as before ± exp(iarg(z)/2) result in the same lines. Further note that a delta or a core at any ζ ∈ C can be obtained by linear transformation. In Fig. 5 the tensor field for a delta and a core at the singular point (0.5, 0.5) are plotted on the unit square [0, 1] 2 .

Numerical methods
In this section, we describe the general setting for investigating the long-time behavior of solutions to the particle model (1), motivated by Burger et al. (2018).
We consider the particle model (1) where the forces are of the form (2) or (3) and investigate the patterns of the corresponding stationary solutions. As in Burger et al.  (15) to make the simulations more efficient.
To solve the N particle ODE system (1) we apply either the simple explicit Euler scheme or higher order methods such as the Runge-Kutta-Dormand-Prince method, all resulting in very similar simulation results. For the numerical simulations we consider Δt = 0.2 for the size of the time step.

Numerical study of the Kücken-Champod model
Using the tensor fields introduced in Sect. 4.1 we consider the interaction model (1) with forces of the form (3) to simulate fingerprint patterns. Here, the repulsion and attraction forces are of the forms (5) and (6) with force coefficients (10) and (11), respectively, and we consider the parameters in (12) to make the simulations as close as possible to the model suggested by Kücken and Champod (2013). It is well known that fingerprints develop during pregnancy and stay the same afterwards provided no fingerprint alterations occur. In order to simulate biologically meaningful fingerprints we aim to model fingerprint patterns as stationary solution to the particle model (1). Based on the analysis of steady states in Sect. 3 it is possible to obtain stationary patterns consisting of multiple roughly parallel ridges along the lines of smallest stress. However, the force coefficients need to be chosen appropriately so that the resulting patterns are also stable. For the simulations in Fig. 6 we consider the tensor field for the delta constructed in Sect. 4.1 and depicted in Fig. 5a.
One can clearly see in Fig. 6 that the particles are aligned along the lines of smallest stress s = s(x) initially, but the patterns dissolve over time and the simulation results have little similarity with fingerprint patterns over large time intervals. Besides, the patterns are clearly no stable steady states in Fig. 6. Hence, the question arises why the patterns simplify so much over time for non-homogeneous tensor fields in contrast to the stationary patterns arising for spatially homogeneous tensor fields, cf. Burger et al. (2018), and how this can be prohibited. To study the long-time behaviour of the numerical solution, it is desirable to have efficient numerical simulations and of course efficient simulations are also necessary to to simulate fingerprints based on cell interactions in practice. In Sect. 2.3 we introduced a cutoff radius for the forces, given by (15), in order to deal with the periodic boundary conditions. Since the forces in the Kücken-Champod model (1) decrease exponentially, the interaction force between two particles is very small if their distance is sufficiently large. This is also illustrated in Fig. 2a for the parameters in (12). Hence, defining the cutoff radius as 0.1 changes the values of the forces only slightly, but it allows us to compute the numerical solution to the Kücken-Champod model (1) by using cell lists (Allen and Tildesley 1989). The idea of cell lists is to subdivide the simulation domain into cells with edge lengths greater than or equal to the cutoff radius of the interaction forces. All particles are sorted into these cells and only particles in the same or neighbouring cells have to be considered for interactions. This results in significantly faster simulations since we only have to consider those particle pairs with relevant sizes of the interaction forces. Note that the cutoff radius has an impact on the number of lines that occur in the solution as shown in Fig. 7 in comparison to a cutoff radius of 0.5 in Fig. 6. In particular the cutoff radius should not be chosen to small because this prevents the accumulation of particles.
The simulation results for the Kücken-Champod model (1) in Figs. 6 and 7 illustrate that the particles align in roughly parallel lines along the lines of smallest stress initially, but the number of roughly parallel lines decreases as time goes on. In particular, the complex patterns that occur initially are not stationary. We can expect a similar behavior (i.e. initial alignment along the lines of smallest stress of the stress tensor field and subsequent accumulation) of the numerical solution if the parameters in the coefficient functions of the repulsion and attraction force in (10) and (11) are slightly changed provided they are repulsive along the lines of smallest stress, as well as shortrange repulsive and long-range attractive along the lines of largest stress. Denoting the directions of smallest and largest stress by s and l, respectively, the transition of  (1) for different cutoff radii for N = 600 and χ = 0.2 at different times t where the stress field represents a delta the initial pattern of multiple lines to fewer and fewer lines along s suggests that the attraction forces are very strong resulting in an accumulation of the particles. Note that this transition is also observed for the long-time behavior of the numerical solution to the Kücken-Champod model (1) for spatially homogeneous tensor fields in Burger et al. (2018) where lines merge over time until finally a steady state of equidistant parallel lines is reached.
In Fig. 8 we show the numerical solution to (1) for a piecewise spatially homogeneous tensor field, randomly uniformly distributed initial data and N = 600, resulting in stationary line patterns along the lines of smallest stress s = s(x). In particular, this tensor field is not smooth. This suggests that smoothness and periodicity are not necessary to obtain stationary solutions aligned along the lines of smallest stress.
The big impact of the choice of the attraction force along the lines of largest stress can be seen by considering Fig. 9. Here, we assume that the total force is given by for δ ∈ [0, 1] for the spatially homogeneous tensor field T = χ s ⊗ s +l ⊗l with l = (1, 0), s = (0, 1) and χ = 1 instead of the definition of F as the sum of F A and F R in (3), i.e. we vary the size of the attraction force and consider a radially symmetric force F. In Fig. 9 the steady states to the interaction mode (1) are shown for different factors δ of the attraction force F A , where N = 600 and initial data distributed equiangularly on a circle with center (0.5, 0.5) and radius 0.005 is considered. One can see in Fig. 9 that δ = 0.1 results in a stationary solution  (1) for different values of δ (i.e. different sizes of the attraction force F A ) and different axis scalings where χ = 1, N = 600 and radially symmetric initial data (equiangularly distributed on a circle with center (0.5, 0.5)) so that the corresponding stationary solutions are also radially symmetric and radius 0.005

Fig. 10
Total force coefficients δ f A + f R along the lines of largest stress for different values of δ and different scaling spread over the entire domain, while ring patterns arise as δ increases. The intermediate state, occurring for δ = 0.3, is of interest in the sequel, as it is an example of a more complex pattern and in particular not all the particles accumulate on one single ring as for δ = 0.5, δ = 0.7 and δ = 0.9 due to too attractive forces.
The forces considered in Fig. 9 and given by δ f A + f R along the lines of largest stress are plotted in Fig. 10 for different values of δ. As observed in the stationary states in Fig. 10, the force along the lines of largest stress is purely repulsive for δ = 0.1, medium-and long-range attractive for δ ≥ 0.5, as well as medium-range attractive and long-range repulsive for δ = 0.3. In particular, the medium-range attractive forces for δ = 0.3 are significantly smaller than for larger values of δ.

A new model for simulating fingerprints
Based on the analysis of stationary states in Sect. 3 as well as the numerical investigation of the Kücken-Champod model in Sect. 4.3 we propose a new model for simulating the formation of fingerprints based on cell interactions where fingerprints are obtained as stationary states to our model. As a next step we propose a bio-inspired model for the creation of synthetic fingerprint patterns which can not only be used to model the formation of fingerprints as stationary solutions but also allows to adjust the ridge distances of the fingerprint lines.

Stationary patterns
In this section we investigate how fingerprints can be obtained as stationary solutions to the Kücken-Champod model (1) where the coefficients of the repulsive and attractive forces are given by (10) and (11), respectively.

Adaptation of the forces in the Kücken-Champod model
Repulsive forces along the lines of smallest stress are an excellent choice to guarantee that the particles form patterns along the lines of smallest stress. Hence we can consider the repulsive coefficient function 0.2 f A + f R for the force along s with the parameter values in (12) where the coefficient functions f A and f R of the attraction and repulsion force are given by (10) and (11), respectively. Short-range repulsion forces along the lines of largest stress prevent collisions of the particles and medium-range attraction forces are necessary to make the particles form aggregates. However, the long-range forces should not be attractive for modeling complex patterns since strong long-range attraction forces prevent the occurrence of multiple roughly parallel lines as stationary solutions. Motivated by the more complex stationary pattern for δ = 0.3 in Fig. 9 and its desired structure of the forces along the lines of largest stress (short-range repulsive, medium-range attractive, long-range repulsive as depicted in Fig. 10) we consider the coefficient function 0.3 f A + f R along the lines of largest stress for the parameters in (24). Hence, the total force F is given by (3) where the repulsion force F R is defined as (10) and the attraction force F A with coefficient function (11) has the new form where we set T (x j ) = 0.3(l · d)l + χ(s · d)s and we consider the parameter values in (24). In Figs. 11, 12 and 13 the numerical solutions for the repulsive force (5), the attractive force (23) and different realistic tensor fields are illustrated. The tensor fields in Fig. 11 are given by a delta and a core, respectively, introduced in Sect. 4.1, while we consider a combination of deltas and cores for the tensor fields in Figs. 12 and 13. As desired the particles align in roughly parallel lines along the vector field s = s(x) Fig. 11 Tensor fields T = T (x) for delta (a-h) and core (i-p) given by s = s(x) and the numerical solution to the extended Kücken-Champod model (1) with attraction force (23) at different times t for χ = 0.2, N = 600, T = T (x) and randomly uniformly distributed initial data and because of long-range repulsion forces these nice patterns are not destroyed over time. Further note that the numerical solution in Figs. 11, 12 and 13 is shown for very large times so that it can be regarded as stationary. In particular, this implies that the adapted forces can be used to simulate fingerprint pattern and more generally any complex patterns is in principal preserved over time.
After this adaptation of the forces it is desirable to use the original definition of the forces (3) with repulsion and attraction force given by (5) and (6), respectively, instead of an attraction force of the form (23). Along l the attraction force (23) can be regarded as 0.3 f A where f A is the attraction force along l in the original definition of the attraction force F A in (6). Note that the parameter γ in the definition of the attractive force coefficient f A in (11) is a multiplicative constant. Hence, we multiply the original value of γ in (12) by 0.3, resulting in and consider the original definition of the forces in (3), (5) and (6). The forces along the lines of smallest and largest stress are plotted for the parameters in (24) in Fig. 14b. Note that they are of the same form as the adapted forces (3), (5) and (23) for the original parameter values (24), shown in Fig. 14a. Because of the same structure of the forces we can expect similar simulation results. In Fig. 15 the numerical solution is shown for two examples, a delta, as well as a combination of a core and a delta. One can clearly see that the particles align along the lines of smallest stress and the resulting patterns are preserved over time. Similarly, one can obtain any complex pattern as stationary solution to the Kücken-Champod model (1) by adapting the underlying tensor field. In particular, this implies that the Kücken-Champod model (1) with forces defined by (3), (5) and (6) for the parameters in (24) can be used to simulate fingerprint patterns which are in principal preserved over time. The long-time behavior of the numerical solutions to the Kücken-Champod model (1) with model parameters (24) is investigated in Fig. 16 where the numerical solution at large times t is illustrated for the tensor field in Example 5 in Fig. 15. Note that the pattern changes only slightly over large time intervals, demonstrating that these patterns are close to being stationary. This slow convergence to steady states, especially for inhomogeneous underlying tensor fields, can also be seen for other pattern forming systems such as the patterns in the SH equation where the time until the steady state is

Pattern formation based on tensor fields from real fingerprints
In this section, we investigate how to simulate fingerprint patterns based on realistic tensor fields. As proposed in Kücken and Champod (2013) the tensor field is constructed based on real fingerprint data. The tensor field is estimated by a combination of the line sensor method (Gottschlich et al. 2009) and a gradient based method as described in Gottschlich and Schönlieb (2012, Section 2.1).
Given some real fingerprint data the aim is to construct the vector field s = s(x) for all x ∈ Ω as the tangents to the given fingerprint lines. This is based on the idea that the lines of smallest stress are given by s and the solution to the interaction model (1) aligns along s. Let θ = θ(x) denote the angle between the vertical axis and the direction of lines of smallest stress s = s(x) at location x, then it is sufficient to consider the principal arguments θ ∈ [0, π) only. Note that for any x ∈ Ω and any given θ(x) we can reconstruct s(x) as (cos(θ (x), sin(θ (x))) since s(x) are defined to be unit vectors. In Fig. 17 fingerprint data, the estimated arguments θ for constructing the tensor field and the lines of smallest stress s = s(x) of the tensor field are shown. Note that the lines of smallest stress s = s(x) of the tensor field and the fingerprint lines in the real fingerprint image coincide.
Considering the tensor field T = T (x) shown in Fig. 17 the associated numerical solution is plotted for two realizations of uniformly distributed initial data in Fig. 18. One can clearly see that the particles align along the lines of smallest stress s = s(x). Besides, Fig. 18 illustrates that we obtain similar, but not exactly the same patterns for different realisations of random uniformly distributed initial data. This is consistent with the well-known fact that everyone has unique fingerprints and even the fingerprints of twins can be distinguished even if the general patterns may seem to be quite similar at first glance (Champod et al. 2016). To quantify the distance to the steady state we consider the change of the positions x j of the particles in successive time steps, given by In Fig. 19 we show the error τ between successive time steps for the numerical solution in Example 6 in Fig. 18 to the Kücken-Champod model (1). After a sharp initial decrease the total change in positions of the particles is approximately 1.0 · 10 −5 , i.e. the movement of the particles is roughly 1.7 · 10 −8 between time steps.  (24) for the non-homogeneous tensor field T = T (x) in Example 6 in Fig. 18 at different times t and randomly uniformly distributed initial data

Interpretation of the pattern formation
In the simulations for spatially homogeneous tensor fields in Burger et al. (2018) as well as for realistic tensor fields in Figs. 11,12,13,15,16 and 18 one can see bifurcations in the solution pattern for certain time steps. More precisely, there exist points where two roughly parallel lines merge with a third roughly parallel line from the other side. These patterns are in the form of the letter 'Y'. The evolution of one of these bifurcations is shown in Fig. 20 for the underlying tensor field in Example 6 in Fig. 18. Note that all these lines are aligned along the lines of smallest stress s of the tensor field and these bifurcations move towards the two neighboring lines over time. This behavior can be explained by attraction forces along the lines of largest stress over medium range distances, i.e. as soon as the distance between the particles along the lines of largest stress l is small enough they attract each other. In particular, the particles close to the bifurcation on the two neighboring lines are the first ones to 'feel' the attraction force along l and the two roughly parallel lines start merging close to the bifurcation. Hence, the single line on the other side of the bifurcation gets longer over time and the bifurcation moves towards the two parallel lines. While the two roughly parallel lines get shorter over time until they are finally completely merged, resulting in one single line. Since the movement of the particles is mainly along l there is a different particle at the bifurcation at each time step. While the particles on the line in the middle roughly remain at the same position apart from realigning along the lines of smallest stress s. This realignment along s is due to the additional number of particles which are aligned along one single line after the merging, as well as due to the repulsive forces along s spreading the particles to make use of the space along s and to avoid high particles densities after merging.

Motivation for a new model
The results in Sect. 5.1.2 illustrate that it is possible to simulate realistic fingerprints with the Kücken-Champod model (1). As seen in the figures, there is some variability in ridge distances and in view of realistic biometric applications, it is of great interest to control them. Note that the total force F in (3), given by the sum of repulsion and attraction force F R and F A of the form (5) and (6), respectively, can be rewritten as by using the definition of the tensor field T in (4) and the definition of the distance vector d(x j , x k ) = x j − x k ∈ R 2 . The coefficient functions of the repulsion and attraction forces (10) and (11), respectively, are plotted along s and l for the parameters in (24) in Fig. 14b. In particular, this motivates us to consider interaction forces of the form (2). We are interested in rescaling the forces now to vary the distances between the fingerprint lines, i.e. we consider F(ηd(x j , x k ), T (x j )) where η > 0 is the rescaling factor. For η = 1 we recover the same solution patterns as in Sect. 5.1.2, while the distances between the fingerprint lines become larger for η ∈ (0, 1) and smaller for η > 1. Note that the force coefficient f A + f R along l is repulsive over long distances. For η = 1, the case that has been considered so far, this is fine for the given parameters in (24). For η > 1, however, the scaling results in repulsive interaction forces along l for particles with shorter distances between each other. Besides, short-range forces have a stronger impact on the interactions. Hence, these short-range repulsive interaction forces prevent the accumulation of particles along l, resulting in several clusters. Note that the forces along s are purely repulsive so that rescaling by any η does not change the nature of the forces.
In order to prevent this behavior and to obtain an interaction model that can be used for different rescalings, the forces need to be changed slightly so that we have very small attractive forces along l for η = 1. This does not influence the pattern formation for η = 1, but for rescaling by η > 1 we can obtain the desired line patterns with smaller distances between each other. In order to achieve this, we consider a straightforward approach first. We consider two cutoffs c 1 and c 2 and define the adapted force F piece-wise such that for |d| < c 1 the force F is of the form (26) as before while for |d| > c 2 we consider an attraction force tending to zero as d → ∞. To obtain a continuous force we consider a linear interpolation of the force on [c 1 , c 2 ]. Setting we consider the force coefficientsf s andf l for interaction forces of the form (2) where the force coefficients are defined as andf  (27) and (28) respectively, for interaction forces of the form (2) and parameter values (24) Here, we consider the parameter values c 1 = 0.06, c 2 = 0.07 and the parameters in the force coefficients (10), (11) are given by (24). The force coefficient f l along l for |d| > c 2 is obtained by multiplying the original force along l by −1. This is based on the fact that the force coefficient f (d, 1) is repulsive for large distances along l for the parameters in (24). In Fig. 21 the force coefficientsf l andf s in (27) and (28), respectively, are shown. In particular, the piecewise definition off l only has a small influence of the form. In Fig. 22, the stationary solution to the particle model (1) for interaction forces of the form (2), force coefficients (27), (28), parameter values (24), the underlying tensor field T = T (x) in Fig. 17 and different rescaling factors η is shown and one can clearly see that η > 1 leads to smaller ridge distances whereas η < 1 results in larger ridge distances. In particular, the interaction model (1) with interaction forces of the form (2) and force coefficients in (27) and (28) can be used to simulate fingerprints with variable ridge distances. Due to the smaller distances between the fingerprint lines for η = 1.2 this leads to a larger number of fingerprint lines on the given domain. Due to this increased number of lines it is desirable to run simulations with larger numbers of particles. However, particle simulations can only be applied efficiently as long as the total particle number is not too large. In order to solve this remedy one can introduce the density ρ = ρ(t, x) associated with the particle positions and consider the associated macroscopic model In future work, advanced numerical methods for solving the macroscopic model (29) with anisotropic interaction forces could be developed for simulating fingerprint patterns.

A bio-inspired model for simulating stationary fingerprints with variable ridge distances
In this section, we consider interaction forces of the form (2) as before with the aim of simulating fingerprints with variable ridge distances based on a bio-inspired approach.  (2), force coefficients (27), (28), parameter values (24), the realistic tensor field T = T (x) in Fig. 17 and N = 2400 particles initially distributed uniformly at random  (14) as well as piecewise defined coefficientsf l andf s in (27), (28) The coefficient functions f l and f s in (27) and (28), respectively, are defined piecewise and it is desirable to obtain a closed form for the coefficient functions. As before we consider exponentially decaying forces describing that short-range interactions between the particles are much stronger than long-range interactions. Since the forces are repulsive and attractive on different regimes, this interplay between repulsion and attraction forces can be regarded as oscillations. Motivated by this, we model the force coefficients f s and f l in (2) as solutions to a damped harmonic oscillator. Note that harmonic oscillators are a common modeling approach in cell biology and the force coefficients f l , f s are given by (13) and are shown in Fig. 23 for the parameters in (14) in comparison with the piecewise defined force coefficientsf l ,f s for the parameters in (24). Note that the parameters (14) are chosen in such a way that the coefficient functions f l , f s of the harmonic oscillator approximate the piecewise defined coefficient functionsf l ,f s in (27), (28), respectively. In Fig. 24 the stationary patterns to (1) for different rescaling factors η are shown. As expected the larger the value of η the smaller the distances between the fingerprint lines and the more lines occur.  (2), force coefficients (13), parameter values (14), the realistic tensor field T = T (x) in Fig. 17 and N = 2400 particles initially distributed uniformly at random  (13), parameter values (14) and N = 2400 particles initially distributed uniformly at random

Whole fingerprint simulations
In Fig. 25 we construct tensor fields from real fingerprint data based on the methods discussed in Sect. 5.1.2. We consider a whole fingerprint image shown in Fig. 25a and determine the underlying tensor field by estimating the arguments θ = θ(x) for every x ∈ Ω. Since we consider the domain Ω = T 2 we extend the tensor field via extrapolation from the original fingerprint image in Fig. 25a, based on Gottschlich et al. (2009). In Fig. 25b, c the arguments θ = θ(x) are shown and the arguments θ are overlayed by the mask of the original fingerprint in black in Fig. 25b. Since s(x) is a unit vector and hence uniquely determined by its argument θ(x) we reconstruct the lines of smallest stress s(x) as (cos(θ (x), sin(θ (x))) in Fig. 25d, e, and overlay the direction field s by the original fingerprint image in black in Fig. 25d. We run simulations for these realistic tensor fields using our new bio-inspired model (1) with interaction forces of the form (2), force coefficients (13) inspired from harmonic oscillators and parameter values in (14) for randomly uniformly distributed initial data and N = 2400 particles. Note that the patterns are preserved over time.
In conclusion, fingerprints with variable ridge distances can obtained as stationary solutions to our bio-inspired model. We consider harmonic oscillators as force coefficients, a well-established modeling approach in biology. Due to lack of experimental data the exact form of the interaction forces, including the parameter choices, cannot be validated with experiments. For this reason, the parameters are chosen such that certain observations are satisfied and the general model formulation of the model allows to consider a large class of models. As part of future work, the numerical results can be tested for realness. The distinction between real and synthetics could be based on Gottschlich and Huckemann (2014) where histograms of minutiae and ridge frequencies are considered. Another procedure for distinguishing real and synthetic fingerprints is based on the underlying stress field only (Imdahl et al. 2018).