\section{Introduction}
A variety of neurological and musculoskeletal diseases affect human gait quality and manifest in specific stride characteristics. Parkinson's disease (PD), for example, is associated with a reduced stride length, shuffling steps or impaired gait initiation. As reduced gait quality can lead to severe reductions in patient mobility and quality of life \cite{Ellis2011}, it is important to quantify, detect and treat gait impairments as early as possible.

Objective quantification of gait impairment is based on stride-specific characteristics such as stride length or stride time. These parameters are commonly extracted with the help of several electronic measurement systems including computerised pressure mats (\cite{Givon2009,Webster2005}), optical motion-capture systems \cite{Kressig2004} or mobile, sensor-based solutions (\cite{Klucken2013,Mariani2013,Rampp2015,Ferrari2015,Trojaniello2014,Aminian2002,Salarian2013}). While the first two require a laboratory environment and are limited in availability, the latter is mobile and inexpensive. This renders mobile, sensor-based solutions the primary choice for unobtrusive gait analysis systems.

Choosing this modality though introduces a conflict between the abstract variables of measurement and the readout parameters requested by the users and is as such entangled with the physical constraints in wearable sensing: For instrumented medical healthcare applications, one might be able to measure accelerations and angular rates at a patient's foot using state-of-the-art inertial sensors. However, the treating physician is not interested in interpreting acceleration signatures for a given stride but rather wants to monitor variables directly related to the situation as for example stride length or heel-strike angle. The efficient translation of abstract data to context-related knowledge thus is the underlying challenge in all applications of wearable sensors and mobile healthcare technologies.

In the extraction of stride parameters, this challenge is addressed from several perspectives. The majority of methods are based on physical and geometric reasoning to extract spatial gait parameters using double-integration of inertial sensor signals (\cite{Mariani2013,Rampp2015,Rebula2013,Ferrari2015,Trojaniello2014}). The main limitation regarding this type of approach is the dependency on a zero-velocity phase within each stride that is needed to re-initialize the integration process. In clinical practice, however, this assumption is easily violated \cite{Hannink2016}. Other approaches aim at driving bio-mechanical models of the lower extremity with sensor data (\cite{Aminian2002,Salarian2013}) or apply machine learning approaches in order to extract the parameters of interest \cite{Hannink2016}.

The underlying problem of efficient data to knowledge translation is recently being addressed very successfully in the field of image understanding. Here, data from images is identified to belong to a certain object class \cite{He2015}, translated to captions that describe the image content \cite{Vinyals2015} or used to identify persons based on face recognition \cite{Taigman2014}. All these applications represent ground-breaking advances in their respective field in terms of recognition rates. The common underlying methodology that allows these achievements is a branch of machine learning called deep learning.

