Interaction Mechanisms Quantified from Dynamical Features of Frog Choruses
Kaiichiro Ota Ikkyu Aihara Toshio Aoyagi
nlin.AO q-bio.NC





Interaction Mechanisms Quantified from Dynamical Features of Frog Choruses



Kaiichiro Ota^1, 2, 
Ikkyu Aihara^3, ∗, 
Toshio Aoyagi^2, 4



1 Cybozu, Inc., Tokyo, Japan

2 JST CREST, Tokyo, Japan

3 Graduate School of Systems and Information Engineering, University of Tsukuba, Tsukuba, Japan

4 Graduate School of Informatics, Kyoto University, Kyoto, Japan

∗ E-mail: aihara@cs.tsukuba.ac.jp

 


Key Words: Acoustic communication, Nonlinear dynamics, Bayesian approach, Synchronization, Japanese tree frogs












§ ABSTRACT


Interaction mechanism in the acoustic communication of actual animals 
is investigated by combining mathematical modeling and empirical data.
Here we use a deterministic mathematical model (a phase oscillator model)
to describe the interaction mechanism underlying the choruses of male Japanese tree frogs (Hyla japonica)
in which the male frogs attempt to avoid call overlaps with each other due to acoustic communication.
The mathematical model with a general interaction term is identified by a Bayesian approach 
from multiple audio recordings on the choruses of three male frogs.
The identified model qualitatively reproduces the stationary and dynamical features of the empirical data, 
supporting the validity of the model identification.
In addition, we quantify the magnitude of attention paid among the male frogs from the identified model, 
and then analyze the relationship between the attention and behavioral parameters by using a statistical model.
The analysis demonstrates the biologically valid relationship 
about the negative correlation between the attention and inter-frog distance, 
and also indicates the existence of a behavioral strategy that 
the male frogs selectively pay attention towards a less attractive male frog 
so as to utilize the advantage of their attractiveness for effective mate attraction.







































































§ INTRODUCTION


Animals show various types of behavior in the form of a swarm.
For instance, fish and birds construct a robust and flexible flock <cit.>.
To maintain  the flock, they need to synchronize their velocity and direction with each other. 
On the other hand, various animals (e.g., mammals, birds, anurans, and insects) 
aggregate in a certain area and utilize acoustic signals for communication and also mate identification 
<cit.>.
Experimental studies demonstrate that these animals tend to alternate their acoustic signals with each other 
<cit.>.
Because such alternating behavior reduces the acoustic interference of their signals,
they are likely to effectively communicate or make themselves more conspicuous within a swarm.
Thus, synchronization and alternation are abundant in an animal swarm, 
and can work as an important indicator that determines the quality of their behavior.


To synchronize or alternate behavior, animals must recognize a specific target in a swarm.
Such selective attention is reported in various systems.
For example, humans pay attention to one of talking people in noisy environment like a party <cit.>; 
fish and birds attend their neighbor in a flock <cit.>;
bats pay attention to specific targets during a prey capture <cit.>;
male insects and frogs call alternately with their neighbor <cit.>.
To understand the roles of the selective attention, it is essential for us 
to quantify their interaction mechanisms 
and compare it with the behavior of an individual animal.



This study aims to quantify interaction mechanisms inherent in the choruses of male Japanese tree frogs (Hyla japonica) 
by identifying a deterministic mathematical model (a phase oscillator model) from empirical data. 
Japanese tree frogs are observed widely in Japan, and breed mainly at a paddy field from April to July <cit.>.
The male frogs form a lek, and produce successive calls to attract conspecific females.
Indoor experiments with multiple male Japanese tree frogs demonstrate
that they avoid call overlaps with each other in the forms of anti-phase synchronization of two frogs,
tri-phase synchronization of three frogs, 
and clustered anti-phase synchronization of three frogs
<cit.> (see Figure <ref> for examples of three frogs).
Moreover, field observations reveal two-cluster synchronization with a larger number of the male frogs 
in which each pair of neighbors tend to call alternately in their natural habitat <cit.>.
Because the temporal overlap of acoustic signals generally mask the information included in each call,
these alternating behavior would be important for the male frogs to 
effectively advertise themselves towards a conspecific female 
<cit.>.