Due to its success in other domains, deep learning is starting to appear in the context of wearable sensing and computing to extract meaningful information from sensor data (\cite{Hammerla2015,Hammerla2016,Ordonez2016,Hannink2016,Zheng2014}). This particular branch of machine learning is said to have large potential in mobile sensing and computing regarding inference accuracy, robustness or class-scaling which are partly missing from state-of-the-art \cite{Lane2015}. {Applications of deep learning in wearable sensing and computing are, however, largely focussing on activity recognition }(\cite{Hammerla2015,Hammerla2016,Ordonez2016,Zheng2014}). {To the} {authors' best knowledge, this and their prior work }\cite{Hannink2016} {are the first applications regarding other topics in the field}.

In this work, we present a framework  based on deep convolutional neural networks and aim at translating the abstract information provided by wearable sensors to context-related expert features requested by the users. The system is trained on a regression task between sensor data and a set of reference output parameters. {Thereby, it extends the authors' prior work} \cite{Hannink2016} {that only addresses a single output parameter.} A pre-requisite for this is a knowledge base, i.e. a collection of wearable sensor data captured in a controlled environment and annotated with the help of a reference system that can directly measure the expert features of interest.

We apply the proposed framework in the context of mobile gait analysis as illustrated in Fig. \ref{fig:flowchart}. In doing so, we focus specifically on the extraction of biomechanical stride parameters with convolutional neural network regression while potential benefits from deep learning approaches to other parts of the pipeline (e.g. segmentation) might be addressed in future work. A total of eight exemplary and stride-specific characteristics are extracted that are clinically relevant as they define gait quality. To this end, two different modelling approaches are compared: A combined model that uses one network architecture to estimate all expert features of interest and an ensemble approach where one neural network is spawned for each output parameter individually. Both models are trained and evaluated on a publicly available and clinically relevant benchmark dataset.

\begin{figure}[!h]
\definecolor{mygreen}{RGB}{103,204,92}
\definecolor{mygray}{RGB}{133,133,133}
\definecolor{myblue}{RGB}{59,180,242}
\centering
\scalebox{1.05}{
\begin{tikzpicture}[node distance=1.6cm]
\tikzstyle{box} = [rectangle, rounded corners, minimum width=3cm, minimum height=1cm,align=center, draw=black, fill=mygreen,anchor=north]
\tikzstyle{process} = [rectangle, minimum width=3.35cm, minimum height=1cm, draw=black, fill=myblue, align=center,anchor=north]
\tikzstyle{arrow} = [thick,->,>=stealth,color=mygray]

\node (gaitsequence) [box] {Gait sequence\\\includegraphics[width=.1\textwidth]{gaitSequence}};
\node (segmentation) [process,below of = gaitsequence,yshift=.5em] {Segmentation};
\node (HSHS) [process,below of=segmentation,yshift=.5em]{Gait event detection\\ HS-based re-adjustment of stride borders\\Extraction of individual strides};

\coordinate[below of=HSHS] (c);

\node (cnn) [process,left,align=left,left of =c,yshift=-3.5em,xshift=-.1cm] {\underline{With CNN regressor:}\\[.5em]
												{\textendash} Stride length\\
												{\textendash} Stride width\\
												{\textendash} Foot angle\\
												{\textendash} Heel contact time\\
												{\textendash} Toe contact time};
\node (timings) [process,align=left,right of = cnn, right,xshift=.1cm] {\underline{With gait events:}\\[.5em]
												 {\textendash} Stride time\\
												 {\textendash} Swing time\\
												 {\textendash} Stance time\\
												 \null\\
												 \null};

\node (params) [process, fit = (cnn)(timings),inner sep=.2em,fill=none,
				label={[xshift=-5.5em]Parameter estimation}] {};
\node (output) [box,below of = params,anchor=north,yshift=-1em] {Stride-by-stride biomechanical parameters};

\draw [arrow] (gaitsequence) -- (segmentation);
\draw [arrow] (segmentation) -- (HSHS);
\draw [arrow] (HSHS) -- (params);
\draw [arrow] (params) -- (output);
\end{tikzpicture}
}
		\caption{Conceptual flowchart for applying the proposed system to mobile gait analysis: An in-comming gait sequence is segmented and gait events are identified within each segment. In a second step, individual strides are defined from heel-strike (HS) to heel-strike and fed to the biomechanical feature extraction routines. These either compute the timings directly related to the identified gait events or estimate spatio-temporal parameters with a convolutional neural network (CNN).}
		\label{fig:flowchart}
\end{figure}

In summary, our main contributions are: (1) A generalisable method for data-driven and integration-free extraction of spatio-temporal gait characteristics, and (2) Technical validation of the proposed method on a clinically relevant and publicly available dataset.

\section{Methods}
\subsection{Data Collection and Setup}

We use a benchmark dataset collected by Rampp et al. \cite{Rampp2015} that is publicly available at \url{https://www5.cs.fau.de/activitynet/benchmark-datasets/digital-biobank/} and briefly described here.

The inertial sensor platform Shimmer2R \cite{Burns2010} consisting of a 3d-accelerometer (range $\pm 6\,\text{g}$) and a 3d-gyroscope (range $\pm 500\, ^\circ$/s) was used for data collection. It was attached laterally below each ankle joint (Fig. \ref{fig:shoes}). In order to avoid gait changes due to different shoe characteristics \cite{Menant2009}, the same shoe model (adidas Duramo 3) was used by all subjects. Data was captured at $102.4$ Hz at a resolution of $12$ bit. Simultaneously, validation data was acquired with the well established pressure mat GAITRite with a spatial resolution of $\pm 1.27$ cm \cite{Webster2005}.
\FloatBarrier
\begin{figure}[!h]
		\centering
		\includegraphics[width=.6\columnwidth]{shoeSetup.pdf}
		\caption{Placement of the inertial sensor and axes definition.}
		\label{fig:shoes}
\end{figure}

In total, 116 geriatric inpatients were assessed at the Geriatrics Centre Erlangen (Waldkrankenhaus St. Marien, Erlangen, Germany). Written informed consent was obtained prior to the gait assessment in accordance with the ethical committee of the medical faculty at Friedrich-Alexander University Erlangen-N\"urnberg (Re.-No. 4208).

For our study, the annotation on the dataset was extended with additional parameters compared to the annotation reported on by Rampp et al. \cite{Rampp2015}. This was based on positions and timings of the patients' heel and toe as measured by the GAITRite reference system. The spatial parameter set was enlarged to cover not only stride length, but also stride width and change in medio-lateral foot angle. Additionally, heel and toe contact times were added to the list of temporal parameters stride, stance and swing time.  Fig. \ref{fig:gaitrite_parameters} gives an overview on the definitions of temporal and spatial parameters. Stride width was defined as shown in Fig. \ref{fig:gaitrite_parameters} b) and positive values were measured towards the lateral side of the shoe.

\begin{figure*}[!t]
\centering
\includegraphics[width=0.98\textwidth]{GaitRite_Parameters}
\caption{Definition of temporal and spatial gait parameters based on heel and toe positions and timings measured with the GAITRite reference system.}
\label{fig:gaitrite_parameters}
\end{figure*}

Patients performed an extensive geriatric assessment, described by Rampp et al. \cite{Rampp2015} in detail. For the scope of this paper, we focused on the free walking test over the GAITRite mat at comfortable walking speed instrumented with the inertial sensors. After excluding datasets from eight patients due to medical reasons (i.e. patients could not complete the measurement protocol), two due to inertial sensor malfunction and additional seven due to measurement errors with the GAITRite system, 99 patients were left for training and evaluation of the proposed method. Compared to Rampp et al. \cite{Rampp2015}, reference values for heel and toe contact times could not be computed for two patients.

Gait disorders or fall proneness were diagnosed in 54\% of the study population. The other top three diagnoses were heart rhythm disorder (70\%), arterial hypertension (69\%) and coronary artery disease (41\%), which are also associated with gait and balance disorders \cite{Salzman2010}. In summary, this dataset constitutes a clinically relevant study population both in terms of the number of patients and the presence of unpredictable gait alterations.

{TABLE} \ref{tab:datasetStats} gives an overview on the extended set of reference parameters on the dataset and their mean value, standard deviation, as well as their minimal/maximal values.
\begin{table}[!h]

	\caption{Extended set of annotation parameters on the dataset and their mean values, standard deviations, as well as minimal / maximal values.}
	\centering
	\scalebox{\tablescale}{
	\begin{tabular}{llrr}
		\toprule
		Output Parameter & Unit & Mean $\pm$ Std.\hspace{1em} & [\hfill Min,\hfill Max\hfill] \\\toprule
Stride length & cm & 80.63 $\pm$ \hspace{0em} 23.23 & [\hspace{.3em}20.01,\hspace{0em} 129.81] \\
Stride width & cm & -1.44 $\pm$ \hspace{0em} 13.29 & [\hspace{0em}-37.52,\hspace{.5em}  33.03] \\
Foot angle & $^\circ$ &  0.07 $\pm$ \hspace{.5em}  3.49 & [\hspace{0em}-11.93,\hspace{.5em}  15.86] \\ \midrule
Stride time & s &  1.23 $\pm$ \hspace{.5em}  0.19 & [\hspace{.58em} 0.74,\hspace{1em}   2.06] \\
Swing time & s &  0.37 $\pm$ \hspace{.5em}  0.08 & [\hspace{.58em} 0.01,\hspace{1em}   1.05] \\
Stance time & s &  0.85 $\pm$ \hspace{.5em}  0.16 & [\hspace{.58em} 0.48,\hspace{1em}   1.65] \\
Heel contact time & s &  0.64 $\pm$ \hspace{.5em}  0.14 & [\hspace{.58em} 0.16,\hspace{1em}   1.52] \\
Toe contact time & s &  0.69 $\pm$ \hspace{.5em}  0.17 & [\hspace{.58em} 0.25,\hspace{1em}   1.57] \\
		\bottomrule
	\end{tabular}
	}
	\label{tab:datasetStats}

\end{table}

\FloatBarrier
\subsection{Preprocessing}
Before the inertial sensor data is fed to the convolutional neural network, we perform a series of preprocessing steps. {The segmentation step mentioned in the overview } Fig. \ref{fig:flowchart} {is already provided by the dataset in our case. Preprocessing therefore includes} extraction of annotated strides from the continuous recordings, calibration from raw sensor readings to physical units, coordinate system transformations to align sensor axes on left and right feet, normalisation w.r.t. sensor ranges and padding to fixed length of 256 samples per stride to ensure fixed size input to the network.
\begin{figure}[h]
		\centering
		\includegraphics[width=.98\columnwidth]{inputSignal.pdf}
		\caption{Exemplary input signal for a stride defined from HS$\to$HS after preprocessing.}
		\label{fig:inputSignal}
\end{figure}

The system is trained on data segments from heel-strike to heel-strike (HS). This choice of stride definition is beneficial since it does not assume a zero-velocity phase that state-of-the-art double-integration approaches need to re-initialize the integration process. In clinical practice, {this assumption is easily violated as e.g. in the case of spastic gait impairments.} {However, initial ground contact can still be detected for e.g. spastic gait patterns and provides a valid segmentation of the signal into strides. This is the type of scenario that we intend to address with stride segments defined from HS$\to$HS.}

In order to provide this kind of input data, we need to detect HS and toe-off (TO) events in the stride segmentation provided by the dataset and adjust the stride borders accordingly. Detection of HS and TO events within the sensor data is done according to Rampp et al. \cite{Rampp2015}. An exemplary input signal for one stride defined from HS$\to$HS is shown in Fig. \ref{fig:inputSignal}.

Furthermore, we can directly compute three expert-features from the list in {TABLE} \ref{tab:datasetStats} based on the HS and TO events in the datastream. Stride time is defined as the HS$\to$HS distance and given a TO event, we can subdivide each stride into its two phases and compute stance and swing time {(see Fig. }\ref{fig:gaitrite_parameters}{ a) or Rampp et al. }\cite{Rampp2015}). This leaves a set of five output parameters to be estimated by the deep convolutional neural networks.

\FloatBarrier
\subsection{Network Architectures}
\begin{figure*}[!t]
		\centering
		\includegraphics[width=.98\textwidth]{ModelAandB}
		\caption{The two network architectures for data to knowledge translation in the context of mobile gait analysis: Model A consists of three convolutional layers with max-pooling followed by three densely connected layers. Model B, however, spawns smaller networks consisting of two convolutional layers with max-pooling and one densely connected layer for each of the output variables. {Additionally, dimensionalities of important layers are indicated.}}
		\label{fig:modelAB}
\end{figure*}

The network architectures used here are based on three elementary building blocks: Convolutional, max-pooling and densely connected layers.

A convolutional connection between layers is defined by a set of $N_i$ kernels $\vector\psi_1 \dots\, \vector\psi_{N_i}$ of length $L_i$ and biases $b_1 \dots \, b_{N_i}$. The index $i$ thereby represents a label for the layer at hand. Given a multi-channeled input vector $\vector{x}_j$ with $j=1...N_{i-1}$, the activation or output of the convolutional connection is computed as
\begin{align}
	\vector{a}_k = ReLU\left(\sum_j \vector\psi_{k,j} * \vector{x}_j + b_k\right)
	\label{eq:convLayer}
\end{align}
with $k=1\dots N_i$. We used a rectifying linear unit ($ReLU$) as the activation function for this type of connection {as this non-linearity is said to model biological neurons better compared to tangent or sigmoid functions} \cite{Glorot2011}.

Convolutional layers are often followed by max-pooling layers to increase robustness of the extracted features \cite{Boureau2010}. This type of connection downsamples the feature maps obtained by the convolutional connection by taking the maximum in temporal windows of length $r$. The downsampling factor is thus $1/r$.

The third type of connection used here is the densely connected layer. This type of connection is defined by a set of weight vectors $\vector{W}_1 \dots \vector W_{N_i}$ and biases $b_1 \dots b_{N_i}$. Given a single-channel input vector $\vector{x}$, the activation of the densely connected layer is computed by matrix-multiplication as
\begin{align}
	a_j = ReLU\left(\sum_i W_{i,j} \, x_i + b_j\right)
	\label{eq:FullLayer}
\end{align}
with $j = 1\dots N_i$. In case the output of the previous layer is a multi-channelled vector, the single-channel input vector $\vector x$ is constructed by concatenation of the individual channels i.e. flattening. Again, we use a rectifying linear unit for activation.

Finally, a readout layer compresses the last dense layer to the number of output variables for the task at hand. The readout layer is identical to a densely connected layer with the identity instead of the $ReLU$ as an activation function. The number of output variables is then encoded in the number of weight vectors for this layer.

Based on these elementary blocks, two models are built:
\begin{itemize}
	\item Model A: Estimating the complete set of output variables with a combined model (Fig. \ref{fig:modelAB}, left).
	\item Model B: Estimating each output variable individually with an ensemble of networks (Fig. \ref{fig:modelAB}, right).
\end{itemize}
Consequently, the individual network architectures in model B can be less complex compared to model A in order to achieve comparable model complexities.

Regarding the application of data to knowledge translation in the context of mobile gait analysis, we decide for a network architecture built from three convolutional layers with max-pooling followed by three densely connected layers and a readout layer for model A. In the convolutional layers we train $N_1 = 32,\, N_2 = 64 $ and $ N_3 = 128 $ kernels of size $30,15$ and $7$ samples respectively as well as the corresponding number of bias terms. Max-pooling is done in non-overlapping, temporal windows of size $r=2$ samples. Given the sampling frequency, the kernel size corresponds to approximately $0.29$ s on all three layers. The three densely connected layers are trained with $N_4=2048,\, N_5=1024$ and $N_6=512$ weight vectors and bias terms respectively. The readout layer has $N_\text{output} = 5$ nodes for model A.

For model B, the individual network architectures are built from two convolutional layers with $N_1 = 16$ and $N_2 = 32$ filters of size $30$ and $15$ samples respectively and one densely connected layer with $N_3 = 1024$ nodes. The max-pooling layers are identical to model A with a downsampling factor of $\nicefrac{1}{2}$. The readout layer, however, has only $N_\text{output} = 1$ node as each individual architecture is responsible for one of the output parameters. {Fig.} \ref{fig:modelAB} {gives an overview on the network architectures used in both models.}

{The theoretical motivation for this choice is to address the most crucial question in network design of global vs. individual modelling with two representative cases. In model A, we only distinguish between different kinds of output parameters at the last level of the network. The features extracted in this architecture therefore have to be general enough to capture information about all of the output parameters. In model B, however, each output parameter has its own feature extraction path that can be optimized to the parameter at hand.}

\subsection{Training}

Training of neural networks is viewed as an optimization problem regarding an error function (implicitly) depending on the network parameters. This error defines a discrepancy measure between predicted output and a ground truth reference on the training dataset or subsets thereof. Weights and biases on all layers are then changed using back-propagation and with the aim to minimize the error. In practice, however, only random subsets of the training dataset (mini-batches), are shown to the optimizer in one iteration of the training loop to speed up the learning phase (stochastic learning) \cite{Goodfellow2016}.

Because of the different numeric ranges and physical units in the output parameters (see {TABLE} \ref{tab:datasetStats}), the network is trained to estimate normalized and dimensionless output variables $\hat y_i$. Therefore, each reference $y_{i,\text{ref}}$ is scaled to the range $[0,1]$ using the minimum/maximum value attained on the entire training set $\mathcal{S}_\text{train}$:

\begin{align}
	\hat{y}_{i,\text{ref}} = \frac{y_{i,\text{ref}} - \min_{\mathcal{S}_\text{train}} y_{i,\text{ref}}}{\max_{\mathcal{S}_\text{train}} y_{i,\text{ref}} - \min_{\mathcal{S}_\text{train}} y_{i,\text{ref}}}
	\label{eq:scaling}
\end{align}
Predictions on the test set are later rescaled to their physical dimensions using the scaling parameters from the training set.

Given predictions $\hat y_i$ on a mini-batch of size $N_\text{batch}$ for each output variable $i=1\dots N_\text{output}$, we define the error as the sum of the individual root-mean-square errors on the mini-batch:

\begin{align}
	E = \sum_i \text{rmsq}\left(\hat y_i - \hat y_{i,\text{ref}}\right)
	\label{eq:error}
\end{align}

For optimization, we use Adam \cite{Kingma2015}, a state-of-the-art optimization method for stochastic learning. On benchmark datasets, it shows faster convergence than other stochastic optimization routines and we use default settings of $\alpha = 1\text{e}^{-3},\, \beta_1=0.9,\,\beta_2 = 0.999$ and $\epsilon=1\text{e}^{-8}$ (for details see \cite{Kingma2015}).
All weights are initialized by sampling a truncated normal distribution with standard deviation 0.01 and biases are initially set to 0.01. We train for a fixed number of 4000 iterations with a mini-batch size of $N_\text{batch}=100$ strides.

To prevent over-fitting, dropout is used on the densely connected layers. This technique effectively samples a large number of thinned architectures on the hidden layers by randomly dropping nodes during training. With this, over-fitting could be significantly reduced in many use-cases and was superior to weight-regularisation methods \cite{Srivastava2014}. We use fixed dropout probabilities of $p^{(4)} = 0.75$, $p^{(5)} = 0.5$ and $p^{(6)} = 0.0$ for model A. For the individual architectures in model B, a dropout probability of $p^{(3)} = 0.5$ is used. Every connection thus has a 50\% chance of being inactive. During testing, however, the full architectures are used and no connections are dropped.

The networks are implemented and trained using google's TensorFlow library \cite{tensorflow2015-whitepaper}.
\FloatBarrier
\subsection{Evaluation Scheme}

Evaluation of the two modelling approaches is based on a 10-fold cross validation scheme. The stride-specific sensor data from 99 patients on the dataset are sorted into training and test partitions depending on the patient identifier to ensure distinct splits of the dataset.
For each of the two models, we iterate over the complete dataset in this fashion and estimate the output variables on the test set in each fold. The estimates from individual folds are then pooled to arrive at average statistics for each output variable and model. As an evaluation statistic, we use average accuracy $\pm$ precision which correspond to the mean and standard deviation of the signed error distribution. The two models are compared based on this statistic and a Levene test of equal variances between the respective distributions to check whether precisions differ significantly. Because the error distributions for the parameter heel contact time are slightly non-gaussian (checked by visual inspection of q-q plots), a Levene test is preferred over e.g. a Bartlett test that is less robust against non-normality.

In order to assess the learning speed and performance of model A and B, we compute the training error for each of the models over the training iterations for an exemplary and patient-wise 90/10\% train/test split of the dataset.

\FloatBarrier
\section{Results}
\subsection{Training}
Fig. \ref{fig:trainingError} shows the error evaluated on the entire training set over the iterations for an exemplary 90/10\% train/test split of the dataset. The error is evaluated for model A and for each of the submodels that constitute model B. In all cases, the fixed number of 4000 iterations is sufficient to reach a stable regime of the error on the entire training dataset and hence we stop the training. {Furthermore, the adaptation of the two models to the training data is comparable w.r.t. the selected error function as $E_{1\dots5} \approx 0.02$ in model B corresponds to a total/summed error of $E \approx 0.1$ for model A} (Fig \ref{fig:trainingError}).
\begin{figure}[!h]
		\centering
		\includegraphics[width=.98\columnwidth]{trainingStats}
		\caption{Error $E$ evaluated on the entire training set as a function of the training iterations for model A, all submodels that constitute model B and
		 an exemplary 90/10\% split of the dataset.}
		\label{fig:trainingError}
\end{figure}

\begin{table*}[!t]

\caption{{Comparison of model A and B regarding average accuracy and precision reached on unseen test data. To compare average precisions, a Levene-test was performed on the respective error distributions at the 0.01 significance level.}}
	\centering
	\scalebox{\tablescale}{
	\begin{tabular}{lrrrrr}
		\toprule
		& \multicolumn{1}{c}{Stride length}& \multicolumn{1}{c}{Stride width} & \multicolumn{1}{c}{Foot angle} & \multicolumn{1}{c}{Heel contact time} & \multicolumn{1}{c}{Toe contact time} \\\toprule
		Model A& $-0.34 \pm 8.10 $ cm \
          & $0.41 \pm 7.79 $ cm\
          & $-0.05 \pm 3.59 $ $^\circ$\
          & $-0.00 \pm 0.08 $ s\
          & $-0.00 \pm 0.12 $ s \\
        Model B  & $\boldsymbol{-0.15 \pm 6.09} $ \textbf{cm} \
          & $\boldsymbol{-0.09 \pm 4.22 }$ \textbf{cm}\
          & $0.13 \pm 3.78 $ $^\circ$\
          & {$0.00 \pm 0.07$ s}\
          & $0.00 \pm 0.12 $ s \\\midrule
       {Levene-test}  & sign. & sign.& n.s. & {n.s.} & n.s.\\
		\bottomrule
	\end{tabular}
	}
	\label{tab:results-modelcomparison}
\end{table*}

\subsection{Stride Parameter Estimation on Unseen Data}
{TABLE} \ref{tab:results-modelcomparison} lists average accuracy and precision {on the unseen test data} achieved by the two models w.r.t. the pooled estimates from each cross-validation fold. The ensemble approach B that spawns one convolutional neural network for each output variable reaches significantly better precision regarding stride length and width while the corresponding mean accuracies also exceed those achieved by model A. On the remaining three parameters foot angle, {heel} and toe contact time, both models perform similarly. Therefore, we consider the ensemble approach B to be superior in this context.

\begin{figure}[!h]
		\centering
		\includegraphics[width=.92\columnwidth]{BA-ModelB}
		\caption{Bland-Altman plots for each of the output variables estimated by model B {on unseen test data}. Additionally, the mean accuracy (solid line) and ${\pm1.96\sigma}$ bounds (dashed lines) are shown.}
		\label{fig:ba-modelB}
\end{figure}

Detailed results for all output parameters on the dataset and the superior model B can be found in Fig. \ref{fig:ba-modelB} and {TABLE} \ref{tab:results-modelB}. Fig. \ref{fig:ba-modelB} includes Bland-Altman plots for each of the output variables as estimated by model B as well as the achieved mean accuracy and precision. {TABLE} \ref{tab:results-modelB} lists the error statistics for each output parameter as well as state-of-the-art results.

\begin{table*}[!t]
	\caption{{Average estimates as achieved by model B on unseen test data and the reference system GAITRite. The error of measurement is reported as mean accuracy and precision. Additionally, results from state-of-the-art double-integration approaches are listed.}}
	\centering
	\scalebox{\tablescale}{
	\begin{tabular}{llrrrl}
		\toprule
		 &  &   \multicolumn{2}{c}{Average estimates}&  \multicolumn{2}{c}{Mean acc. $\pm$ prec.}
		 \\\cmidrule(l){3-4}\cmidrule(l){5-6}
		Output Parameter&Unit&\multicolumn{1}{c}{Model B}&\multicolumn{1}{c}{GAITRite}&\multicolumn{1}{c}{Model B}&\multicolumn{1}{c}{State-of-the-art}\\\toprule

		Stride length & cm& 80.78 $\pm$ \hspace{0em} 21.82 & 80.63 $\pm$ \hspace{0em} 23.23  &-0.15 $\pm$  6.09 &-0.26 $\pm$ 8.37 $^{\text{\cite{Rampp2015}}}$\\
		&&&&&-0.16 $\pm$ 7.02 $^{\text{\cite{Ferrari2015}},\ast}$\\
		&&&&&\hspace{.33em}0.10 $\pm$ 1.90 $^{\text{\cite{Trojaniello2014}},\ast}$\\

		Stride width & cm& {$-1.34 \pm \hspace{0em}  12.49$} & {$-1.44 \pm \hspace{0em}  13.29$}  &  -0.09 $\pm$  4.22&\null\hspace{2.45em}--$^\text{\tiny(see text)}$  \\
		Foot angle & $^\circ$& -0.06 $\pm$ \hspace{.5em}  2.90 &  0.07 $\pm$ \hspace{.5em}  3.49  &  0.13 $\pm$  3.78 &\hspace{.33em}0.12 $\pm$ 3.59 $^{\text{\cite{Mariani2013}},\ast}$\\\midrule
		Heel contact time & s&  0.64 $\pm$ \hspace{.5em}  0.12 &  0.64 $\pm$ \hspace{.5em}  0.14  &  0.00 $\pm$  0.07 &\null\hspace{2.45em}--\\
		Toe contact time & s&  0.68 $\pm$ \hspace{.5em}  0.14 &  0.69 $\pm$ \hspace{.5em}  0.17  &  0.00 $\pm$  0.12 &\null\hspace{2.45em}--\\
		Stride time & s&  1.23 $\pm$ \hspace{.5em}  0.19 &  1.23 $\pm$ \hspace{.5em}  0.19  & -0.00 $\pm$  0.07 &\hspace{.33em}0.00 $\pm$ 0.07$^{\text{\cite{Rampp2015}}}$\\
		Swing time & s&  0.37 $\pm$ \hspace{.5em}  0.07 &  0.37 $\pm$ \hspace{.5em}  0.08  &  0.00 $\pm$  0.05 &-0.00 $\pm$ 0.05$^{\text{\cite{Rampp2015}}}$ \\
		Stance time & s&  0.86 $\pm$ \hspace{.5em}  0.16 &  0.85 $\pm$ \hspace{.5em}  0.16  & -0.00 $\pm$  0.07 &\hspace{.33em}0.00 $\pm$ 0.07$^{\text{\cite{Rampp2015}}}$ \\
		\bottomrule
		\multicolumn{6}{l}{\hfill\scriptsize Model B: $n=99$ geriatric patients and 1185 individual strides; $^\ast $  different evaluation dataset}\\
	\end{tabular}
	}
	\label{tab:results-modelB}

\end{table*}

\section{Discussion}
We present a method for data to knowledge translation in the context of sensor-based gait analysis to extract a total of 8 spatio-temporal stride characteristics. Within the proposed framework, we compare two different approaches in the case of vectorial knowledge: A) A combined modelling approach estimating the complete set of output parameters and B) An ensemble approach where individual, less complex models are spawned for each output parameter. The resulting model complexities are thereby designed to be of similar magnitude.
With the superior ensemble approach B, the spatial parameters stride length, stride width and foot angle are estimated with mean accuracy and precision of $-0.15 \pm 6.09$ cm, $-0.09 \pm 4.22$ cm and $0.13 \pm 3.78^\circ$ respectively. The temporal stride characteristics stride, swing and stance time are predicted with average precision of $\pm0.07$, $\pm 0.05$ and $\pm 0.07$ s respectively. Additionally, the heel and toe contact times are determined up to $\pm 0.07$ and $\pm 0.12$ s respectively. Thereby, we provide technical validation on a clinically relevant dataset {containing 1185 individual strides from 99 geriatric patients}.

The estimation of stride length is outperforming the double-integration result by Rampp et al. \cite{Rampp2015} by 2.3 cm (27\%, {statistically significant}) in precision. Compared to the deep learning approach presented by Hannink et al. \cite{Hannink2016} which estimates this parameter up to $-0.27\pm5.43$ cm based on HS$\to$HS strides, our results are slightly worse. This might be due to reduced network complexity in the model presented here (Hannink et al. \cite{Hannink2016} use twice as much kernels and biases on each of the two convolutional layers) or the fact that the underlying datasets are not completely identical (101 vs. 99 patients here due to technical reasons). Regarding results on other datasets, Ferrari et al. \cite{Ferrari2015} report a measurement error of $-0.16 \pm 7.02$ cm for the parameter stride length based on double-integration and a dataset of 1314 strides captured from ten elderly PD patients. Trojaniello et al. even reach an average accuracy and precision of $0.1\pm1.9$ cm and thereby almost resolve their reference precision of $\pm 1.3$ cm \cite{Trojaniello2014}. However, this result is only evaluated on a small dataset of 532 strides from 10 elderly PD patients and points out the limits of this comparison: The results achieved have to be seen as a function of the variability across subjects captured in the evaluation dataset. A different evaluation dataset does not necessarily ensure a fair comparison of methods. And it shows the need for a unified evaluation of stride length estimation methods presented here and in the literature (\cite{Rampp2015,Hannink2016,Ferrari2015,Trojaniello2014}) on the same, large cohort study as e.g. the publicly available dataset described by Rampp et al. \cite{Rampp2015}.