This paper is organized as follows:  
(1) a deterministic mathematical model (a phase oscillator model <cit.>) is identified 
for each pair of male Japanese tree frogs from our empirical data according to a Bayesian approach, 
(2) the attention paid among the male frogs is quantified on the basis of the identified model,
and (3) the relationship between the quantified attention and behavioral parameters is analyzed by using a statistical model.






§ RESULTS


To investigate interaction mechanisms inherent in the acoustic communication of actual animals, 
we analyze empirical data of male Japanese tree frogs that were obtained 
from our previous experiments and  data analysis <cit.> 
(see Sec. <ref> for details). 
In each experiment, we randomly captured three male frogs at a field site, and then placed them along a straight line at the interval of 50 cm in a stationary room. 
Spontaneous calling behavior of the male frogs was recorded by three microphones that were placed in the vicinity of each frog.
The audio data were analyzed by the method of independent component analysis,
and were separated into call signals of respective frogs.
Here we use four datasets of the audio data in which all the frogs stably called more than 1400 times 
in four hours, 
allowing us to precisely estimate a phase oscillator model by utilizing the large sample size of call timing.




Figure <ref> explains how we identify a phase oscillator model from the empirical data.

A phase oscillator model is a well-known deterministic mathematical model 
derived from simple assumptions about periodicity and interaction <cit.> 
(see Sec.<ref> for the details of a phase oscillator model)
that are essential for the acoustic communication of animals.
Various experimental and theoretical studies have shown that 
a phase oscillator model can qualitatively reproduce the synchronization phenomena in actual systems 
<cit.> including alternating chorus patterns of male Japanese tree frogs <cit.>.

In this study, we first calculate a phase ϕ_n(t_i) (n=1, 2, 3) with discrete time t_i for the nth frog 
from the separated audio data according to Eq. (<ref>).
Then, a Bayesian method <cit.> is applied to the time series data of ϕ_n(t_i) 
so as to estimate the parameters of a phase oscillator model 
(i.e., an interaction term Γ_nm(ϕ_n(t_i)-ϕ_m(t_i)), a natural frequency ω_n, 
and a noise intensity σ_n). 
Figure <ref> shows a representative result of the model identification that is obtained 
from a single dataset with three frogs.
Green bold lines represent the interaction terms with the 95% confidence interval. 





To confirm the validity of the model identification, 
we perform numerical simulation by using the identified model and compare it with empirical data.
Note that we set only one parameter σ_n to be slightly larger than its estimated value throughout the following analysis  
because this parameter was very likely to be estimated at a smaller value 
because of the intermittency of the chorus over a long time scale (see Sec. <ref> for details) 
that is out of the scope of the phase oscillator model.
Figure <ref>A represents the scatter plot of phase differences 
that are obtained from numerical simulation of the identified model.
By contrast, Figure <ref>B shows the scatter plot of phase differences 
that are directly calculated from our empirical data. 
Here we plot the phase differences only when one of the phases hits 0, 
which is consistent with our method for calculating a phase difference from discrete call timing <cit.>.
The comparison of Fig. <ref>A and B demonstrates that the identified model can qualitatively reproduce the experimental result, 
supporting the validity of our model identification.
In addition, our empirical data of Fig. <ref>B shows a complex transition 
among clustered anti-phase synchronization and tri-phase synchronization, 
which is consistent with our previous study <cit.>.
To  investigate the mechanism of such a transition, 
we further analyze the relationship between a phase difference ϕ_n-ϕ_m 
and its differential d(ϕ_n-ϕ_m)/dt by using the identified model (see Sec.<ref>).
Figure <ref> demonstrates that d(ϕ_1 - ϕ_3)/dt and d(ϕ_2 - ϕ_3)/dt have critical states  
while d(ϕ_1 - ϕ_2)/dt has equilibrium states.
Because the transitions among multiple states can be driven by added noise under the existence of critical states 
in various dynamical systems, 
this property about critical states is likely to be the origin of the transitions 
observed in the empirical data of the frog choruses.