Regarding stride width, there is little related work. Horak et al. \cite{Horak2013} even consider it "{difficult to obtain with body-worn sensors}". Nevertheless, Rebula et al. \cite{Rebula2013} report an intraclass correlation coefficient (ICC) of 0.88 between a motion-capture reference and a sensor-based estimation of stride width. Our result corresponds to an ICC of 0.95 and is thus outperforming state-of-the-art double-integration methods w.r.t. stride width.

Mariani et al. \cite{Mariani2013} determine the foot or turning angle up to $0.12 \pm 3.59^\circ$ on data from 10 elderly subjects and a study protocol that included a 180$^\circ$ turn as opposed to our straight walking data. In this respect, our results are comparable.

For the parameters stride, stance and swing time, our results are identical to Rampp et al. \cite{Rampp2015} due to identical methods. Estimation of heel and toe contact times has not been reported in literature to the best of the authors' knowledge. Sabatini et al. \cite{Sabatini2005} propose detection of heel-off and toe-on/foot-flat events by thresholding the angular velocity in the sagittal plane. However, this approach is rather heuristic and the precision regarding these events was not evaluated. Thus, our result constitutes the state-of-the-art regarding sensor-based estimation of heel and toe contact times.

Based on these contact times and HS/TO events that are detectable with state-of-the-art methods, the current work enables detection of heel-off and toe-on events {(Fig.} \ref{fig:gaitrite_parameters}{) }that don't manifest that clearly in the sensor signals. Thereby, each gait cycle can be further sub-divided into loading response, mid-stance phase, terminal stance and pre-swing as defined by Perry et al. \cite{Perry1992} and their dependency on disease state or speed could be evaluated in future work. The latter would extend the work by Hebenstreit et al. \cite{Hebenstreit2015} that was based on motion capture data to a mobile setting.

The entire processing pipeline presented in the current work is based on stride segments defined from heel-strike to heel-strike and therefore independent of the zero-velocity assumption. In everyday clinical practice and the presence of impaired gait, this assumption is easily violated and limits applicability of state-of-the-art double-integration approaches \cite{Hannink2016}. However, there are no theoretical considerations that prohibit application of the proposed method to populations experiencing severe gait disturbances as in the case of spasticity. Thereby, the proposed system is a suitable substitute for the assumption-governed double-integration techniques and could enable mobile gait analysis in these clinically critical cases.

The main limitation of the proposed method is that the resulting data to knowledge translation is as good as the knowledge base. This is because the knowledge base is used to learn the non-linear relationship between the input sensor data and the output parameters and this mapping implicitly depends on the samples collected in the knowledge base. We thus strongly stress the importance of sharing benchmark datasets within the community and create larger, community-maintained knowledge bases.

The implementation of the framework is generalisable and flexible. Information from other wearable sensing modalities (e.g. barometric pressure) can be introduced by adding additional channels to the input signal. Application to other data to knowledge translation problems in this field can be done by exchanging the knowledge base. {TensorFlow }\cite{tensorflow2015-whitepaper} {generically supports model quantisation and lower level arithmetics that are needed for inference with deep convolutional neural networks on mobile devices. Although this is not a necessity in mobile gait analysis, where the emphasis lies on the mobility of the sensing technology, it might be needed in future work.}

Future work includes the application of the proposed framework to other data to knowledge translation problems and thereby the establishment of a generally applicable system. In this respect, especially the end of training needs to be addressed in a data-adaptive manner. {As this work did not include a rigorous exploration of the parameter space (number and dimensionality of kernels, etc.), this part is left for future work.} In the context of mobile gait analysis, individualisation or domain adaptations aspects as well as sequential modelling approaches that account for across-stride context will be investigated. Additionally, the aforementioned benchmark-evaluation of several biomechanical parameter estimation methods for a fair comparison of methods will be covered in future work.

\FloatBarrier
\section{Acknowledgements}
This work was supported by the FAU Emerging Fields Initiative (EFIMoves). The authors would like to thank Samuel Sch\"ulein and Jens Barth for their effort in compiling the benchmark dataset as well as all participants of the study for their contributions.
\bibliographystyle{IEEEtran}


\section{Conclusions}

\label{sec:discussion}
As the only truly wild horse in the world, the Przewalski horses, now mostly living in relatively easily accessible nature reserves, have drawn considerable attention.  Both their collective movements \cite{bourjade2009} and the formation of their harems have attracted interest \cite{ozogany2014}.
However, the unique collective motion displayed by this species, as a large herd consisting of cohesive harems moving together in a coordinated way, has not been modelled to date.

Our model, adapted from a model designed for cells, is  able to qualitatively reproduce the motion of a wild horse herd moving on an open plain, along with formation of groups consisting of one leader and some followers and bachelor groups (group of leaders without followers), with a roughly adequate correspondence of dimension scales. During the analysis of the behaviour of the model we found three interesting phenomena.

First, the herd in our model will at times rotate around its center of mass, while we have not observed the horses to circle, many animals do. This is the direct effect of the fact, that in our model the individuals are unable to stop. Indeed, animals that do rotate along a common axis are usually also unable to stop (e.g.\ flying animals) or is infeasible or dangerous for them to stop. Although some efforts have already been made to model the stopping of a group of animals \cite{kunal2010,Ferdinandy2012}, we suggest further investigations into a model, that would allow for not only the stopping of, but also for the resuming of locomotion.

Second, the translational order parameter is decreased when we introduce followers among the leaders, thus the considered grouping process within the herd effectively reduces locomotion efficiency. In many systems the interactions during motion that give rise to collective motion is for the sake of more efficient locomotion of the group as a whole, but the harem formation within a herd is first and foremost due to reproductive reasons. Thus it is not surprising that the reproductive benefits might outweigh the slight decrease in locomotive efficiency.

Third, the results obtained from our model are not in agreement with the observed group size distribution of the herd that motivated our work (the latter being a lognormal while the former being a normal distribution). Our simple model operates solely with interactions based on spatial distances, while group-forming processes in real societies have many complex attributes, thus deviations from the exact features of the empirical population is expected. On the other hand, the collective motion in many species can be described by purely distance-based interactions, making the exact nature of these deviations non-trivial. Consequently, we propose further investigations of collectively moving systems to find the properties that allow for the spatial formulation of interactions within the system. It can be supposed that in systems where individuals are interchangeable (in the meaning that individual recognition during the motion is not feasible), like a group of cells, ants or a flock of starlings, considering only distance-based interactions is enough to reproduce the observed collective motion pattern, but in animals living in structured social systems (and maintaining an individual recognition), like horses, social factors are much more important during interactions than actual distances, thus interchangeability might be one such property.

\section{Acknowledgements}
We are grateful to the Directorate of Hortobágy National Park for providing us with their data as well as authorizing our observations. This work was supported by the Hungarian Academy of Sciences (grant numbers MTA 01 031 and 2011TKI552).

\section{Introduction and motivation}
\label{sec:introduction}

Living in social structures with multiple levels of hierarchy is widespread in the animal kingdom \cite{Grueter2012a,Grueter2012}. Examples range across several taxa, beginning with humans and primates \cite{Abegglen1984,Kummer1968}, through elephants \cite{Wittemyer2005}, to whales \cite{Baird2000,Whitehead2012} and equids \cite{Feh2001,Rubenstein2004}. There are numerous examples of subgroups  forming around a single individual. For example groups may emerge around a matriarch from her descendants, like in african elephants \cite{Wittemyer2005}, sperm whales \cite{Whitehead2012}, and
killer whales \cite{Baird2000}. Alternatively a reproductive unit may form around a breeding male with several breeding females and their young as in Przewalski
horses \cite{Boyd1994} and plains zebras \cite{Rubenstein2004}. These breeding units can sometimes also include non-breeding males as well, like in hamadryas baboons \cite{Kummer1968} or geladas \cite{Dunbar1975}.

Our aim in the current study was to examine the way in which such a two-level hierarchy may spontaneously emerge in a group and what implications that hierarchy might have for the collective motion of the group. Our motivation and empirical basis was the collective motion of a Przewalski horse herd in light of group formation within the herd, aided by observations made in \cite{ozogany2014} at the Hortobágy National Park in Hungary. As mentioned before, the Przewalski horse herd is split into harems, organized around a breeding male, with several breeding females and their young offspring. So-called bachelor groups, which consist of males that do not have there own harem are also present \cite{hartmann2012}. It should be noted, that although zebra harems form herds in the wild and have a very similar social structure to the Przewalski horse, the Przewalski herd at Hortobágy is only semi-wild as it lives in a bounded environment, which may force them into a herd. Although this has not been studied thoroughly, park officials reported, that the initial population did not form a herd, which only appeared after the growth of population density.

Both the collective motion of several different species of animals \cite{Vicsek2012}, and the emergence of hierarchy within the social system of the Przewalski horse \cite{ozogany2014} have already been modelled. Conversely, the collective motion of animals that are hierarchically organized into subgroups within a larger group have not been modelled. Thus, we aimed to construct a model of group formation and collective motion of a herd composed of sub-groups as a self-propelled particle (SPP) model in two dimensions, where we identified leaders forming harems, and followers making up these harems.

Leadership is a common concept invoked to explain coordinated group movements. In ungulates this is often attributed to a single individual. A recent study \cite{bourjade2015} raises interesting questions about the validity of such a concept, based on observation of two groups of 12 and 6 Przewalski horses. Using different definitions of leadership (moving first, moving in front, or eliciting joining to movement), no individuals that could be consistently classified as a leader were identified. Some limitations to that study are that several types of movements were not measured. In addition movements in the breeding season were also not measured; this was deemed problematic because in the breeding season the stallions directly elicit movements of their harems away from other stallions. Also, due to methodological reasons, only short periods of the day were observed. It has been shown in some cases that in the same group different type of leadership hierarchies might arise in different contexts \cite{Watts150518,nagy2013context}, and there are examples in nature where certain individuals in animal groups consistently act as leaders, for example in zebras and dolphins \cite{larissa2009,fischhoff2007}. These imply that although attributing leadership to a single individual might not be applicable in all circumstances, it does have explanatory power in a wide range of scenarios. As such it stands to reason that conceptualizing the division between leaders and followers dichotomously helps simplify modelling at a minor cost. Simplifying modelling is helpful in the initial understanding of the type of collective movement we analysed in the present article. Nonetheless this recent study has implications that warrant further field studies of leadership in animals. It would be particularly illuminating to have detailed and continuous data on the movement of large groups of Przewalski horses, which is not easy task, since we do not know of any herd, where attaching measurement devices to the animals is allowed by officials. A substitute for real wild horses could be domestic and feral horses, on which studies have also been carried out \cite{briard2015,krueger2014}. Interestingly, these studies concentrate on movement initiation and not the collective motion itself, showing a somewhat different point-of-view in physicists and biologists.

Herein we consider an earlier SPP model of collective cell movement \cite{szabo2006} and extend it with a two level hierarchy by introducing two distinct types of particles (i.e.\ leaders  and followers) while simultaneously attempting to limit the increase in the number of parameters. In contrast with \cite{ozogany2014} the group formation is not driven by the environment of the herd, but by interactions dynamically evolving during the collective motion of the individuals. While formulating the model, we concentrated on mimicking the movements of Przewalski horses. While this specificity adds some complexity to our model, relative to what is usual in statistical physics, it is mostly related to nuances in movement and does not play a major factor in group formation, which was our main focus.

Our study could have potential implications for understanding how and why group formation occurs in nature, how group formation affects the system in which it is happening and the rules governing collective motion in a two-level system. Inferring the universalities and the particulars of the different kind of mechanisms, could potentially be used to artificially control both living and human-made systems, such as domestically kept horse herds or flocks of drones.

\section*{References}

\section{Model}

\label{sec:model}

The model is based on \cite{szabo2006} in which a model was developed to depict the collective motion of cells. We modified this model to accommodate two types of SPP-s (leaders and followers), asymmetric interactions and group formation rules. While extending this model we aimed at minimizing the number of extra parameters. Compared with the usual SPP models the model of \cite{szabo2006} gives smoother results due to intrinsic relaxation times.  We choose parameters that allow the development of motion that resembles the movements of a herd made of harems as close as possible within the framework of the model. We provide a graphical overview of the model in Figure \ref{fig:modeldescription} and an introduction here.

The movements of the horses in the model are confined to a square area, large compared to the size of the herd, representing the herding area available to them (Figure \ref{fig:modeldescription} boundary). Periodic boundary conditions were not considered, first, because it is not realistic, and second, because it does not make sense in a co-moving herd to conceptualize that the front may interact with the rear. Also, we introduce a tendency for horses that stray too far from the herd to head back while still going in the general direction of the herd's (Figure \ref{fig:modeldescription} a)).

All horses may follow all other horses, but the strength of the interaction depends on the types and orientations of the SPP-s in question. Given, that it is plausible that leaders must also pay attention to followers, they will follow followers too, but to a much smaller extent than the other way around. Although the interactions taking place are based on metric distances, we introduce a directedness, meaning that a horse will follow the ones in front of it more than the ones behind it. Several types of interaction modes have been suggested in modelling collective motion. Early models used a simple metric distance, e.g., interacting with anybody nearer than a given distance \cite{vicsek1995}. Later topological distances were introduced, e.g., interacting with a fixed number of nearest neighbours \cite{ballerini2008}. Recently it has been proposed that the most biologically correct interaction ranges should be based on visual perception \cite{strandburg2013}. In our case, vision plays little part as equine vision is near 360° \cite{brooks1999} and neither the distances within the herd nor the density of the herd imply that occlusion would have a major effect on interactions. As such, the effect of following the ones in front, rather than the ones behind is related more to the logic of not turning around if there are others heading in the same direction as oneself.

Leaders who acquire followers (i.e., a harem), will stay farther away from other leaders than if they were without followers (Figure \ref{fig:modeldescription} b) and c)). Harems are established based on spatial distance, but followers will gradually belong more and more to the leader they follow, making it easier for them to stay close, because of the stronger and slightly longer distance interactions with their leader than with another leader (Figure \ref{fig:modeldescription} d)).