Next, we quantify the magnitude of the attention paid among the male frogs based on the result of the model identification.
The stationary distribution of each phase difference is first calculated by numerically solving the Fokker-Planck equation
of the identified model (see Sec. <ref>).
Figure <ref> demonstrates that some distributions have a sharp peak.
For instance, the distribution of ϕ_2 - ϕ_1 has a sharp peak around π, 
meaning that the 2nd frog attempted to call alternately with the 1st frog.
Then, we calculate the Kullback-Leibler divergence of the calculated distribution 
from uniform distribution for each pair of the male frogs 
(see Sec. <ref>), 
which corresponds to the quantification of the attention.
The right panel of Fig. <ref> shows the result of the quantification 
in which line width represents the magnitude of the attention. 
Consequently, an asymmetric structure is observed in the attention paid among the male frogs.
For example, the 1st frog pays strong attention to the 2nd frog 
while the 2nd frog pays just weak attention to the 1st frog.



To further examine the validity of the model identification, 
we analyze the relationship between the attention and behavioral parameters 
by using a generalized linear mixed model (see Sec. <ref>). 
In this analysis, the magnitude of the attention is treated as a response variable,
and three behavioral parameters (i.e., an inter-frog distance, the probability for taking a chorus leader (leader probability), and an inter-call interval)
are treated as explanatory variables of fixed factors. 
These variables are calculated from four datasets with 12 frogs 
in which all the frogs stably produced calls. 


Figure <ref> shows the posterior mean and the 95% confidence interval of the coefficients of respective fixed factors. 
This analysis demonstrates that 
(1) the inter-frog distance has a negative effect on the attention, 
meaning that male frogs pay more attention to their neighbor,
and also indicates that (2) the leader probability and the inter-call interval have negative and positive effects on the attention, respectively.





§ DISCUSSIONS


In this study, we quantify interaction mechanisms in the chorus of male Japanese tree frogs 
by identifying a phase oscillator model from empirical data.
The identified model qualitatively reproduces the stationary and dynamical features of the frog choruses, 
which supports the validity of our model identification.
Then, the magnitude of the attention paid among the male frogs is quantified on the basis of the identified model.


To our knowledge, this is the first study that shows the evidence of selective attention inherent in an animal chorus  
by combining empirical data with a deterministic mathematical model.



The relationship between the attention and the behavioral parameters (Figure <ref>)
gives perspectives on the choruses of male Japanese tree frogs.





  * Inter-frog distance: 
Our analysis using a statistical model demonstrates the negative relationship between the attention and the inter-frog distance, 
which means that male Japanese tree frogs pay more attention to their neighbor.
This is consistent with previous studies reporting that a neighboring pair of males alternate their calls 
in various species of frogs and insects <cit.>.
Because sounds attenuate depending on distance, the calls of a neighboring pair should arrive at a higher intensity than the calls of a distant pair.
We speculate that such a sound attenuation depending on distance 
is the origin of the negative relationship between the attention and the inter-frog distance. 



  * Inter-call interval and leader probability:
In various species of frogs and insects, females prefer a conspecific male 
that leads their chorus <cit.> and produces calls at a higher repetition rate <cit.>.
Our analysis indicates that male Japanese tree frogs pay more attention to a male that rarely leads their chorus and produces calls at a lower repetition rate.
These results suggest that the male frogs pay more attention to a less attractive male and then call alternately with him. 
Given that alternating chorus patterns can reduce the acoustic interference of their calls <cit.>, 
this feature would be important for male frogs to effectively advertise themselves towards conspecific females
by utilizing the advantage of their attractiveness over a neighboring male.



Thus, our analysis is likely to show the important features in frog choruses, 
which is valid in the contexts of the acoustic communication of male frogs as well as the strategy for mating.
However, it should be noted that 
the 95% confidence intervals of the coefficients of the inter-call interval and leader probability do not include but very close to 0
while the 95% confidence interval of the coefficient of the inter-frog distance is obviously beyond 0 (Figure <ref>). 
This suggests that the effects of the leader probability and inter-call interval on the attention are marginal compared to the effect of the inter-frog distance.