Our model starts from randomized initial positions and velocities, without followers being assigned to any leader, thus all followers find groups and leaders at the same time. Our model forgoes the introduction of complex social rules by using only spatial interactions as described above and not taking into consideration that in reality, a new horse would be introduced to a herd already split into harems. On the other hand, taking the latter into consideration would not allow for the study of emergent group formation.

\begin{figure}[ht!]

\centering
        \includegraphics[scale=0.5]{figure1.pdf}

\caption{Graphical overview of the model depicting a small herd inside the boundary with  various parts of $\bm{F}_\text{int,r}(\bm{r}_i,\bm{r}_j)$ and $\bm{F}_\text{com}$ shown. Radii are drawn to scale (cf. Table \ref{tab:parameters} for actual values), and the herd is magnified from within the boundary to show the forces. Solid arrows depict direction of forces, dashed arrows depict actual velocities. The following details are included: a) a horse farther from the center of mass than the given boundary (large green circle centred on the center of mass) will move towards the herd but also in the direction the herd is going, b) \& c) leaders without groups can go closer to each other than to a leader with a group, while followers can go even closer to a leader, d) the attraction radius of the follower-leader interactions is generally smaller than that of the leader-leader interactions, but it is increased when interacting with the leader of the follower's group.
}%
\label{fig:modeldescription}

\end{figure}
\subsection{Formal model description}
\label{sec:formal model}
We have $N_\text{L}$ number of leaders and $N_\text{F}$ number of followers (the list of parameters can be found in Table \ref{tab:parameters}). The 2-dimensional motion of the horse $i \in \lbrace 1$,$ N = N_\text{L} + N_\text{F} \rbrace$ is described by the overdamped dynamics
\begin{equation}
\diff{\bm{r}_i(t)}{t} = v_i^0\bm{n}_i(\theta_i) +
\sum\limits_{\substack{j=1\\j\neq i}}^N\bm{F}_{\text{int}}(r_{ij},\varphi_{ij})
+ \bm{F}_\text{com}(\bar{\bm{r}}-\bm{r}_i,\bar{\bm{v}}) +
\bm{F}_\text{wall}(\bm{r}_i,\bm{v}_i) + \boldsymbol\xi
\end{equation}

where $t$ is time, $\bm{r}_i$ is the position of and $\bm{v}_i$ is the velocity of horse $i$, $v^0$ is a preferred speed which differs for leaders and followers, $\bm{n}_i$ is a unit vector characterized by the angle $\theta_i$,  $\bm{F}_{\text{int}}$ is a pairwise interaction with $r_{ij} = |\bm{r}_i - \bm{r}_j|$ and $\varphi_{ij}$ being the angle between $\bm{r}_i - \bm{r}_j$ and $\bm{v}_i$, $\bm{F}_\text{com}$ is a global force dependant on the position ($\bar{\bm{r}}$) and the velocity ($\bar{\bm{v}}$) of the center of mass of the herd, $\bm{F}_\text{wall}$ is the force acting at the boundaries and $\boldsymbol\xi$ is a vector whose components are delta-correlated white noise terms with zero mean.

The direction of the self-propelling velocity $\bm{n}_i(t)$, described by the angle $\theta_i(t)$, attempts to relax to $\bm{v}_i(t) = \operatorname{d}\!\bm{r}_i(t)/ \operatorname{d}\!t$ with a relaxation time $\tau_i$:

\begin{equation}
\frac{\operatorname{d}\!\theta_i(t)}{\operatorname{d}\!t} = \frac{1}{\tau_i}\operatorname{arcsin}\!\left[\left( \bm{n}_i(t) \times \frac{\bm{v}_i(t)}{|\bm{v}_i(t)|}\right)\cdot \bm{e}_z\right],
\end{equation}
where $\bm{e}_z$ is a unit vector orthogonal to the plane of motion, and $\tau_i$ differs for leaders and followers. This relaxation provides smooth transitions of the $\bm{n}_i(t)$ desired velocities. The value of $\tau$ was chosen larger for leaders than followers, implying that leaders are harder to ''convince'' than followers to change directions, but our results are not sensitive to changes in $\tau$.

The $\bm{F}_{\text{int}}(r_{ij},\varphi_{ij})$ force that carries the direct interaction between the horses can be split into the product of a spatial part ($\bm{F}_{\text{int,r}}(r_{ij})$), and a coefficient part ($F_{{\text{int,}\varphi}}(\varphi_{ij}$), the latter being dependent on the angle of the direction of horse $j$ from horse $i$ and the direction of the velocity of horse $i$. The spatial part consists of a pair-wise, asymmetrical force, the direction of which lies on the line passing through the center of masses of the interacting horses and the magnitude of which is the function of the distance $r_{ij}$ between the horses \cite{szabo2006}. The actual form of the force depends on the type of horses involved:

\begin{equation}
\bm{F}_\text{int,r}(r_{ij}) =
\begin{cases}
\bm{F}_\text{LL}(r_{ij}), & \mbox{if } i \mbox{ and } j \mbox{ are both leaders,}  \\
\bm{F}_\text{FL}(r_{ij}), & \mbox{if } i \mbox{ is a follower and } j \mbox{ is leader,}\\
\bm{F}_\text{LF}(r_{ij}), &\mbox{if } i \mbox{ is a leader and } j \mbox{ is a follower,} \\
\bm{F}_\text{FF}(r_{ij}),& \mbox{if } i \mbox{ and } j \mbox{ are both followers.}
\end{cases}
\end{equation}

For all four cases there are two radii defined, $R^\text{AT}$ which is the range of attraction, and a smaller radii $R^\text{EX}$, which is the range of repulsion, and also a distance $L$, which defines a distance inside $R^\text{AT}$ but outside of $R^\text{EX}$, splitting the force into four parts depending on distance, namely a repulsive, an attractive and two non-interacting regimes, with different coefficients for all four types of interaction in both the interacting regimes ($F^\text{AT}$ for the attractive and $F^\text{EX}$ for the repulsive), thus having 8 radii with 8 coefficients and 4 distances. On the example of $\bm{F}_\text{LF}(r_{ij})$ the equations look like this (leader-leader and follower-leader interactions are slightly different):

\begin{equation}
\label{eq:sampleforce}
\bm{F}_\text{LF}(\bm{r}_i,\bm{r}_j) = \bm{e}_{ij} \times
\begin{cases}
F^\text{EX}_\text{LF} \frac{r_{ij}-R^\text{EX}_\text{LF}}{R^\text{EX}_\text{LF}}, & r_{ij} < R^\text{EX}_\text{LF},\\
0, & R^\text{EX}_\text{LF} < r_{ij} < R^\text{EX}_\text{LF} + L_\text{LF}, \\

F^\text{AT}_\text{LF} \frac{r_{ij}-R^\text{EX}_\text{LF}}{R^\text{AT}_\text{LF} - R^\text{EX}_\text{LF}-L_\text{LF}}, & R^{EX}_\text{LF} + L_\text{LF}\le r_{ij} \le R^\text{AT}_\text{LF}, \\

0, & R^\text{AT}_\text{LF} < r_{ij},

\end{cases}
\end{equation}
where $\bm{e}_{ij} = (\bm{r}_i -\bm{r}_j)/r_{ij}$. The non-interacting part between $R^\text{EX}$ and $R^\text{AT}$ was chosen to be very small its only function being is to remove some ''vibrations'' that arise at such low densities, when a horse is on the edge of the attractive and repulsive regimes. The form of the force is one of the simplest ways to define gradually growing forces based on distances and the values of the specific parameters were chosen to imitate that leaders with harem wish to protect their followers from other leaders, while bachelor leaders themselves can create groups.

In the cases of leader-leader ($\bm{F}_\text{LL}$) and follower-leader ($\bm{F}_\text{FL}$) interaction this picture is slightly changed due to the formation of groups. Followers will develop a certain amount of affinity to leaders who are close by, that increases in strength when they are close to the leader and decreases when they are farther away from the leader. Each follower keeps track of time spent near each leader with the quantities $D_{ij} \in [0,\infty]$, which follow the simple dynamics
\begin{equation}
\label{eq:maledist}
\frac{\operatorname{d}\!D_{ij}}{\operatorname{d}\!t} =
\begin{cases}
+1, & r_{ij} \leq R^\text{AT}_\text{LF}, \\
-1, & r_{ij} > R^\text{AT}_\text{LF} \text{ and } D_{ij} > 0, \\
\quad 0, & r_{ij} > R^\text{AT}_\text{LF} \text{ and } D_{ij} \leq 0.
\end{cases}
\end{equation}
This is then translated into an affinity
\begin{equation}
\label{eq:affinity}
A_{ij} = 2A\left(\frac{1}{1+\exp\left(\frac{-D_{ij}}{\tau_A}\right)}-0.5\right)+1,
\end{equation}
where $\tau_A$ is the characteristic time of affinity increase and $A$ is a constant. The form of Eq. \ref{eq:affinity} was chosen so that $A_{ij}$ goes smoothly from $1 \rightarrow A+1$ as $D_{ij}$ goes from $0 \rightarrow\infty $. This effectively changes the parameters in Eq. \ref{eq:sampleforce} (but not in Eq. \ref{eq:maledist}!) for the $\bm{F}_\text{FL}$ case from $F^\text{AT}_\text{FL} \rightarrow A_{ij}F^\text{AT}_\text{FL}$ and from $R^\text{AT}_\text{LF} \rightarrow A_{ij}R^\text{AT}_\text{LF}$. This allows a follower to split farther from the leader it belongs to, without leaving the harem, thus introducing more consistency into the group compositions.

The definition of groups is based on the values $D_{ij}$. Every follower is considered to be in the group of the leader for which the value of $D_{ij}$ is largest for the given follower. The leader--leader interaction differs in one aspect if either of the participating leaders have a group, by effectively increasing the repulsive radius $R^\text{EX}_\text{LL}$ of both leaders fivefold when interacting with each other. As such two leaders can be close to each other only if they don’t each have their own groups. This is reminiscent of the distinction between bachelor groups, where males are close together and harems, where the males are farther apart.

The velocity dependent part is the same for both leaders and followers:
\begin{equation}
F_{\text{int,}\varphi}(\varphi_{ij}) = \frac{-1}{1+\exp(-4(\varphi_{ij}-\frac{\pi}{2}))} + 1,
\end{equation}
which effectively means, that a horse will pay more attention to horses that are in front of it, rather than those that are behind it. The form was chosen because of the saturation properties. The total interaction is thus
\begin{equation}
\bm{F}_{\text{int}}(r_{ij},\varphi_{ij}) = F_{\text{int,}\varphi}(\varphi_{ij})\bm{F}_\text{int,r}(r_{ij}).
\end{equation}

The force $\bm{F}_\text{com}$ keeps the herd roughly together, since if one strays farther than $R_\text{com}$ from the center of mass of the herd it will experience the force
\begin{equation}
\bm{F}_\text{com}(\bar{\bm{r}}-\bm{r}_i,\bar{\bm{v}}) =
F_\text{com}\frac{|\bm{r}_i - \bar{\bm{r}}|-R_\text{com}}{R_\text{com}}\left( \frac{\bar{\bm{r}}-\bm{r}_i}{|\bar{\bm{r}}-\bm{r}_i|} +
\beta\frac{\bar{\bm{v}}}{|\bar{\bm{v}}|}
\right),
\end{equation}
where $\beta$ is parameter that tunes how much the horse is guided in the direction the center of mass is heading and $F_\text{com}$ is the overall strength of the force. Since $R_\text{com}$ is relatively large this force is usually inactive, but will smoothly guide a lost horse back into the herd (adopted from \cite{kunal2010}).

The force $\bm{F}_\text{wall}$ sets the boundary conditions. The herd is confined to a square area defined by the length $D$. This box is impenetrable and horses cannot leave it. For the herd to approach this hard boundary in a realistic way, there is a characteristic distance $R_\text{wall}$ where the force $\bm{F}_\text{wall}$ is turned on:
\begin{equation}
\bm{F}_\text{wall}(\bm{r}_i,\bm{v}_i) = \frac{F_\text{wall}}{2}\left( \sin \left[\pi\left( \frac{R_\text{wall}-d_{iw}}{R_\text{wall}}-\frac{1}{2} \right)\right] + 1 \right)
\begin{pmatrix}
\bm{v}_i\cdot\bm{n}_w \\
\bm{v}_i\cdot\bm{t}_w
\end{pmatrix},
\end{equation}
where $d_{iw}$ is the distance of the horse from the boundary, $\bm{n}_w$ is the normal vector of the boundary and $\bm{t}_w$ is the tangent vector of the boundary, driving the horses smoothly along the wall (adopted from \cite{viragh2014} and \cite{viragh2014corr}).

Initially both leaders and followers are evenly distributed over a square with a linear size of 500, with velocities also randomly distributed.

\subsection{Parameters}
Going, in a na\"{\i}ve way, from the one-type-particle model of \cite{szabo2006} to the two-type-particle model would increase the number of required parameters from 14 to 30 (some parameters are doubled and some are increased fourfold given every possible combination of the particles). By considering that some of these are unnecessary to duplicate (or make four of) our model has 23 parameters. Of these only 7 are relevant in the sense that the formation of meaningful groups is sensitive to their value (parameters that would destroy cohesion even in a one-type-particle model were not taken into account), not considering the size of the herd. A parameter was considered relevant if an increase by twofold or a decrease by half resulted in 0.1\% of followers not being in a group on average (this is less than one per a realization of the model). For a complete list of parameters see Table \ref{tab:parameters}. Parameters were chosen so that cohesive movement occurs and that group formation happens. Except for cases where there was a reason to do otherwise, parameters that could be different for leaders and follower were kept the same. The distances were chosen based on observations, the coefficients of the various forces were chosen so that the phenomenology of the movements resembles that of a real herd. The leaders are slightly faster than followers so that they are able to stay  in front of their harem. It must be noted, that in many cases, leaders in real-life examples may not be at the front of their group, but rather at the side or behind; we elected to use the leading-from-front paradigm for the purpose of simplicity. Other choices pertaining to parameter value selection have been mentioned in the previous section describing the model.

\begin{table}[]
\centering

\begin{tabular}{@{}llll@{}}
\toprule
variable & description & default value & approx. dimensions   \\ \midrule

\multicolumn{4}{c}{relevant variables} \\ \midrule
 $A$ & affinity of followers for leaders & 1.3  & 1.3\\
 $\tau_\text{A}$ & characteristic time of affinity & 500 & \unit[218]{s}\\
  $F_\text{FL}^\text{AT}$ & strength of F-L attraction & 0.03 & \unitfrac[0.0125]{m}{s} \\

   $R_\text{LL}^\text{AT} $ & radius of L-L attraction & 200 & \unit[36]{m} \\
 $R_\text{LL}^\text{EX} $ & radius of L-L repulsion & 15 & \unit[2.7]{m}\\
 $R_\text{LL}^\text{*EX} , 5 R_\text{LL}^\text{EX}$ & -- & 75 & \unit[13.6]{m}\\
 $F_\text{LL}^\text{AT} , F_\text{LF}^\text{AT}$ & strength of L-L and L-F attraction & 0.01 & \unitfrac[0.0042]{m}{s} \\

$N_\text{L}$ & number of leaders & 25 & 25 \\
$N_\text{F}$ & number of followers & 175 & 175 \\

\midrule
\multicolumn{4}{c}{irrelevant variables} \\ \midrule
  $F_\text{FF}^\text{AT}$ & strength of F-F attraction & 0.0002 & \unitfrac[0.000083]{m}{s} \\
   $F_\text{LL}^\text{EX}$&  strength of L-L repulsion  & 2 & \unitfrac[0.83]{m}{s} \\
$v_\text{L}^{0}$ & velocity of leaders & 1 & \unitfrac[0.416]{m}{s}\\
$v_\text{F}^{0}$ & velocity of followers & 0.9 & \unitfrac[0.375]{m}{s}\\
$\tau_\text{L}$ & L velocity relaxation time & 3 & \unit[1.31]{s}\\
$\tau_\text{F}$ & F velocity relaxation time & 1 & \unit[0.44]{s}\\
 $\xi$& strength of the noise  & 0.5 & \unitfrac[0.21]{m}{s}  \\
 $L_\text{LL} , L_\text{LF} ,L_\text{FL} ,L_\text{FF}$& non-interaction distances & 1 & \unit[0.18]{m}\\
 $R_\text{LF}^\text{AT} , R_\text{FL}^\text{AT}  , R_\text{FF}^\text{AT}$ & radii of attraction & 50 & \unit[9.1]{m}\\
 $R_\text{LF}^\text{EX}  , R_\text{FL}^\text{EX} , R_\text{FF}^\text{EX}$ & radii of repulsion & 5 & \unit[0.9]{m}\\
$F_\text{LF}^\text{EX}  , F_\text{FL}^\text{EX} , F_\text{FF}^\text{EX}$&  strength of repulsion   & 5 & \unitfrac[2.08]{m}{s} \\
  $R_\text{com}$ & radius of the cohesion force & 250 & \unit[45.5]{m}\\
 $F_\text{com}$ & strength of the cohesion force & 2.5 &\unitfrac[1]{m}{s} \\
 $\beta$ & cohesion force parameter & 0.01 & 0.01\\
 $F_\text{wall}$ & strength of boundary repulsion & 3 & \unitfrac[1.25]{m}{s}\\
 $R_\text{wall}$ & distance of boundary repulsion & 200 & \unit[36.4]{m}\\
 $D$ & linear size of bounding box & 10000 & \unit[1800]{m}\\
 \bottomrule
\end{tabular}

\caption{Table of the parameters of the model grouped according to relevancy in group formation. L and F abbreviate leader and follower respectively. The approximate proper dimensions are based on a comparison with observed horses, see Section \ref{sec:dimscales} for details.}\label{tab:parameters}
\end{table}

\section{Results}

\label{sec:results}

Our model, with the given parameters, produces a cohesive and ordered motion of the entire herd, while forming groups around leaders and also bachelor groups from group-less leaders. This is in qualitative agreement with the actual observed herd moving on an open plane and as an interesting extra phenomenon, our model also includes ``fights'' between leaders for followers. By ``fights'' we mean a situation where two or more leaders without groups get extremely close to one or more followers and after a short time, one of the leaders ``wins'', i.e. a follower is ascribed to be in the leader's group for long enough for it to chase away the other leaders (see Supplementary video 1).

We find that the forming of groups within the herd causes cohesiveness to drop compared to a case without groups. We also find, that in accordance with but with a greater precision than the previous study, the group size distribution of the horses living in the Hortobágy National Park is lognormal. In contrast to this, the current model, based solely on spatial interactions, gives a normal distribution, which implies that spatial interactions alone are not enough to produce the observed group structure.

\subsection{Cohesiveness of the herd}
\begin{figure}
\centering
\includegraphics[scale=1]{figure2.pdf}%
\caption{Starting from a uniform random distribution of positions and velocities (left side) the herd forms groups and exhibits ordered motion (right side). Blue dots represent leaders and red dots represent followers (see Supplementary video 2 for a video example).}
\label{fig:screenshot}
\end{figure}

Starting from uniform random initial positions and velocities of the individuals, after sufficient time, the model develops ordered motion throughout the herd while forming groups and thus arriving at a structured and co-moving herd (see Figure \ref{fig:screenshot} and Supplementary video 2). We assume that during collective migration the horses cannot stop, thus there are two phases of ordered movement:  translational movement, and collective rotation about the -- otherwise slowly moving -- center of mass. Indeed, when it is not possible to stop (e.g. due to fear), but is not feasible or desirable to move the herd as a whole, herding mammals have been observed to rotate around a common point. To measure translational cohesiveness we use the following translational order parameter

\begin{equation}
\Phi_t = \frac{1}{N} \left |\sum_{i = 1}^{N} \frac{\bm{v}_i}{|\bm{v}_i|} \right |,
\end{equation}

and to measure the rotational cohesiveness we introduce the following rotational order parameter

\begin{equation}
\Phi_c = \frac{1}{N}\sum_{i = 1}^N \mathcal{P}\left(\frac{\bm{v}_i}{|\bm{v}_i|}\right),
\end{equation}

where $\mathcal{P}$ denotes projection onto the normal of the line going through $\bm{r}_i$ and $\bar{\bm{r}}$.

Going from a totally disordered translational movement to totally ordered translational movement $\Phi_t$ will grow from 0 to 1, while $\Phi_c$ will move from -1 to 1, as the system moves from a totally ordered rotation around the center of mass in one direction, through no collective rotation to totally ordered rotation in the other direction.

We find that the system, with parameters given in Table \ref{tab:parameters} switches between two modes, one of ordered rotation and one of ordered translational motion (see Figure \ref{fig:motiontypes} and Supplementary video 3 for an example of a transition from rotational to translational motion). Since the horses in the model do not have the capacity to stop, in an event of indecision about the direction to move they must rotate about a common axis, namely the center of mass. By averaging over the full length of 1000 runs in total, we find that the rotations have no specific direction, as expected (Mann-Whitney U-test, $p=0.96$ on left-right similarity).

\begin{figure}
\centering
\includegraphics[scale=0.65]{figure3.pdf}%
\caption{The herd as a whole either exhibits an ordered translational motion or rotates around a slowly drifting center of mass. This two different type of motion can be distinguished due to the values of the translational (a) and rotational (b) order parameters. The plots are from the same specific run of the model, with the curves smoothed by a window of $\Delta t = 1000$. The spikes during the translational phase are caused by the confining wall (see Supplementary video 4 for a sample of an interaction with the wall).}
\label{fig:motiontypes}
\end{figure}

\begin{figure}[ht]
\centering
        \includegraphics[scale=1]{figure4.pdf}

\caption{%
The pair-correlation for a herd composed of leaders and followers, but only calculated on the leaders (solid line), and in the case where only leaders are present (striped line). In the latter case the distances are scaled with $R_{LL}^{*EX}/R_{LL}^{EX}$, to compensate for the effect of no leader having a group (leaders without followers can be closer to each other than ones with followers). The main structure of the herd, even with followers, is set by the leaders, but the presence of followers slackens the rigidity of this structure.
}
\label{fig:densitycorr}
\end{figure}

Calculating the pair-correlation function

\begin{equation}
\rho(\bm{r}) = \langle \delta(\bm{r}-\bm{r}_i) \rangle,
\end{equation}
for the leaders in the normal scenario (i.e. where followers are present) and in the scenario where followers are missing, we find that the main structure of the herd is given by the leaders, and introducing followers only slightly loosens this (aside from the fact that it increases the distances between the leaders, see Figure \ref{fig:densitycorr}). We also investigated the effect of introducing followers among the leaders on the order parameter of the translational movement. Comparing $\Phi_t$ (calculated using only the velocities of the leaders in two cases, one where there are only leaders and one where there are also followers) we find that order is decreased when allowing for followers and forming of groups (see Table \ref{tab:orderparam}). This loss in the efficiency of the movement of the herd as a whole points to benefits gained from social groupings outside the paradigm of simple locomotion.

\begin{table}[ht!]
\centering

\begin{tabular}{@{}lllll@{}}

\toprule
&  $\langle\Phi_t\rangle$ & $\langle\Phi_t\rangle$ (only leaders) & $\langle\Phi_c\rangle$ &$\langle \Phi_c\rangle$ (only leaders) \\
\midrule
w/o followers & $0.866 \pm 0.016$ & -- & $-0.002 \pm 0.011$ & -- \\
with followers & $0.608 \pm 0.018$ & $0.633 \pm 0.017$& $0.025 \pm 0.020$ & $0.026 \pm 0.021$  \\
\bottomrule

\end{tabular}

\caption{The translational and rotational order parameters averaged over 120 simulations with standard errors. The duration of the runs were many times longer than the stabilization of groups. The first row is from simulations where only leaders were present, the second row is the full model with followers. In this case the averages were calculated on the whole herd as well as on the leaders only. Adding followers and thus moving in groups decreases the order of translational movement, implying that group formation has  benefits other than increased herd cohesion. Although the herd would rotate often, as  expected, there is no specific direction of the rotation (Mann-Whitney U-test on 1000 runs, where the simulations was terminated at a time not long after stabilization of groups yields a $p = 0.96$ on left-rigth similarity).}\label{tab:orderparam}

\end{table}

\subsection{Group size distribution}
Starting from a uniform random spatial distribution and group-less state, the model, after
sufficient time, will produce co-moving groups based on the relative positions of leaders and followers. The emerging group size distribution is normal, although some leaders and followers may not belong to a group. We define groups by the highest (non-zero) $D_{ij}$ values of the followers, i.e a group consists of the leader and the followers with their highest $D_{ij}$ rating corresponding to this leader. This effectively means that groups are formed by followers spending the most time with a specific leader.  The group size distribution rapidly reaches a close-to-final state and after some time relaxes to the final state (see Figure \ref{fig:haremdistributiontime}). We show the transition by creating a histogram of the group sizes at regular intervals during a simulation and taking the sum of the differences of each respective bin of the histogram in two consecutive measurements, averaged over 1000 independent simulations .

\begin{figure}[ht]
\centering
        \includegraphics[scale=1]{figure5.pdf}

\caption{%
The groups size distribution quickly stabilizes as it is shown by the plot of $\Delta$. To calculate $\Delta$ we create a histogram of the group sizes at regular intervals during a simulation and take the sum of the differences of each respective bin of the histogram in two consecutive measurements. Each point is averaged over 1000 independent simulations.
}
\label{fig:haremdistributiontime}
\end{figure}
On Figure \ref{fig:haremdistribution} we show a comparison of the simulated distribution with real data obtained from a Przewalski horse herd (see \cite{ozogany2014} for details). Since harem sizes gradually change over time among the horses, the real data has been improved by taking into account historical harem size distributions, showing a more clear lognormal distribution than in the previous study of \cite{ozogany2014}. In this previous study a network model was formulated to account for the lognormal distribution of the group sizes, while the current model, based on purely spatial interactions was not able to reproduce this. This indicates that at this level of complexity, it is not possible to reduce social interactions to spatial interactions.

\begin{figure}[ht]
\centering
        \includegraphics[scale=1]{figure6.pdf}

\caption{Comparison of the group size distribution in the model with an empirical one (the group size distribution of Przewalski horses living in the Hortobágy National Park). The empirical distribution follows a clearer lognormal distribution than in \cite{ozogany2014} due to the incorporation of historical data. The distribution obtained from the model is close to normal and is mean-fitted to the empirical distribution. We attribute the difference to the fact that the social interactions of horses are too complex to capture in purely spatial interactions.
}
\label{fig:haremdistribution}

\end{figure}

\subsection{Dimension scales}
\label{sec:dimscales}
Horses usually travel by walking, which is roughly around $\unitfrac[1.5]{km}{h}$ based on our aerial observations averaged over several minutes. In this model $v^0_L$ is the corresponding parameter of the walking speed. To compare our model's length scale with that of reality we have calculated the pair-correlation function of the of the wild horses by using aerial pictures of the real herd and that of the herd in our model and compared the first peaks. This roughly equates the arbitrary length unit of our model to $\unit[0.18]{m}$ in reality (c.f. $R_\text{LF}^{EX}$ in Table  \ref{tab:parameters}). From this we can calculate that the arbitrary time unit of our model is roughly equal to $\unit[0.44]{s}
. This puts $\tau_\text{L}$ and $\tau_\text{F}$ at about reaction time ($\unit[0.5-1.5]{s}$), $\tau_\text{A}$ to about 3 and a half minutes, and the emergence of a coherent collective motion, with stable harems to slightly less than $\unit[10]{minutes}$. Since $\tau_\text{L}$ and $\tau_\text{F}$ both characterize a fast cognitive process it is not unrealistic that the characteristic times are on the scale of reaction times. Since in wild horses the groups do not form from randomly distributed individuals spontaneously, but rather evolve in an already laid down social context, the time needed for group formation is not readily comparable to that of the real herd. On the other hand, for a group of 200 unfamiliar individuals, where leaders are already appointed and everybody is actually already moving, the 10 minutes seems like a reasonable time for group formation (the authors' personal experience with spontaneous group formation in human groups of comparable sizes would allow for even longer times).