The present methodology is widely applicable to the analysis on various types of communication in animals because of the following reasons: 
(1) a phase oscillator model is derived from the simple assumptions about periodicity and interaction 
both of which are ubiquitous in the communication of animals relying on various signals 
(e.g., sounds <cit.>, lights <cit.>, visual display <cit.>, and electric fields <cit.>), 
(2) our methodology only requires the timing of signal emissions so as to identify the mathematical model, 
(3) our methodology allows us to separately identify the unidirectional interaction term of a phase oscillator model,
and (4) the identified model qualitatively reproduces the stationary and dynamical features of empirical data.
An important point is that the interaction term of the phase oscillator model varies 
depending on the value of the phase difference that corresponds to the change in the inter-signal interval among individual animals.
This study demonstrates that such a dynamical property of the phase oscillator model can precisely capture not only the stationary distribution of the phase difference 
but also the dynamical feature of the transition among multiple synchronization states. 
We believe that this property of a dynamical model is advantageous compared to traditional methods (e.g., the calculation of the histogram of the inter-signal interval) 
when studying selective attention in animal communication that shows complicated deterministic dynamics. 













§ MATERIALS AND METHODS





 §.§ Estimation of phase dynamics from empirical data



In this study, we use the empirical data of male Japanese tree frogs that were obtained from our previous experiments <cit.>.
In each experiment, we randomly captured three individuals of the male frogs that were calling in a paddy field of Kyoto University (35°01' 57.2"N, 135°47'00.0"E).
Then, we put each frog in a small mesh cage, and placed three cages along a straight line at the interval of 50cm in a stationary room.
After sunset, the spontaneous calling behavior of the male frogs was recorded in the stationary room by three microphones that were placed in the vicinity of each cage.
After the recording, the male frogs were released at the paddy field where they were captured.
The experiment was carried out 44 times in total between 2008 and 2009.
In most trials of the experiments, some of three frogs rarely or did not call at all. 
Because this study aims to estimate the parameters of a mathematical model due to a statistical approach
and the precision of the parameter estimation generally depends on a sample size of the empirical data, 
we choose four datasets in which all the three frogs stably produced calls more than 1400 times in four hours;
the four datasets were obtained from the experiments performed on 26th May, 16th June, 17th June in 2008, and 29th May in 2009, respectively.


All the experiments and collection of Japanese tree frogs were carried out within the facility of Kyoto University 
in accordance with the guideline approved by the Animal Experimental Committee of Kyoto University.
Then, the audio data of the four datasets were analyzed by the method of independent component analysis, 
and were separated into call signals of respective frogs <cit.>.
In addition, we carefully checked the quality of the sound-source separation 
and excluded just three choruses from 143 choruses in which the separation did not work well
<cit.>.
From the separated audio data, we estimated the call timing of respective frogs 
according to the method of Ref.<cit.>.



For each dataset, the call timing of three male frogs is described as T_n, k that represents the timing of the kth call vocalized by the nth frog (n=1, 2 or 3).
By using the call timing T_n, k, we calculate a phase ϕ_n(t_i) for the nth frog at discrete time t_i as follows <cit.>:


    ϕ_n(t_i) = 2πt_i-T_n, k/T_n, k+1-T_n, k.


This calculation is performed for the set of T_n, k and T_n, k+1 
that satisfy the conditions of T_n, k≤ t_i < T_n, k+1 and T_n, k+1-T_n, k < 0.9 sec. 
The first condition (i.e., T_n, k≤ t_i < T_n, k+1) is necessary 
for restricting the phase in the range of 0 ≤ϕ_n(t_i) < 2π.
The second condition (i.e., T_n, k+1-T_n, k < 0.9 sec) is assumed 
because of the property of the frog chorus 
in which male Japanese tree frogs intermittently start and stop their periodic calling behavior over a long time scale.
Namely, each frog periodically calls almost at the interval of 0.3 sec for several tens of seconds, stays silent for several minutes, and then repeats this cycle <cit.>.
Because such an intermittency over a long time scale is out of the scope of the phase oscillator model, we assume the second condition 
and then calculate the phase ϕ_n(t_i) only during the periodic calling behavior. 
Consequently, the phase ϕ_n(t_i) increases from 0 to 2π, 
and then reset to be 0 when the nth frog produces a call during the periodic calling behavior. 









 §.§ Phase oscillator model



To reproduce the stationary and dynamical features in the choruses of male Japanese tree frogs, 
we use a phase oscillator model <cit.> with additive noise as follows:


    dϕ_n(t)/dt   =   ω_n + ∑_m=1, m ≠ n^NΓ_nm(ϕ_n(t)-ϕ_m(t)) + ξ_n(t),


where ϕ_n(t) ∈ [0, 2π) (n = 1, 2, ..., N) is the phase of the nth frog. 
ω_n  is a positive parameter that describes the intrinsic angular velocity of the nth frog.
We then assume that the nth frog produces a call at ϕ_n(t) = 0, which is consistent with the definition of Eq. (<ref>).
Consequently, 2π/ω_n gives the intrinsic inter-call interval of the nth frog
<cit.>. 
Γ_nm(ϕ_n(t)-ϕ_m(t)) is the unidirectional interaction term of the phase oscillator model, 
and is defined as a 2π-periodic function of ϕ_n(t)-ϕ_m(t) <cit.>.
In the context of acoustic communication amomg male frogs, this term represents how the nth frog controls his call timing by responding to the calls of the mth frog 
<cit.>.
ξ_n(t) is the term of additive noise.
We assume that this term is given by independent white Gaussian noise with magnitude σ_n 
that satisfies the relationship <ξ_n(t)ξ_n(s)>=σ_nδ(t-s).
Note that we utilize this term ξ_n(t) so as to better fit the model parameters from empirical data that contain noisy component. 






 §.§ Identification of a phase oscillator model



Given the time series data of the phase ϕ_n(t_i) with discrete time t_i 
that was obtained from our experiments using actual frogs (see Sec. <ref>),
we estimate the unknown parameters of the phase oscillator model according to the Bayesian method of Ref. <cit.>.
In the method, we first define the following likelihood function from the time series data of the phase ϕ_n(t_i):


    L_n = ∏_i 𝒩[ ϕ_n(t_i+1) - ϕ_n(t_i)/t_i+1 - t_i| ω_n 
        + ∑_m=1,m ≠ n^3 Γ_nm(ψ_nm(t_i)), σ_n^2/t_i+1 - t_i],



where 𝒩(x|m, s^2) represents a Gaussian function with mean m and variance s^2, 
and ψ_nm denotes the phase difference ϕ_n - ϕ_m.
Because each interaction term is defined as a 2π-periodic function of the phase difference ψ_nm <cit.>,
it can be expanded into a Fourier series as follows:


    Γ_nm (ψ_nm) = ∑_k=1^M(a_nm^(k)cos kψ_nm + b_nm^(k)sin kψ_nm).








We apply a standard model comparison method to determine the maximum order M in Eq.(<ref>).
That is, we choose M such that a marginal likelihood function ℒ_M is maximized.
Technically, we computed and compared ℒ_M for M = 1, 2, …, 30.
For our data, the maximum order M was chosen between 1 and 9.



Then, the parameters to be estimated are
ω_n, a_nm^(1), …, a_nm^(M), b_nm^(1), …, b_nm^(M) and σ_n,
which are denoted by a shorthand notation c_n.
To estimate c_n in a Bayesian framework, 
we use a reasonable conjugate prior distribution p_prior(c_n),
which is a Gaussian-inverse-gamma function.
Bayes' theorem then gives the posterior parameter distribution 


as follows:

    p_post(c_n) ∝ L_n(c_i) p_prior(c_n).

Because the prior distribution p_prior(c_n) is conjugate to the likelihood L_n(c_i),
we can easily calculate the posterior distribution p_post(c_n) 

(see Ref.<cit.> and its supplemental material for detail).







 §.§ Examination of critical states



Empirical data on the choruses of male Japanese tree frogs demonstrate the complicated transition among multiple synchronization states 
(Fig.<ref>B) <cit.>.
Here we examine the origin of this transition on the basis of the estimated model.
Equation (<ref>) without a noise term 
yields the time evolution of the phase difference ψ_nm≡ϕ_n - ϕ_m as follows:


    G_nm(ψ_nm) ≡d(ϕ_n - ϕ_m)/dt = ω_n - ω_m + Γ_nm(ψ_nm) - Γ_mn(-ψ_nm),


where Γ_nm(ψ_nm) is the maximum a posteriori (MAP) estimation of the interaction term

that was obtained from the analysis of Sec. <ref>.


This function G_nm(ψ_nm) quantifies how the phase difference changes in time 
under the mutual interaction between the nth and mth frogs without the effect of another frog.
If G_nm(ψ_nm) has a zero-crossing (a stable equilibrium) at ψ_nm = ψ_nm^* , 
the phase difference approaches ψ_nm^* and then stays near the value forever.
Such stationary dynamics typically come from strong interaction.
By contrast, if G_nm(ψ_nm) has no zero-crossing but is very close to zero at a certain point ψ_nm = ψ_nm^*  (a critical state), 
the phase difference stays near ψ_nm^* for a certain amount of time and eventually departs the point.
Such critical dynamics typically come from moderate interaction.






 §.§ Quantification of selective attention



To assess selective attention among male frogs, we evaluate the stochastic feature of the estimated model 
by taking the effect of noise into consideration.

For example, to evaluate the attention from the 1st frog to the 2nd frog,
we analyze the following model only with the unidirectional interaction term Γ_12(ϕ_1(t) - ϕ_2(t)):


    dϕ_1/dt = ω_1 + Γ_12(ϕ_1(t) - ϕ_2(t)) + ξ_1(t),
    	
    dϕ_2/dt = ω_2 + ξ_2(t).


These equations describe the situation that the 1st frog exposed to the noise term ξ_1(t) 
pays attention to the 2nd frog according to the interaction term Γ_12(ϕ_1(t) - ϕ_2(t)) (Eq. (<ref>))
but the 2nd frog exposed to the noise term ξ_2(t) does not pay any attention to the 1st frog (Eq. (<ref>)).
Hence, this is a concise mathematical model capturing the attention from the 1st frog to the 2nd frog 
under the effect of the noise ξ_1(t) and ξ_2(t).
Subtracting Eq. (<ref>) from Eq. (<ref>) yields the time evolution of the phase difference ψ_12≡ϕ_1(t) - ϕ_2(t) as follows:



    dψ_12/dt = ω_1 - ω_2 + Γ_12(ψ_12) + ξ_1(t) + ξ_2(t).


Then, we calculate the stationary distribution of the phase difference ψ_12
by numerically solving the Fokker-Planck equation of Eq.(<ref>)
that gives the time evolution of the distribution of the phase difference (i.e., f(ψ_12, t)),


    ∂ f(ψ_12, t)/∂ t
        = -dΓ_12/dψ_12(ψ_12)f
        - [ω_1 - ω_2 + Γ_12(ψ_12)]∂ f/∂ψ_12
        + σ_1^2 + σ_2^2/2∂^2 f/∂ψ_12^2,


until it converges. 
Subsequently, we can calculate the stationary distribution of the phase difference for all the pairs of male frogs, 
and then describe it as f̂(ψ_nm).



The stationary distribution of a phase difference (f̂(ψ_nm)) allows us to quantify the degree of attention.
If there is no attention from the nth frog to the mth frog, the distribution of ψ_nm is almost uniform; 
on the other hand, if the nth frog calls synchronously with the mth frog as the result of paying strong attention to him, 

the distribution has a sharp peak. 

In this study, we quantify such a sharpness of the distribution f̂(ψ_nm) by using the Kullback-Leibler divergence 
from uniform distribution u(ψ_nm) ≡1/2π as follows:


    D_KL(f̂||u) = ∫_0^2πf̂(ψ_nm)logf̂(ψ_nm)/u(ψ_nm) dψ_nm.


Consequently, a nearly uniform distribution is represented by D_KL(f̂||u) ∼ 0
while a sharply peaked distribution is represented by D_KL(f̂||u) that is much larger than 0.






 §.§ Relationship between selective attention and behavioral parameters



To further examine the validity of the model identification, 
we analyze the relationship between selective attention and behavioral parameters.
Here we focus on the following behavioral parameters of male Japanese tree frogs:
(1) an inter-frog distance, (2) an inter-call interval, and (3) leader-follower relationship.
An inter-frog distance represents the distance between each pair of male frogs that was measured in our experiments <cit.> (see Sec. <ref>).

An inter-call interval is calculated as δ T_n, k = T_n, k+1 - T_n, k using the sequences of the call timing T_n, k
only when the condition T_n, k+1 - T_n, k≤ 0.6 sec is satisfied. 
The leader-follower relationship is determined according to the following definition: 
the leader, the 1st follower, and the 2nd follower are defined as the frogs that start calling first, second, and third
within the same chorus, respectively.
Note that male Japanese tree frogs start calling with low-intensity sound, 
making it difficult for us to automatically determine the leader-follower relationship.
Hence, we manually determined the leader-follower relationship by carefully looking at all the separated audio data 
and calculated the probability for taking a chorus reader for each frog.




Next, we analyze the relationship between the magnitude of attention and the three behavioral parameters 
by using a generalized linear mixed model (GLMM).
GLMM is a well-known statistical model that is used in various research areas 
so as to analyze the effects of multiple explanatory variables on a response variable <cit.>.
Here, we treat the magnitude of the attention paid from the nth frog to the mth frog (see Eq.(<ref>)) as a response variable Y_att.
We then treat the three behavioral parameters 
(i.e., the inter-frog distance between the nth and mth frogs, the inter-call interval of the mth frog, and the leader probability of the mth frog) 
as explanatory variables of fixed factors, 
and describe them as X_dis, X_int and X_prob, respectively.
In addition, we treat frog index and experimental date as explanatory variables of random factors 
(ξ_frog and ξ_date) 
because these factors are difficult to be quantified but very likely to affect the calling behavior of male frogs. 
Here, we assume that the random factors ξ_frog and ξ_date follow 
a normal distribution of zero mean value 
with standard deviations of σ_frog and σ_date, respectively.
By using these variables, we construct the following GLMM:


    logα = β_0 + β_disX_dis + β_int X_int + β_prob X_prob + ξ_frog + ξ_date,
    	
    
    	Y_att∼Gamma(αβ, β).


Here, the response variable Y_att is always positive
because of its definition (see Eq.(<ref>)).
To reproduce this feature, we assume that (1) the response variable follows a gamma distribution that always takes a positive value,
and (2) the parameter α of the gamma distribution
(this parameter gives the mean value of the distribution) 
is linked to the explanatory variables with a log function;  
this framework is consistent with a traditional gamma regression with multiple explanatory variables.
We confirmed that there is no multicollinearity among any pairs of the explanatory variables of fixed factors
(the absolute value of Pearson's correlation coefficient is less than 0.36).
Posterior distributions of all the unknown parameters (i.e., β_0, β_dis, β_int, β_prob, β, σ_date, and σ_site)  
were estimated from MCMC samples generated by R Statistical Software (ver. 3.4.2) and Stan (ver. 2.17.2).
Note that we normalized the explanatory variables of fixed factors from 0 to 1
prior to the calculation of the MCMC samples, 

and confirmed the convergence of the MCMC samples by using R̂ with a threshold of 1.01 <cit.>. 
The posterior mean and the 95% confidence interval of β_0, β_dis, β_int, and β_prob
are shown in Figure <ref>.







§ ACKNOWLEDGEMENT

We thank K. Aihara, H.G. Okuno, R. Takeda, T. Mizumoto, H. Awano, K. Itoyama, and Y. Bando 
for their valuable comments on this study.











§ AUTHOR CONTRIBUTIONS

K.O., I.A. and T.A. designed the research;
K.O. and I.A. analyzed the data;
K.O. performed simulation; 
K.O., I.A. and T.A. wrote the paper.




§ COMPETING INTERESTS

We declare we have no competing interests.




§ FUNDING


This study was partially supported by JSPS Grant-in-Aid for Challenging Exploratory Research (No. 16K12396) 
and Grant-in-Aid for Young Scientists (No. 18K18005) to I.A.







































