\section{Introduction}
Human interactions are of particular importance in the dynamics of disease spreading and social contagions. Complex networks and human dynamics provide two different perspectives on the understanding of such interactions. Many appealing topological structures, differing from random graphs, have been uncovered in social interactive networks, such as small world \cite{Watts-Strogatz-1998-Nature}, scale-free degree distribution \cite{Onnela-Saramaki-Hyvonen-Szabo-Menezes-Kaski-Barabasi-Kertesz-2007-NJP}, community \cite{Palla-Barabasi-Vicsek-2007-Nature, Ahn-Bagrow-Lehmann-2010-Nature}, weak tie \cite{Onnela-Saramaki-Hyvonen-Szabo-Lazer-Kaski-Kertesz-Barabasi-2007-PNAS}, to list a few. The bursty characteristics of the temporal patterns are revealed in various types of human interacting activities from the perspective of human dynamics, such as email communications \cite{Barabasi-2005-Nature,Malmgren-Stouffer-Motter-Amaral-2008-PNAS}, short-message correspondences \cite{Hong-Han-Zhou-Wang-2009-CPL,Wu-Zhou-Xiao-Kurths-Schellnhuber-2010-PNAS,Zhao-Xia-Shang-Zhou-2011-CPL}, cellphone conservations \cite{Candia-Gonzalez-Wang-Schoenharl-Madey-Barabasi-2008-JPAMT,Karsai-Kaski-Barabasi-Kertesz-2012-SR}, letter correspondences \cite{Oliveira-Barabasi-2005-Nature,Li-Zhang-Zhou-2008-PA,Malmgren-Stouffer-Campanharo-Amaral-2009-Science}, and so on. Because static networks are limited to the description of sequences of instantaneously interacting links, temporal networks have received considerable research attention recently \cite{ZhaoJohCrSor10,Holme-Saramaki-2012-PR,Pan-Saramaki-2011-PRE}. Indeed, temporal networks have the advantage that temporal patterns of interacting activities can be accounted for on each individual node. This suggests that we should pay more attention to the analysis of temporal patterns of human activities at the individual level rather than at the population level.

Overwhelming evidence shows that human activities differ from Poisson processes and exhibit burst patterns in which the high activity periods are alternated with the low activity periods. Such burst patterns will give rise to fat tailed distributions of inter-event intervals. Several mechanisms have been proposed to explain such fat tailed distributions, including priority-queuing processes driven by human rational decision making \cite{Barabasi-2005-Nature,Vazquez-2005-PRL,Vazquez-Oliveira-Dezso-Goh-Kondor-Barabasi-2006-PRE,Maillart-Sornette-Frei-Duebendorfer-Saichev-2011-PRE}, Poisson processes modulated by circadian and weekly cycles
\cite{Malmgren-Stouffer-Motter-Amaral-2008-PNAS,Malmgren-Stouffer-Campanharo-Amaral-2009-Science}, adaptive interests \cite{Han-Zhou-Wang-2008-NJP}, heterogeneity in the size distribution of social influences \cite{Saichev-Sornette-2013-PRE}, and so on. However, not all of the activities in our daily life can be interpreted by one of the proposed mechanism alone. For instance, the temporal patterns of  short message communications are dominated by both Poissonian task initiations and decision making for task execution \cite{Wu-Zhou-Xiao-Kurths-Schellnhuber-2010-PNAS}. Moreover, the identification of the burst patterns plays a crucial role in understanding the memory behavior in human activities.  Karsai et al. utilized an inter-event time threshold to locate the burst clusters and demonstrated that the deviation from an exponential distribution for the number of events in a burst cluster serves as an indicator of the memory behavior in the timing of events \cite{Karsai-Kaski-Barabasi-Kertesz-2012-SR,Karsai-Kaski-Kertesz-2012-PLoS1,Wang-Yuan-Pan-Jiao-Dai-Xue-Liu-2015-PA}. By defining that each event inside a burst, except the first and last ones, must have at least $k$ neighbor events within $\Delta t$, \cite{Quadri-Zignani-Capra-Gaito-Rossi-2014-PLoS1} proposed a density-based algorithm to detect the burst clusters in the multi-event series. We also need to mention that, in the field of neuroscience, there are many methods to detect the burst patterns in the spike trains, such as Poisson Surprise method \cite{Legendy-Salcman-1985-JN}, Rank Surprise method \cite{Gourevitch-Eggermont-2007b-JNM,Gourevitch-Eggermont-2007a-JNM}, Robust Gaussian Surprise method \cite{Ko-Wilson-Lobb-Paladini-2012-JNM}, and so on.

The enormous increase of popularity and the use of cellphones not only provides convenient high-efficiency real-time communication between people; it also facilitates information spreading and social contagion. One of the crucial ingredients needed to understand the dynamics of information diffusion is the set of temporal patterns of individual cellphone conversation activities.As an illustration of the motivation for testing how population wide analyses inform on individual level studies and vice-versa, consider the financial and economic crisis of 2008, which led to strong criticisms of the economic profession for its use of the concept of the representative agent. While convenient mathematically, this assumption of homogeneisation is grossly insufficient to capture the heterogeneity of economic agent preferences, which is underpinning the complex nonlinear dynamics of the socioeconomic world. We believe that a number of high profile analyses have sinned in the same direction by analyzing human activities at the population level and claiming the existence of non-Poisson behaviors with supposed long-memory and fat-tailed statistics at the individual level.

In order to understand in greater depth the temporal patterns of cellphone conversation activities at the individual level and to decipher their immediate consequences for the simulation of individual call activities, we study the interacting activities of cellphone voice communications for each of 73,339 cellphone users. This research is a continuous work based on the results of \cite{Jiang-Xie-Li-Podobnik-Zhou-Stanley-2013-PNAS}. In that paper, an automatic fitting technology is developed to find the distribution of the waiting time between consecutive outgoing calls of each individual in our selected sample, which is made of 100,000 individual cellphone users each having at least 997 outgoing calls. It is found that there are 3,464 individuals with a power-law distribution of inter-call duration and 73,339 individuals with a Weibull distribution of inter-call duration. The users with a power-law duration distribution exhibit extreme calling behaviors and can be classified as robot-based callers, telecom fraud, and telephone sales. The users with a Weibull duration distribution correspond to the ordinary cell phone customers. Our analysis here is trying to understand the calling behaviors of the 73,339 individuals more deeply and find the key ingredients to model their individual call dynamics. We find that the calling events of the ordinary users can be separated into independent bursts. And we also find strong evidence that the bursts and the call events within bursts are dominated by two independent Poisson processes. We propose a minimal model by modulating the two Poisson processes by a Markov chain, which is in good agreement with our empirical results for about half of the users.

Our results show that population level analyses are misleading, quite analogously to the misguided
use of a representative agent in economics, as mentioned above. Mixing the dynamics of a heterogeneous population of agents gives the appearance of long-memory dynamics, which is just an artefact of the mixing process. In other words, people are different and cannot be represented by a single representative law. This is what we demonstrate in the present article. Our finding underlines the need to fully account for heterogeneous preferences and behaviors in human population to understand and model adequately the burst patterns in individual activities.

\section{Data Sets}

Our data include the call activities of 5,921,696 individuals in Shanghai. The records cover two separated periods, 28 June 2010 to 24 July 2010 and 1 October 2010 to 31 December 2010. By excluding the days in which there are missing records, we have a total of 109 days. Each record entry includes the information on caller number, callee number, call starting time, and call length. The caller and callee number is encrypted in order to protect personal privacy.  Note that the calls with zero durations are discarded . 73,339 individuals in the top 100,000 most active individuals have been selected for our study based on the following criteria: these 73,339 chosen users have the same truncated Weibull distribution and exhibit normal behaviors of phone call activities \cite{Jiang-Xie-Li-Podobnik-Zhou-Stanley-2013-PNAS}. It would be interesting to see whether the call activities of such phone users are governed by the same mechanism.

We only take the outgoing calls to calculate the inter-call durations, because the incoming calls must be initiated by someone else and are byproducts of the outgoing calls. Note also that we are able to replicate the call events in the phone communication network only with respect to the process of individual outgoing calls. The key ingredients to model individual outgoing calls are: (1) when to initiate a call and (2) whom to call. Investigating the waiting time between consecutive outgoing calls is the first step towards understanding the first key ingredient when to initiate a call. We define the inter-call durations as follows. For a given series of outgoing call events $\{c(t_i):i=1,\cdots,n\}$, where $c(t_i)$ is the event of the $i$-th out-going call and $t_i$ is its starting time, the inter-call duration is defined as the time that elapses between two consecutive outgoing calls and it is calculated via $d_i = t_{i} -t_{i-1}$. The discontinuous recording days and the inactive periods in the night result in very large inter-call durations. Therefore, we restrict the durations to a period of one day (the typical human circadian rhythm). Each day is partitioned at 4:00 A.M., at which the lowest call volume in a 24-hour period is observed. This also allows us to take into
account the people who go out and stay awake later as well. Our consideration on the intraday inter-call durations are further supported by the following three arguments. (1) Human activities exhibit strong circadian rhythms \cite{Song-Qu-Blumm-Barabasi-2010-Science}. On each day, people repeat their daily activities with strong regularity, typically going to work in the morning and coming back home in the evening. (2) 99.2\% of the intercall durations from the last call before 4 AM to the first call after 4 AM are greater than 4 hours, which supports our conjecture that the intercall durations crossing the time point 4:00 AM is very large. Obviously, in more than 99\% of the cases, this inactive period has no impact on separating the call events into burst sizes and non-burst sizes under the situation where we do not cut the call periods into daily pieces. (3) The likely reason for observing such large values of intercall durations is that most of the cellphone users are sleeping at that time. Apparently, the underlying dynamics here is different from that in the daytime.

\section{Empirical Results}

We previously reported that the inter-call durations of each individual in our analyzing sample conform to the Weibull distribution \cite{Jiang-Xie-Li-Podobnik-Zhou-Stanley-2013-PNAS}. This implies that their call dynamics may be governed by the same mechanism. In order to uncover the underlying mechanism, we focus our efforts on the understanding of burst patterns in the individual calling process, which reveals the switching between the periods of high activity and low activity. Due to the limited usefulness of the Hurst exponent and of the autocorrelation function to describe dependence in such processes with a fat-tail inter-event distribution \cite{Yang-Jiang-Wang-Wang-Fang-2012-AHN}, we use the distribution of the number of events in high activity and low activity periods to indicate the presence of the correlated behaviors in the event series for each individual. Specifically, a series of call events of one given individual is defined as $\{c(t_i):i=1,\cdots,n\}$, where $c(t_i)$ is the event of the $i$-th out-going call and $t_i$ is its starting time. A high activity period is defined as a cluster of consecutive calls
delineated by the $j$-th and the $k$-th ones,
$s_h=\{c(t_i): 1 \le j < i \le k \le n\}$, such that each such $i$-th call follows the pevious call within a short time interval $\Delta t$.
In other words, all events $i$ in that period must meet this requirement $t_i - t_{i-1} \le \Delta t$. Note that $\Delta t$ is a parameter of the empirical analysis. \Fref{Fig:Show:CallEvents} shows the sequence of outgoing calls for one typical cellphone user on June 28, 2010. The red vertical lines correspond to the outgoing calls that fall within $\Delta t =500$ seconds of previous calls. The number of events in each high activity period $e_h$ is
obtained by counting the total number of calls belonging to the same high activity period. Comparing with the burst size $E$ defined in Ref.~\cite{Karsai-Kaski-Barabasi-Kertesz-2012-SR}, we have $e_h = E-1$. According to the definition \cite{Karsai-Kaski-Barabasi-Kertesz-2012-SR}, the burst size $E$ should conform to the exponential distribution if the call events are independent and the deviation from an exponential distribution diagnoses the existence of a correlation between consecutive events. This conclusion remains valid for our $e_h$ defined in the present paper.

\begin{figure}[htbp]
 \centering
 \includegraphics[width=16cm]{Fig1.eps}
 \caption{\label{Fig:Show:CallEvents} (color online). Sequence of outgoing calls for one typical cellphone user on June 28, 2010. The red (respectively black) vertical lines correspond to the outgoing calls that follow their previous outgoing calls
 within a waiting time less (respectively greater) than $\Delta t = 500$ seconds.}
\end{figure}

Complementarily, the low activity period is defined as the set of events $s_l=\{c(t_i): 1 \le j < i \le k \le n\}$ in that period separated from their nearest previous call by a duration greater than $\Delta t$, which requires all the events in the same low activity period to obey $t_i - t_{i-1} > \Delta t$. These events are represented by black vertical lines in \fref{Fig:Show:CallEvents}. The number of events in each low activity period $e_l$ is obtained by counting the total number of events belonging to the same low activity period. Again, deviations of the distribution of $e_l$ from an exponential would imply the existence of a dependence in the timing of these events.

\begin{figure}[htbp]
  \centering
  \includegraphics[width=16cm]{Fig2.eps}
  \caption{\label{Fig:PDF:NumEvents:AllUsers} (color online).Probability distribution of the number $e_h$ of call events in high activity period at the population level. (A) Probability distribution of $e_h$ calculated with threshold $\Delta t = 100$, $300$, $600$, and $1000$ seconds. (B) Probability distribution of $e_h$ calculated with threshold $\Delta t = q \langle d \rangle$, in which $q=0.5$, $1$ and $3$ and $\langle d \rangle$ is the mean inter-call duration. Note that, for each specific user, we use his/her mean inter-call duration to calculate $\Delta t$. (C) Probability distribution of the average inter-call duration $\langle d \rangle$ for the whole sample. }
\end{figure}

\Fref{Fig:PDF:NumEvents:AllUsers} A shows that the probability distributions of $e_h$ corresponding to different values of $\Delta t$ for individuals in our sample all exhibit clean power-law behaviors with an exponent $\approx 3.6 \pm 0.4$, which is in agreement with the power-law distribution of burst size $E$ with an exponent around $4.1$ shown in Fig.~2 A in Ref.~\cite{Karsai-Kaski-Barabasi-Kertesz-2012-SR}. Because the call frequency varies a lot for different individuals (see the distribution of individual mean inter-call duration in \fref{Fig:PDF:NumEvents:AllUsers} C), using the same value of $\Delta t$ to calculate $e_h$ for different users may give different value of $e_h$ at different scales, which could be a possible source of the power-law distribution. In order to remove such influence, we express the threshold $\Delta t$ in units of $\langle d \rangle$, $\Delta t = q \langle d \rangle$, where $\langle d \rangle$ is the average value of inter-call durations for each specific cell phone user. The probability distributions of $e_h$ corresponding to different values of $q$ for all individuals in our sample are plotted in \fref{Fig:PDF:NumEvents:AllUsers} B. Again, we see very nice power-law behaviors with an exponent around 5.3 in the tail. \Fref{Fig:PDF:NumEvents:AllUsers} C plots the probability distribution of the mean inter-call durations for the individuals in our sample. The mean duration spans a range from 50 to 10,000 seconds, which indicates a strong heterogeneous call frequency for the individuals in our sample.
\begin{figure}[htbp]
  \centering
  \includegraphics[width=16cm]{Fig3.eps}
  \caption{\label{Fig:ICD:Memory} (color online). Complementary cumulative distribution of the number of events $e_h$ in high activity period (A-C) and the number of events $e_l$ in low activity period (D-E), for three individuals, 6493248, 674842, and 961831 (one column for each individual). Each marker corresponds to a different $q$, defined such that $\Delta t = q \langle d \rangle$, where $\langle d \rangle$ is the average value of individual inter-call durations and $\Delta t$ is the threshold used to define high activity period and low activity period.}
\end{figure}

For all our different specifications of the threshold $\Delta t$,  we find that the power-law distributions of the number $e_h$ of events in high activity are robust. Since our focus here is to understand the individual calling patterns, it would be interesting to check whether the distributions of $e_h$ for each individual exhibit the same pattern as that of the aggregated sample. We further scan the values $q = 0.5$, $1$, and $2$ to investigate the influence of $\Delta t$ on the distributions of the number $e_h$ of events of in high activity periods and the number $e_l$ of events in low activity periods for each individual in our sample. \Fref{Fig:ICD:Memory} shows the distributions of $e_h$ and $e_l$ for three users. Surprisingly, we find that both the number of events in high activity periods and low activity periods conform to an exponential distribution for different values of $q$. These results are in striking contrast to the power law distributions of $e_h$ in \fref{Fig:PDF:NumEvents:AllUsers} B at the population level. Such exponential distributions of $e_h$ and $e_l$ at the individual level supports strongly the hypothesis of no temporal dependence between individual calls in both high activity periods and low activity periods. Our results also provide a caution that building individual model based on empirical results at the population level may be misleading.

\begin{figure}[htbp]
  \centering
  \includegraphics[width=16cm]{Fig4.eps}
  \caption{\label{Fig:PDF:BurstyPeriods} (color online). Complementary cumulative distributions of inter-call durations within bursts  (A-C) and time intervals between consecutive bursts (D-E) for three individuals (one column for each individual). Each color corresponds to a different $q$, defined such that $\Delta t = q \langle d \rangle$, where $\langle d \rangle$ is the average value of individual inter-call durations and $\Delta t$ is the threshold used to define high activity period and low activity period.}
\end{figure}

As illustrated in \Fref{Fig:Show:CallEvents}, given a threshold $\Delta t$, we have classified all outgoing call events into two groups, the events in high activity period (red lines) and the events in low activity period (black lines), according to the waiting time from their previous nearest neighbors. Using this classification, we define a bursty cluster as a set of successive calls, where the first one belongs to the low activity class (black line) and all the others are part of the high activity class (red lines). The cluster ends at the last call in the bursty class that is followed by an event in the low activity period. Such a cluster can be interpreted as a sequence of high activity calls somehow triggered or related to the first one that acts as an initiating event. The clusters found for the call sequence shown in \fref{Fig:Show:CallEvents} are depicted by shadow areas. Note that not all black lines are following by red lines, indicating that not all calls initiate other calls.

For a given configuration of individual call bursts with a specific threshold $\Delta t = q \langle d \rangle$, we investigate the probability distributions of inter-call durations $d_b$ within bursts and of the time intervals $d_n$ between consecutive bursts. \Fref{Fig:PDF:BurstyPeriods} shows these two complementary cumulative distributions for three individuals. One can observe that the distribution of $d_b$ conforms approximately to a right-truncated exponential distribution for different values of $q$ which, together with the exponential distribution of $e_h$ in high activity periods, provides strong supports to the hypothesis that the calls in burst clusters are generated by a Poisson process. And the distributions of the time intervals $d_n$ between consecutive bursts are close to left-truncated exponential distributions, the more so, the larger $q$ is, which suggests that the bursts are also generated by a Poisson process. The Poisson initiations of bursts are also observed in individual short message communications \cite{Wu-Zhou-Xiao-Kurths-Schellnhuber-2010-PNAS}.

Motivated by the observations in \Fref{Fig:PDF:BurstyPeriods}, we conjecture the existence of an optimal threshold $\Delta t_{\rm{opt}} = q_{\rm{opt}} \langle d \rangle$, which can be determined as the value that minimizes the residuals of the fits of the empirical distributions by truncated exponential distributions. To determine $\Delta t_{\rm{opt}}$, we perform a Maximum Likelihood Estimation (MLE) of the empirical distributions using the two following expressions:
\begin{itemize}
  \item[(i)] the left-truncated exponential distribution of the time intervals $d_n$ between consecutive bursts is given by
    \begin{equation}
       \label{Eq:BurstyPeriods:LTE}
       p(d_n,\lambda_n) = \lambda_n \exp(-\lambda_n d_n)/ \exp(-\lambda_n \Delta t)~;
    \end{equation}
  \item[(ii)] the right-truncated exponential distribution of the inter-call durations $d_b$ within bursts is
    \begin{equation}\label{Eq:BurstyPeriods:RTE}
       p(d_b, \lambda_b) = \lambda_b d_b \exp(-\lambda_b d_b)/[1-\exp(-\lambda_b \Delta t)]~.
  \end{equation}
\end{itemize}
We use the complementary cumulative distributions $C_{\rm fit}$ derived from the probability distributions (\eref{Eq:BurstyPeriods:LTE} and \eref{Eq:BurstyPeriods:RTE}), to define the fitting residual, as proposed in Ref. \cite{Wu-Zhou-Xiao-Kurths-Schellnhuber-2010-PNAS}, between the empirical and theoretical complementary distribution functions:
\begin{equation}\label{Eq:Separate:Period:Residual}
 R = \frac{ \left\{ \sum_i^g \left[\frac{C_{i, {\rm{fit}}}^n - C_{i,
{\rm{emp}}}^n}{C_{i, {\rm{fit}}}^n + C_{i, {\rm{emp}}}^n}\right]^2 + \sum_i^m
\left[\frac{C_{i, {\rm{fit}}}^b - C_{i, {\rm{emp}}}^b}{C_{i, {\rm{fit}}}^b +
C_{i, {\rm{emp}}}^b}\right]^2  \right\}^{\frac{1}{2}} }{\sqrt{g+m}} ~£¬
\end{equation}
where $g$ and $m$ are the number of observations in the sample of $d_n$ and $d_b$, which are generated by the threshold $\Delta t$.

By varying $q := \Delta t / \langle d \rangle$ from 0.2 to 4 with a step of 0.01, we estimate the fitting residuals $R$ for each $q$ and obtain the optimal threshold as $\Delta t_{\rm{opt}} = q_{\rm{opt}} \langle d \rangle$, with $q_{\rm{opt}}$ being the value that is associated with the minimum fitting residuals $R_{\min}$, as illustrated in \fref{Fig:PDF:BestBurstyPeriods}. We find that the distributions of $d_b$ and $d_n$ can be well described by truncated exponential distributions at the optimal $\Delta t_{\rm{opt}}$. This confirms that individual cellphone call sequences can be generated by two Poisson processes, in which one initiates call events, and the other represents the stochastically activated bursts.

\begin{figure}[htbp]
  \centering
  \includegraphics[width=16cm]{Fig5.eps}
  \caption{\label{Fig:PDF:BestBurstyPeriods} (color online). Results of the optimal threshold for three individuals, 6493248, 674842, and 961831 (one row for each one). (A-C) Plots of the fitting residual \eref{Eq:Separate:Period:Residual} as a function of $q = \Delta t / \langle d \rangle$, where $\langle d \rangle$ is the average value of inter-call durations and $\Delta t$ is the threshold used to define high activity period and low activity period. The marker $\times$ represents the value $(q_{\rm{opt}}$ for which the fitting residual \eref{Eq:Separate:Period:Residual} attains its minimum $R_{\min})$. (D-E) Probability distributions of inter-call durations within bursts, using $\Delta t_{\rm{opt}}$. (G-H) Probability distributions of the time intervals between consecutive bursts for $\Delta t_{\rm{opt}}$.}
\end{figure}

We separate the call events into different bursts according to the best $q_{\rm{opt}}$ and perform four additional KS tests. The corresponding hypothesises are that ($H_1$) the number of events in the high activity periods conforms to an exponential distribution, ($H_2$) the number of events in the low activity periods obeys an exponential distribution, ($H_3$) the inter-call duration inside bursts follows to a right-truncated exponential distribution, and ($H_4$) the waiting time between consecutive bursts complies a left-truncated exponential distribution. At the significant level of 0.01, the hypothesis $H_1$ and $H_2$ cannot be rejected for 66,918 out of 73,339 (about 91.2\%) the individuals, which provides very strong evidence in favor of the Poisson nature of the call events in the high and low activity periods. Furthermore, it is found that the hypothesis $H_3$ and $H_4$ cannot be rejected for 27,387 individuals (about 37.3\% of the whole population) at the significant level of 0.01, in which 26,631 individuals (about 36.3\% of the whole population) also accept the hypothesis $H_1$ and $H_2$ at the same level. Such results support the our conjecture that the call dynamics are dominated by two alternated Poisson processes for about 36.3\% of the individuals in our sample.

\section{Model}

Our empirical results provide a very clear foundation for conceptualising our model, which is constructed on three key ingredients: (1) a Poisson process to trigger the calling burst clusters, (2) a Poisson process to activate the call events inside burst clusters, and (3) an alternation of the two Poisson processes. In some sense, our model is similar to the cascading nonhomogeneous Poisson model of individual email activities \cite{Malmgren-Stouffer-Motter-Amaral-2008-PNAS}, which has two essential ingredients: (1) generating the bursts according to the circadian and weekly cycles of human activity and (2) assigning the number of events into a burst based on a distribution. However, our model has the advantages of simplicity of its components,
smaller computation power needs  for calibration, and simple implementation in simulations. Our model also
benefits from a strong empirical foundation concerning the Poisson nature in individual phone call activities.
In addition to providing an important way for understanding how normal cell phone users activate their personal call events,
our model also offers a wide range of applications for the simulation of phone call traffic and of social contagion within human communication networks.

Based on the model conception, the individual call dynamics is modelled by combining the two Poisson processes within a two-state first-order Markov chain. One Poisson process $\mathcal{P}_n$ with intensity $\lambda_n$ models the activation of bursts. The other Poisson process $\mathcal{P}_b$ with intensity $\lambda_b$ models the flow of call events within bursts. We define $s_n$ as the state in which calls are generated by the process $\mathcal{P}_n$ and $s_b$ as the state in which calls are dominated by the process $\mathcal{P}_b$. The two states swap from one to the other according to a transition (or conditional) probability $p_{ij}$, which can be estimated empirically,
\begin{equation}\label{Eq:TranProb}
    p_{ij} = p(j=s_b|i = s_n) = \frac{p(i=s_n, j = s_b)}{p( i = s_n)}~.
\end{equation}
Here, $p(i = s_n)$ is the probability that a call is in state $s_n$, which can be estimated as the ratio between the total number of calls in state $s_n$ to the total number of calls. Then, $p(i=s_n, j = s_b)$ is the joint probability for the occurrence of two consecutive calls that the first one is in state $s_n$ and the second one is in state $s_b$.

To be able to generate synthetic call sequences, we add a minimum call duration $t_p$ that reflects the characteristic time of a human conversation. We simulate our model as follows:
\begin{itemize}
  \item[(a)] Assuming that the initial state is $s_n$, we generated a call according to the process $\mathcal{P}_n$, in which the call is initiated with the probability of $\lambda_n t_p$;
  \item[(b)] Once a call occurs, the state for generating the next call is determined according to the transition probability $p_{ij}$. For example, if the current state is $s_n$, the next state will be $s_b$ with the probability $p(j=s_b|i=s_n)$, or will remain $s_n$ with probability $p(j=s_n|i=s_n)$;
  \item[(c)] When a new state is determined, the new call is generated according to the Poisson process in that state. If the new state is $s_n$, the new call is initiated with probability of $\lambda_n t_p$. If the new state is $s_b$, the new call is initiated with probability of $\lambda_b t_p$. When the new call is generated, we go back to step (b);
  \item[(d)] When the simulation time reaches our target, the program stops.
\end{itemize}

We estimate the parameters $t_p, \lambda_b, \lambda_n$ and the transition probabilities from the empirical data for each cellphone users and inject them in the above model to simulate their specific call dynamics. \Fref{Fig:ICD:CDF:Data:Model} compares the obtained distributions of inter-call durations obtained by our simulations of three individuals with their empirical counterparts. We find an excellent agreement. We further use the log-residual between the two complementary cumulative distribution, formulated as $\ln(C_M/C_E)$, to quantify the agreement between the model and the real data. \Fref{Fig:ICD:CDF:Data:Model} D plots the log-residual at each inter-call time $d$ for three users. It is observed that our model exhibit a very nice agreement with the real data up to $d<10^4$. We further estimate the probability of log-residual conditioned on the inter-event time for all 73,339 users under consideration. This conditional probability is illustrated in \fref{Fig:ICD:CDF:Data:Model} E. We also plots the average log-residual of all users as a solid line. One can see that our model seems to overestimate the data when $d>10^4$. The possible reason for this could be that our optimized objective function used in the parameters estimation does not aim at minimizing the residuals of the cumulative distribution of the inter-call time between the simulating data and real data, but aiming at minimizing the fitting residuals of two truncated exponential distributions (see \Eref{Eq:Separate:Period:Residual}).

\begin{figure}[htbp]
 \centering
  \includegraphics[width=16cm]{Fig6.eps}
 \caption{\label{Fig:ICD:CDF:Data:Model} (color online). Model comparisons. (A-C) Comparison of the distributions of inter-call waiting times in our simulations and in the real data for three individuals, 6493248, 674842, and 961831. {color{red}{The $p$-value are shown under the individual ID.}} (D) Plot of the log-residual between the model data and the empirical data ($\ln(C_M/C_E)$) with respect to the inter-call duration $d$ for three users. (E) Plots of the average log-residual at each inter-call time for all 7,3339 users. The color bar represents the probability of log-residual conditioned on the inter-call time. All the model data are averaged over 20 simulations.}
\end{figure}

\section{Statistical tests}

If our model can explain the daily individual call activities, at least, the intercall durations from the model should have the same distribution as the empirical data. We perform a strict statistical test to check the agreement of the intercall duration distributions between the empirical data and our proposed model for each mobile phone user. Rather than using the Kolmogorov-Smirnov test, which is a point test strongly influenced by the central part of the distribution \cite{Malevergne-Sornette-2006}, we adopt an integral-based test of the goodness of fit of the empirical distributions by our two-state first-order Markov chain combining the two Poisson processes, defined in terms of the statistics,
\begin{equation}\label{Eq:Model:StatTest:Area}
    A = \int_x |C_E(x) - C_M(x)| {\rm{d}}x~,
\end{equation}
in which $x= \ln d$. This statistic $A$ quantifies the distance between empirical and model distributions as given by the (absolute value of the) area between the two complementary cumulative distributions \cite{Malmgren-Stouffer-Motter-Amaral-2008-PNAS}.

The null hypothesis of the Monte Carlo test is that the intercall durations of the real individual call events and the modelling call events, which is generated by the two-state first-order Markov-chain combining the two Poisson processes, are drawn from the same distribution. The rejection of the null hypothesis implies the failure of our model in the explanation of the individual call activities. If the null hypothesis cannot be rejected, this equals to being against the alternative hypothesis, that both samples of intercall durations do not come from the same distribution, and giving signals that our model can explain the individual call activities at least in the aspect of the duration distribution. Such kind of decision tests, including the KS test and the CvM test, is widely used to determine that one sample of data come from a reference distribution or two samples of data are drawn from the same distribution. Our model is rejected if the absolute difference between $A$ and $\langle A_s \rangle$ is greater than the critical value associating with the probability $1-\alpha/2$. In other words, we can reject the null hypothesis if the $p$-value is less than $\alpha$. The $p$-value is defined by ${\rm{Pr}}(|A_s - \langle A_s \rangle| > |A - \langle A_s \rangle|)$, where $A_s$ is the same as $A$ but obtained from the analysis on the model data, which is generated by our proposed model with the best parameters estimated from the real data.

The procedure to calculate the $p$-value of this statistic $A$ is proceeded as follows:
(a) Generation of the model data with the best parameters estimated from real data;
(b) Calculation of the statistic $A$ between the obtained model distribution and its empirical counterpart;
(c) Regarding the model data as playing the same role as real data, generation of synthetic data with the best parameters estimated from the model data;
(d) Calculation of the statistic $A_s$ between the distribution obtained from the model and its synthetic replicas. One thousand replicas are generated;
and
(e) This allows to evaluate a two-tailed $p$-value, defined as ${\rm{Pr}}(|A_s - \langle A_s \rangle| > |A - \langle A_s \rangle|)$,  where $\langle A_s \rangle$ is the mean value of the synthetic area statistics.

 The probability distribution of $p$-values is illustrated in \fref{Fig:ModPar:PDF:Parameters} A. At the significant level $\alpha = 0.01$, we find that 35'703 individuals have a $p$-value larger than $0.01$. This implies that, for this population, the null hypothesis that our model is the generating process of the call sequences can be rejected at the $99\%$ confidence level.

\begin{figure}[htbp]
  \centering
  \includegraphics[width=16cm]{Fig7.eps}
  \caption{\label{Fig:ModPar:PDF:Parameters} (color online). Results of statistical tests. (A) Cumulative distribution of $p$-value. (B-F) Comparison of the distributions of model parameters and transition probabilities for those individuals who pass the $p$-value test and the others who do not pass. Panels B-F correspond to the distribution of $t_p$, $\lambda_b$, $\lambda_n$, $p(j = s_n|i=s_b)$, and $p(j = s_b|i=s_n)$, respectively. }
\end{figure}

In \fref{Fig:ModPar:PDF:Parameters} B-F, we compare the distributions of the calibrated model parameters $t_p$, $\lambda_b$ and $\lambda_n$ and of the transition probabilities, for the 35'703 individuals who pass the $p$-value test and the other 37'636 individuals who do not. The two distributions of $t_p$, $\lambda_b$ and $\lambda_n$ are statistically identical for the two groups of individuals. This supports the two-state Markov-chain Poisson model of the inter-call durations at the individual level. The broad distributions of parameters highlight the heterogeneity in the individual calling behaviors and support our message about the importance of calibrating at the individual rather than at the population level. Specifically, one can observe that the tails of the distributions of  $\lambda_b$ and $\lambda_n$ are approximately power laws (linear in the log-log representation of panels C and D of \Fref{Fig:ModPar:PDF:Parameters}). The least-square method gives the power-law exponents 4.7 for $\lambda_b$ and 7.6 for $\lambda_n$. This means that the rate of activity of each individual is well-approximated by an exponential distribution but the distribution of the exponential rates across the population is itself a power law. Another intriguing observation is that the power-law exponent of $\lambda_b$ is close to the power-law exponent 5.6 of $e_h$ in \fref{Fig:PDF:NumEvents:AllUsers} B. This can be explained by the mechanism of ``Superposition of Distributions'' described in chapter 14.4.1 of Ref.~\cite{Sornette-2004}, one obtains power law distributions of activities simply by the superposition of many exponential distributions with power law distributions of their characteristic scales. Therefore, we are led to suggest that the power law distribution of the number of events in bursty period at the population level in \fref{Fig:PDF:NumEvents:AllUsers}B is simply the result of this mechanism of a ``Superposition of Distributions''. Our findings thus provide an alternative explanation of the population-level results, which is different from Karsai et al. \cite{Karsai-Kaski-Barabasi-Kertesz-2012-SR,Karsai-Kaski-Kertesz-2012-PLoS1} in which the advanced explanations in terms of long memory behavior and  reinforcement process is proposed at the individual level.

Panels E and F of \fref{Fig:ModPar:PDF:Parameters} identify significant differences between the distributions of transition probabilities. Individuals who pass the test have a relatively smaller probability $p(j = s_n|i=s_b)$ to transit from the bursty state to the non-bursty state and a relatively large transition probability $p(j = s_b|i=s_n)$ to cross from the non-bursty state to the bursty state.

\section{Discussion}

Our results on the statistics of calls differ from that concerning short text messages. In Ref.~\cite{Wu-Zhou-Xiao-Kurths-Schellnhuber-2010-PNAS}, it was found that the initiation of bursts in short messages is Poisson, as for the calls found here. However, the inter-message durations within bursts are distributed as a power law, while we reported an exponential distribution of inter-call durations within bursts. The likely reason is that text messages often follow a ``sending-response'' pattern while calls are self-contained so that a caller is answered usually in the same call to her query. In contrast, the ``sending-response'' nature of many text messages introduces an inherent delay that is in general broadly distributed as a result of priority-queuing processes \cite{Barabasi-2005-Nature,Vazquez-2005-PRL,Vazquez-Oliveira-Dezso-Goh-Kondor-Barabasi-2006-PRE,Maillart-Sornette-Frei-Duebendorfer-Saichev-2011-PRE}, circadian and weekly cycles \cite{Malmgren-Stouffer-Motter-Amaral-2008-PNAS,Malmgren-Stouffer-Campanharo-Amaral-2009-Science}, adaptive interests \cite{Han-Zhou-Wang-2008-NJP}, and so on.

In summary, we have found that individual call activities can be separated into independent bursts, and that the bursts and the call events within bursts are driven by two independent Poisson processes. We found the existence of a memoryless behavior of individual call activities, which are diagnosed by the exponential distribution of the number of call events in both high activity and low activity periods. We also documented exponential distributions of the time intervals between consecutive bursts and of the inter-call duration within bursts. We have proposed a minimal model to replicate the call activities for each individual, based on three ingredients: Poisson initiation of bursts, Poisson activation of calls within bursts and Markov-chain state shift between the two Poisson processes. We have integrated these three ingredients in a two-state first-order Markov-chain combining the two Poisson processes. By using the parameters estimated from real data, our model is able to reproduce the
empirical results very well. The Monte Carlo tests confirm that the validity of our model, which cannot be rejected at the significant level of 0.01 for 48.68\% of the total population of 73'339 cellphone users that we have analyzed. Our model is able to account for the diversity of individual behaviors within a simple universal structure. It avoids the caveat of model misidentification resulting from the use of population level data.

\ack
This work was partially supported by the NSFC Grant (11205057 and 11375064), the Ph.D. Programs Foundation of Ministry of Education of China Grant 20120074120028, China Scholarship Council (201406745014), and the Fundamental Research Funds for the Central Universities.

\section*{References}
\providecommand{\newblock}{}


\section{Introduction}
The atom-centered point charge (PC) model of molecular electrostatics has been a mainstay of biomolecular simulations for decades.\cite{Cox1981, Singh1984, Bayly1993, Hinsen1997, Simmonett2005,Dupradeau2010,Huang2013, Gotz2014, Dickson2014, Arnautova2015, Sprenger2015, Mukhopadhyay2015, Maier2015, Gabrieli2015}
While chemically intuitive and straightforward in technical implementation, this model does not provide a sufficiently detailed description of the anisotropic features of the molecular electrostatic potential (MEP), such as lone pairs, $\pi$-systems, and $\sigma$-holes, etc. which are mostly governed by higher-order multipole terms.\cite{Cardamone2014,Kramer2014}
These anisotropic effects, however, can be described within the PC approximation by moving beyond the atom-centered paradigm, i.e. by adding non-atom centered PCs/extended points.\cite{Karamertzanis2004, Cerutti2013, Gao2014, Devereux2014}
Although increasing the number of PCs per atom improves the quality of the electrostatic model, it also can exacerbate well-known ill-conditioning and redundancy problems\cite{Francl1995, Sigfridsson1997, Jakobsen2014} of the PC fitting procedures, leading to numerically unstable solutions.\cite{Tschampel2007, Dixon1997, Bayly1993}

These numerical instabilities are usually related to a large variation of the PC values for atoms in the interior of the molecule, so-called buried atom effect.\cite{Stouch1992,Stouch1993,Bayly1993, Hinsen1997}
The buried atom (usually methyl and methylene carbons) charges can dramatically change due to trivial changes in the PC fitting problem (the probe grid sampling, spatial orientation of the molecule, etc.), and/or have  inconsistent values across very similar molecules or even conformers of the same molecule.\cite{Breneman1990,Woods1989}
As the inclusion of non-atom centered PCs into the model produces even more buried centers, it should also increase the numerical instabilities of the PC fitting problem.

In fact, these numerical problems are rooted in the mathematical nature of the PC derivation---the least squares (LS) fitting to the reference MEP:\cite{Cox1981, Singh1984}
\begin{equation} \label{eq:chi2}
\chi^2 \left( \mathbf{q} \right) = | \mathbf{\Phi} - \mathbf{A}\mathbf{q} |^2 = |\mathbf{\Phi}|^2 + \mathbf{g}^\intercal \cdot \mathbf{q} + \mathbf{q}^\intercal \mathbf{H} \mathbf{q},
\end{equation}
\begin{equation}
\mathbf{g} = -2\mathbf{A}^\intercal \mathbf{\Phi},
\end{equation}
\begin{equation} \label{eq:hessian}
\mathbf{H} = \mathbf{A}^\intercal \mathbf{A},
\end{equation}
where the LS sum $\chi^2$ is the subject of minimization and the solution satisfies normal equations:\cite{lawson74}
\begin{equation} \label{eq:normal}
    \mathbf{A}^\intercal \mathbf{A} \mathbf{q} = \mathbf{A}^\intercal \mathbf{\Phi}.
\end{equation}
Here,
the elements of the LS matrix $\mathbf{A}$ correspond to the inverse distance $1/r_{ij}$
between the PC $i$ and the grid point $j$;
$\mathbf{\Phi}$ is $T$-dimensional vector of the reference values of MEP;
$\mathbf{q}$ is $N$-dimensional vector of the PC values;
$\mathbf{g}$ is the gradient of the function $\chi^2$ at the origin ($\mathbf{q}=0$);
$\mathbf{H}$ is the Hessian matrix of LS sum $\chi^2$.

While the ill-conditioning is common to many LS fitting problems,\cite{Hansen2010, Machta2013, Transtrum2015}
numerical difficulties associated with PC fitting are further compounded by commonly used total charge constraint using Lagrange multiplier.\cite{Besler1990,Francl1995, Cieplak1995,Sigfridsson1997}

One of the most widely used techniques to alleviate the numerical instabilities of PC fitting is to add artificial restraints to the PC values of the buried atoms.\cite{Bayly1993, Dupradeau2010, Zeng2013,Burger2013} Although this method can be extended to models with off-center PCs/extended points, one may wonder if it would be possible to overcome these difficulties in a more elegant way, based on better physical understanding of the problem.

For instance, an important insight can be gleaned from the eigendecomposition of the LS sum Hessian matrix (eq. \ref{eq:hessian}):
\begin{equation} \label{eq:kappa}
\mathbf{H} \mathbf{u}_i = \kappa_i \mathbf{u}_i.
\end{equation}
Indeed, the ill-conditioned nature of the LS matrix $\mathbf{A}$ can be related to the significant differences in the eigenvalues $\kappa_i$, i.e. the LS sum curvatures along the directions defined by the eigenvectors $\mathbf{u}_i$.\cite{Ivanov2015,Hansen2013} Because of the 2--3 order of magnitude variation of the $\kappa_i$ values, different sets of PCs can produce essentially the same MEP, as these solutions have the same positions along the high-curvature directions, although the positions along the low-curvature directions could be quite different.\cite{Ivanov2015} Importantly, the eigenvectors with the largest curvatures usually correspond to the total charge and dipole moment components of the molecule, while the lower-curvature eigenvectors do not seem to be associated with particular multipole moments.\cite{Ivanov2015, Laio2002}

However, the exact physical origin of the correspondence between the large curvature eigenvectors and the first terms of the multipole expansion is unclear, along with the nature of the low-curvature eigenvectors.
Particularly, it is not clear if the presence of the low-curvature modes of the $\mathbf{H}$ matrix and thus the ill-conditioning of the LS problem is solely because of the nature of the PC fitting problem, or due to some numerical factors, e.g. an incomplete sampling of the reference MEP grid.

To address these questions, we here revisit the PC fitting problem from the first principles. While the atom-centered PC model traces back to the intuitive chemical concept of the atomic charge, we consider a general PC model as a case of the inverse problem, where one seeks to recover the source charge distribution from its effect, i.e. electrostatic potential distribution. Based on the properties of the Coulomb law, we construct a best-case electrostatic model for which the inverse problem can be solved exactly, both in the continuous case, as well as in the case of a discrete (non-atom centered) PC approximation.

Using this model, we investigate the nature of the eigenvectors $\mathbf{u}_i$ and their eigenvalues $\kappa_i$, and dissect the factors responsible for the ill-conditioning of the LS fitting problem, and discuss how these insights can be used to improve and simplify the existing PC derivation procedures.

\section{Point Charge Fitting as an Inverse Problem}
A problem where given an effect (in this case the MEP $\Phi$) defined in the region $V_\Phi$, its cause (a charge distribution $\rho$) defined in the region $V_\rho$ needs to be determined belongs to a general class of inverse problems and can be described by the Fredholm integral equation of the first kind:\cite{Groetsch2007}
\begin{equation}\label{eq:integral}
    \int_{V_\rho} k(\mathbf{r}, \mathbf{r'})  \rho \left( \mathbf{r}'  \right) d \mathbf{r'} = \Phi \left( \mathbf{r} \right),
\end{equation}
where kernel $k(\mathbf{r}, \mathbf{r'})$
specifies the evolution of the cause $\rho(\mathbf{r}')$ into the effect  $\Phi(\mathbf{r})$, that in this case corresponds to the Coulomb law:
\begin{equation} \label{eq:kernel}
    k(\mathbf{r}, \mathbf{r}') = \frac{1}{|\mathbf{r}-\mathbf{r'}|}.
\end{equation}

The integral equation can also be represented as an operator equation:
\begin{equation} \label{eq:operator}
    K\rho=\Phi,
\end{equation}
where $K:U \rightarrow V$ is a linear operator defined on space $U=\mathrm{range}(K^\ast) \in L^2$ of square integrable functions,
and takes values in space $V=\mathrm{range}(K) \in L^2$;
$K^\ast: V \rightarrow U$ is adjoint of $K$. This equation can be solved exactly if and only if $\Phi \in V$.
However, in general it is not the case, so a function $\rho$ that minimizes the residual norm $| \Phi-K\rho |$
is considered as the LS solution and thus satisfies the normal equation:\cite{Groetsch2007, Roman2008}
\begin{equation}
    K^\ast K \rho = K^\ast \Phi.
\end{equation}

This LS solution can be obtained as the linear combination of the basis vectors $u_i \in U$: \cite{Groetsch2007}
\begin{equation} \label{eq:inverse}
  \rho = K^\dagger \Phi=\sum_{i=1}^\infty \frac{\langle\Phi,v_i\rangle}{\mu_i}u_i,
\end{equation}
where $K^\dagger$ is the Moore-Penrose inverse, $\mu_i$ is a singular value, $v_i$ and $u_i$ are left and right singular vectors, respectively and the inner product $\langle\Phi,v_i\rangle$ is defined as
\begin{equation}
	\langle\Phi,v_i\rangle = \int_{V_\Phi} \Phi(\mathbf{r}) v_i(\mathbf{r}) d \mathbf{r}
\end{equation}

The orthogonal bases $\{u_i\}_{i=1}^\infty$ and $\{v_i\}_{i=1}^\infty$ also form the eigenbases of $K^\ast K$ and $K K^\ast$ with eigenvalues $\mu_i^2$:
\begin{equation} \label{eq:K*K}
K^\ast K u_i = \mu_i^2 u_i,
\end{equation}
\begin{equation} \label{eq:KK*}
K K^\ast v_i = \mu_i^2 v_i.
\end{equation}

To obtain a numerical solution to the integral equation (eq. \ref{eq:integral}), the regions over which the MEP and charge distribution are defined are sampled using a numerical quadrature. Given $N$ quadrature nodes over the charge distribution and $T$ nodes over the MEP region the integral equation is transformed into a system of $T$ linear equations:
\begin{equation} \label{eq:linear}
    \mathbf{K} \mathbf{q} = \mathbf{\Phi},
\end{equation}
where the $T\times N$ matrix $\mathbf{K}$ is identical to the LS matrix $\mathbf{A}$ from eq. \ref{eq:chi2} and contains the kernel elements $k_{ij}$, as this matrix originates from the kernel $k(\mathbf{r}, \mathbf{r'})$ in the integral equation (eq. \ref{eq:integral}). It will be further referred to as $\mathbf{K}$ in order to highlight its mathematical origin.

Then, the PC value at the node $i$ is the product of the charge density $\rho_i$ and the quadrature weight $w_i$:
\begin{equation}
	q_i = \rho_i w_i
\end{equation}

Since the number of the reference values $T$ is usually larger than the number of the unknown PC values $N$, the system of linear equations is overdetermined. Then, a solution that minimizes the LS sum $\chi^2 \left( \mathbf{q} \right)$ (eq. \ref{eq:chi2}) and satisfies normal equations (eq. \ref{eq:normal}) is considered as the numerical solution to the integral equation (eq. \ref{eq:integral}). This solution can be obtained using singular value decomposition (SVD) of matrix $\mathbf{K}$: \cite{lawson74, Hansen2013, Roman2008}
\begin{equation} \label{eq:svd}
  \mathbf{q} = K^\dagger \mathbf{\Phi}=\sum_{i=1}^{r} \frac{ \mathbf{\Phi} \cdot \mathbf{v}_i } {\mu_i} \mathbf{u}_i,
\end{equation}
where  $K^\dagger$ is the Moore-Penrose pseudoinverse;
$\mu_i$ are singular values of matrix $\mathbf{K}$;
vectors $\mathbf{v}_i$ and $\mathbf{u}_i$ are left and right singular vectors.
If the rank $r$ of matrix $\mathbf{K}$ is less than the dimension of $\mathbf{q}$ ($r<N$), then the matrix $\mathbf{K}$ is rank deficient.

Similarly to the continuous case (eqs. \ref{eq:K*K}-\ref{eq:KK*}), the orthogonal bases $\{v_i\}_{i=0}^r$ and $\{u_i\}_{i=0}^r$ form eigenbases for $\mathbf{K} \mathbf{K}^\intercal$ and $\mathbf{K}^\intercal \mathbf{K}$:
\begin{equation}
    \mathbf{K} \mathbf{K}^\intercal \mathbf{v}_i = \mu_i^2 \mathbf{v}_i,
\end{equation}
\begin{equation} \label{eq:curvature}
    \mathbf{K}^\intercal \mathbf{K} \mathbf{u}_i = \mu_i^2 \mathbf{u}_i,
\end{equation}
where $\mathbf{K}^\intercal \mathbf{K}$ is also a Hessian matrix (eq. \ref{eq:kappa})
and $\mu_i^2$ is identical to its eigenvalue $\kappa_i$, which is the $\chi^2$ curvature along the direction $\mathbf{u}_i$:\cite{Brown2003}
\begin{equation}
	 \mu_i^2 =\kappa_i
\end{equation}

In many LS problems, PC fitting included, the singular values vary in a wide range, revealing the underlying ill-conditioning.\cite{Francl1995, Sigfridsson1997, Hansen2010, Machta2013}
As a singular value $\mu_i$ is a denominator in the LS solution (eq. \ref{eq:svd}), the smaller the singular value, the larger the effect of the corresponding singular vector $\mathbf{u}_i$ on the LS solution. Thus, even small variations along $\mathbf{u}_i$ with small singular value lead to a significant variations of the LS solution, although these variations do not lead to significant change in the quality of the fit $\chi^2$.\cite{Ivanov2015}
To understand the origins of the ill-conditioning in PC fitting, we next consider a system for which the inverse electrostatic problem can be analytically solved.

\begin{figure}[H]
\includegraphics{fig1
\caption{\label{fig:twospheres}
Schematic representations of the probe $S_R$ and charged $S_a$ spheres in the continuous (A) and discrete (B) forms.
Operators $K$ (eq. \ref{eq:K}) and matrix $\widetilde{\mathbf{K}}$
are represented schematically.}
\end{figure}

\section{The Two-Sphere Model}

The Coulomb kernel (eq. \ref{eq:kernel}) can be conveniently expanded in terms of spherical harmonics so the source $\mathbf{r}'$ and the observation $\mathbf{r}$ coordinates are separated but share the same origin:\cite{Stone,Jackson1999}
\begin{equation} \label{eq:expansion}
\begin{gathered}
k(\mathbf{r}, \mathbf{r}') = \frac{1}{|\mathbf{r}-\mathbf{r'}|} \\
= \sum_{l=0}^\infty \sum_{m=-l}^l \frac{4\pi}{2l+1} \frac{r_<^l}{r_>^{l+1}} Y_{lm} (\widehat{ \mathbf{r}'}) Y_{lm} (\widehat{ \mathbf{r} }),
\end{gathered}
\end{equation}
where $\widehat{ \mathbf{r}} = \mathbf{r}/r$ denotes the unit vector
defined by the polar $\varphi$ and azimuthal $\theta$ angles;
$r_<$ is the smaller and $r_>$ is the larger of $r$ and $r'$;
$Y_{lm}$ are orthogonal real-value spherical harmonics:\cite{f1}
\begin{equation} \label{eq:ortho}
    \int_S Y_{lm}(\widehat{\mathbf{r}}) Y_{l'm'}(\widehat{\mathbf{r}}) d \Omega = \delta_{ll'} \delta_{mm'},
\end{equation}
where $d \Omega$ is the differential of the solid angle.

Then, in the region beyond the divergence sphere where the charge density vanishes,
the MEP can be expanded in a multipole series:\cite{Jackson1999, Stone}
\begin{equation} \label{eq:mexp}
    \Phi(\mathbf{r}) = \sum_{l=0}^\infty \sum_{m=-l}^l \sqrt{ \frac{4\pi}{2l+1}  } r^{-l-1}  Q_{lm}^{mol} Y_{lm} (\widehat{ \mathbf{r}  }),
\end{equation}
where a molecular multipole moment $Q_{lm}^{mol}$ is given by
\begin{equation} \label{eq:Qlm}
    Q_{lm}^{mol} = \sqrt{ \frac{4\pi}{2l+1} } \int r^l \rho(\mathbf{r}) Y_{lm} (\widehat{\mathbf{r}}) d^3r.
\end{equation}

The form of the kernel expansion (eq.\ref{eq:expansion}) suggests that if the radii $r=R$ and $r'=a$ are fixed, the kernel $k(\mathbf{R}, \mathbf{a})$ can uniquely map a charge density over a spherical surface $S_a$ to the corresponding potential $\Phi(\mathbf{R})$ on a sphere $S_R$ and vice versa.
Thus, for a probe sphere $S_R$ with the radius $R$ greater than the radius of divergence sphere, the MEP can be reproduced exactly by a sphere $S_a$ with surface charge density $\sigma(\mathbf{a})$ such that the multipole moments of the sphere $Q_{lm}^{S_a}$ are equivalent to the multipole moments of the molecule $Q_{lm}^{mol}$:
\begin{equation} \label{eq:equiv}
    Q_{lm}^{S_a} \equiv Q_{lm}^{mol},
\end{equation}
where multipole moments of the sphere are:
\begin{equation}\label{eq:Qlm-sphere}
    Q_{lm}^{S_a} = \sqrt{ \frac{4\pi}{2l+1} } a^l \int_{S_a} \sigma(\mathbf{a}) Y_{lm}(\widehat{\mathbf{a}}) d\Omega.
\end{equation}

In this case, the original integral eq. \ref{eq:integral} is transformed into a surface integral equation:
\begin{equation} \label{eq:surface}
  \int_{S_a} k(\mathbf{R}, \mathbf{a}) \sigma(\mathbf{a}) d\Omega = \Phi(\mathbf{R}),
\end{equation}
or, equivalently, in an operator form
\begin{equation}
    K \sigma = \Phi,
\end{equation}
where $K: L^2(S_a) \rightarrow L^2(S_R)$ is a compact infinite-rank operator (Fig. \ref{fig:twospheres}A):
\begin{equation} \label{eq:K}
    K \sigma = \sum_{l=0}^\infty \sum_{m=-l}^l \mu_l \langle \sigma, Y_{lm}^{S_a} \rangle Y_{lm}^{S_R},
\end{equation}
where subscripts $S_a$ and $S_R$ denote the spheres, on which the corresponding spherical harmonics are defined;
the projection $\langle \sigma, Y_{lm}^{S_a}  \rangle$ is the inner product on the $L^2(S_a)$ space:
\begin{equation} \label{eq:sigmaY}
    \langle \sigma, Y_{lm}^{S_a}  \rangle = \int_{S_a} \sigma(\mathbf{a}) Y_{lm}^{S_a}(\widehat{\mathbf{a}}) d \Omega
\end{equation}
and for each degree $l$  there is a singular value $\mu_l$ in the form of the distance-dependent factor from the MEP expansion (eq. \ref{eq:expansion}):
\begin{equation} \label{eq:mu}
\mu_l = \frac{4\pi}{2l+1} \frac{a^l}{R^{l+1}}.
\end{equation}

Accordingly, the spherical harmonics $Y_{lm}^{S_R}$ and $Y_{lm}^{S_a}$ are left and right singular vectors and thus the eigenfunctions of the operators $K^\ast K$ and $K K^\ast$, while the squares of the singular values $\mu_l$ are their eigenvalues (eqs. \ref{eq:K*K}-\ref{eq:KK*}).
Since the singular values $\mu_l$ and spherical harmonics $Y_{lm}^{S_a}$ and $Y_{lm}^{S_R}$ form a singular system of the operator $K$, the solution to integral equation (eq. \ref{eq:surface}) can be expressed as:
\begin{equation} \label{eq:K+}
    \sigma = K^\dagger \Phi = \sum_{l=0}^\infty \sum_{m=-l}^l \frac{\langle \Phi, Y_{lm}^{S_R}  \rangle}{\mu_l} Y_{lm}^{S_a}.
\end{equation}

According to the multipole expansion (eq. \ref{eq:mexp}), the inner  product $\langle \Phi, Y_{lm}^{S_R} \rangle$ depends on the radius $R$ of the probe sphere and the multipole moments of the molecule:
\begin{equation} \label{eq:phiY}
    \langle \Phi, Y_{lm}^{S_R}  \rangle = \int_{S_R} \Phi(\mathbf{R}) Y_{lm}^{S_R}(\widehat{\mathbf{R}}) d \Omega =
    \sqrt{\frac{4\pi}{2l+1}} \frac{1}{R^{l+1}}Q_{lm}^{mol}.
\end{equation}

The dependence on the radius $R$ cancels out, so the charge density depends only on the radius $a$ of the sphere $S_a$ and the molecular multipole moments:
\begin{equation} \label{eq:sigma}
    \sigma(\mathbf{a}) = \sum_{l=0}^\infty \sum_{m=-l}^{l} \sqrt{ \frac{2l+1}{4\pi} } a^{-l}
    Y_{lm}^{S_a}(\widehat{\mathbf{a}})Q_{lm}^{mol},
\end{equation}
and the charged sphere $S_a$ exactly reproduces the MEP $\Phi(\mathbf{R})$.

\section{Analytical Lebedev Grid Point Charge Model}

We can construct an approximate discrete analog of the two-sphere model (eqs. \ref{eq:surface}-\ref{eq:sigma}, Fig. \ref{fig:twospheres}) using a quadrature that exactly integrates spherical harmonics $Y_{lm}$ over a sphere up to a given $l$ (eqs. \ref{eq:sigmaY} and \ref{eq:phiY}), e.g. the widely used\cite{Leang2012,Steinmann2012,Parrish2013} Lebedev quadrature,\cite{Lebedev1999} that defines $N$ quadrature nodes (Table I in the supplementary material\cite{SI}) with predetermined angular coordinates $\theta_i$, $\varphi_i$, and integration weights $w_i$:
\begin{equation}
    \int_S Y_{lm}(\theta, \varphi) d \Omega = \sum_i^N Y_{lm}(\theta_i, \varphi_i) w_i.
\end{equation}
Then, given the surface charge density $\sigma_i$ the corresponding point charge is:
\begin{equation}
	q_i = \sigma_i w_i.
\end{equation}

Due to the orthogonality of the spherical harmonics $Y_{lm}$ (eq. \ref{eq:ortho}), the $N$-node Lebedev quadrature that exactly integrates spherical harmonics over the sphere $S_a$ up to $l=2n$
\begin{equation} \label{eq:ortho2}
    \sum_i^N Y_{lm}^{S_a}(\theta_i, \varphi_i)
    Y_{l'm'}^{S_a}(\theta_i, \varphi_i) w_i^{S_a} =
    \widetilde{\mathbf{Y}}_{lm}^{S_a} \cdot \widetilde{\mathbf{Y}}_{l'm'}^{S_a} =
    \delta_{ll'} \delta_{mm'},
\end{equation}
defines an orthonormal basis:
\begin{equation} \label{eq:Ybasis}
    \widetilde{\mathbf{Y}}_{S_a} = \{\ \widetilde{\mathbf{Y}}_{lm}^{S_a},\ -l\leq m \leq l\ \}_{l=0}^n,
\end{equation}
where the $\widetilde{\mathbf{Y}}_{lm}^{S_a}$ vectors have $d_n = (n+1)^2$ elements defined as:
\begin{equation}
    \widetilde{Y}_{lmi}^{S_a} = Y_{lm}(\theta_i, \varphi_i) \sqrt{w^{S_a}_i}.
\end{equation}
Similarly, the probe sphere $S_R$ can be represented by a $T$-node Lebedev grid that integrates spherical harmonics up to $l=2t$ and  defines an orthogonal basis $\widetilde{\mathbf{Y}}_{S_R}$ of dimension $d_t = (t+1)^2$.

In this discrete representation, the operator $K$ (eq. \ref{eq:K}) then becomes a $T \times N$ matrix $\widetilde{\mathbf{K}}$:\cite{Ahrens2009}
\begin{equation} \label{eq:Ktilde}
    \widetilde{\mathbf{K}} \ \widetilde{\pmb{\sigma}}\ = \widetilde{\mathbf{\Phi}},
\end{equation}
where the elements of $\widetilde{\mathbf{K}}$, $\widetilde{\pmb{\sigma}}$, and $\widetilde{\mathbf{\Phi}}$ are:
\begin{gather}
    \widetilde{K}_{ij} = \sqrt{w_i^{S_a} w_j^{S_R}} / r_{ij},\\
    \widetilde{\sigma}_i = \sigma_i \sqrt{w_i^{S_a}},\ \
    \widetilde{\Phi}_j = \Phi_j \sqrt{w_j^{S_R}}.
\end{gather}

Since usually the probe grid has more points than the source grid, i.e. $T>N$, the matrix equation (eq. \ref{eq:Ktilde}) is a LS problem (eq. \ref{eq:chi2}) that can be solved using SVD of the matrix $\widetilde{\mathbf{K}}$ (eq. \ref{eq:svd}), giving a discrete analog of eq. \ref{eq:sigma}:
\begin{equation} \label{eq:svd2}
    \widetilde{\pmb{\sigma}} = \sum_{l=0}^n \sum_{m=-l}^l
    \frac{\widetilde{\mathbf{\Phi}} \cdot \widetilde{\mathbf{Y}}_{lm}^{S_R}} {\mu_l} \widetilde{\mathbf{Y}}_{lm}^{S_a},
\end{equation}
where
$\widetilde{\mathbf{Y}}_{lm}^{S_R}$ and $\widetilde{\mathbf{Y}}_{lm}^{S_a}$ are left and right singular vectors, and the corresponding singular values $\mu_l$ are the same as in the continuous case (eq. \ref{eq:mu}).

Since we use the Lebedev quadrature, the dot product $\widetilde{\mathbf{\Phi}} \cdot \widetilde{\mathbf{Y}}_{lm}^{S_R}$ corresponds to exact numerical integration and gives a result identical with the continuous case (eq. \ref{eq:phiY}):
\begin{equation}
    \widetilde{\mathbf{\Phi}} \cdot \widetilde{\mathbf{Y}}_{lm}^{S_R} = \sum_{j=0}^T \Phi_j Y_{lmj}^{S_R} w_j =
    \sqrt{\frac{4\pi}{2l+1}} \frac{1}{R^{l+1}}Q_{lm}^{mol},
\end{equation}
so the solution to eq. \ref{eq:Ktilde} depends only on the radius $a$ and the multipole moments $Q_{lm}^{mol}$:
\begin{equation} \label{eq:analytical_sigma}
\widetilde{\pmb{\sigma}} =
\sum_{l=0}^{n} \sum_{m=-l}^l \sqrt{ \frac{2l+1}{4\pi} } a^{-l} Q_{lm}^{mol} \widetilde{\mathbf{Y}}_{lm}^{S_a}.
\end{equation}

The corresponding PC values $q_j$ can be obtained using the quadrature weights $w_j^{S_a}$:
\begin{equation}
  q_i=\sigma_i w^{S_a}_i = \widetilde{\sigma}_i \sqrt{w_i^{S_a}},
\end{equation}
or, in a vector form:
\begin{equation} \label{eq:analytical}
\mathbf{q} =
\sum_{l=0}^{n} \sum_{m=-l}^l \sqrt{ \frac{2l+1}{4\pi} } a^{-l} Q_{lm}^{mol} \mathbf{Y}_{lm}^{S_a} \odot \mathbf{w}^{S_a},
\end{equation}
where $\mathbf{w}^{S_a}$ is the vector of the quadrature weights for the sphere $S_a$. Therefore, we can use Lebedev grid that shares the origin with a molecule to construct an analytical PC model that exactly reproduces molecular multipole values up to the degree $n$.

From this model, we can see that the ill-conditioning of the PC fitting due to the decay of the singular values is intrinsic to the inverse electrostatic problem, as the singular values $\mu_l$ decrease with increasing $l$ (eq. \ref{eq:mu}).
Indeed, the higher the multipole moment, the smaller its contribution to the overall electrostatic potential. Also, this contribution  gets smaller as we move the probe further away from the source, and the singular values get smaller with the increasing radius of the probe sphere $R$, or decreasing radius of the source sphere $a$.

The ill-conditioning problems become even more severe as we switch from modeling the MEP using the Lebedev quadrature, which is the best suited to reproduce the molecular multipoles, to an irregular atom-centered quadrature, as shown on a numerical example below.

\begin{figure}[h]
\includegraphics{fig2
\caption{\label{fig:models}
Cross-section representations of the quadratures used for two-sphere model (A) ($n=1$, $N=6$ and $t=11$, $T=194$ for spheres $S_a$ and $S_R$, respectively) as compared with the traditional atom-centered model (B). Green circles correspond to the point charges; blue circles correspond to the reference grid points}
\end{figure}

\section{Lebedev Grid vs. Atom-Centered Model: a Numerical Example}

First, we consider an electrostatic PC model of a methanol molecule with PCs placed at the nodes of the Lebedev quadrature over the sphere $S_a$ ($a=2$ au) (Fig. \ref{fig:models}A). In this case, the PC values can be obtained analytically from the reference multipole moments (eq. \ref{eq:analytical}) or by numerical fitting to the reference MEP over the probe sphere $S_R$ ($R=8$ au, $t=11$, $T=194$):
\begin{equation} \label{eq:numerical}
    \widetilde{\bm{\sigma}} =\sum_{i=1}^{r} \frac{ \mathbf{\widetilde{\Phi}} \cdot \mathbf{\widetilde{v}}_i } {\mu_i}
    \mathbf{\widetilde{u}}_i,
\end{equation}
where the PC value can be found as $q_j = \widetilde{\sigma}_j \sqrt{w_j^{S_a}}$ and
the maximum rank $r$ is the number $N$ of quadrature nodes/PCs over the sphere $S_a$.
The quality of the fit is measured using the root mean square deviation (RMSD) calculated over the $T$ nodes of the probe grid:

\begin{equation}
\mathrm{RMSD} = \sqrt{ \frac{\chi^2}{T}  }.
\end{equation}

Naturally, the analytical PC values from eq. \ref{eq:analytical} exactly reproduce the molecular multipole moments up to the degree $n$ defined by the quadrature (Table \ref{tab:multipoles}).
For each degree $l$ there are $2l+1$ values of order $m$, so overall $(n+1)^2$ multipole moments are reproduced, which matches the dimension $d_n$ of the corresponding basis $\widetilde{\mathbf{Y}}_{S_a}$ (eq. \ref{eq:Ybasis}).
As the dimension $d_n$ increases, more multipole moments are reproduced and the RMSD rapidly approaches zero (Fig. 4 in the supplementary material\cite{SI}).

Since the dimension $d_n$ does not match the number of quadrature nodes $N$ (Table I in the supplementary material\cite{SI}),\cite{Lebedev1976, Sloan1995} we can obtain numerical solutions with eq. \ref{eq:numerical} that are equivalent to the analytical results (eq. \ref{eq:analytical}) by setting the rank $r$ to the dimension of the grid, $d_n=(n+1)^2$ (Table \ref{tab:multipoles}).\cite{f2}

\begin{figure}[t]
\includegraphics{fig3
\caption{\label{fig:singulars}
Normalized singular values $\mu_i/\mu_1$ obtained using the exact analytical expression eq. \ref{eq:mu} (green circles) as compared with the numerical values obtained from SVD of the LS matrix for the two-sphere model (red stars) and for atom-centered model (black circles).
Lebedev quadratures with $n=1$, $N=6$ and $t=11$, $T=194$ were used for the charged $S_a$ ($a=2$ au) and probe $S_R$ ($R=8$ au) spheres, respectively.}
\end{figure}

As the first $d_n$ multipole moments $Q_{lm}^{mol}$ are reproduced by the PC model, the first $d_n$ numerical singular values $\mu_i$ exactly match the radius-dependent part (eq. \ref{eq:mu}) from the inverse distance expansion (Fig. \ref{fig:singulars}), and the corresponding right singular vectors $\mathbf{\widetilde{u}}_i$ match the basis $\widetilde{\mathbf{Y}}_{S_a}$ (Fig. \ref{fig:bases}):
\begin{equation}
    \{ \mathbf{\widetilde{u}}_i \}_{i=1}^{d_n} = \{ \widetilde{\mathbf{Y}}_{lm}^{S_a},\ -l\leq m \leq l  \}_{l=0}^{n}.
\end{equation}

If we do not restrict the rank $r$  to the dimension of the grid $d_n$, numerical SVD of the LS matrix $\widetilde{\mathbf{K}}$ (eq. \ref{eq:numerical}) produces $N$ singular vectors/values. While this slightly improves the RMSD (Table \ref{tab:multipoles}), the additional $N-d_n$ singular vectors cannot be described analytically (Fig. \ref{fig:bases}), as they go beyond the dimension $d_n$ of the corresponding basis $\widetilde{\mathbf{Y}}_{S_a}$.
However, in the fortuitous case of the quadrature with $n=1$ and $N=6$, the remaining $6-4=2$ vectors resemble the basis vectors
$\widetilde{\mathbf{Y}}_{2-2}$ and $\widetilde{\mathbf{Y}}_{2-1}$, so the corresponding quadrupole moments $Q_{2-2}$ and $Q_{2-1}$ are accurately reproduced, although the exact numerical integration of the spherical harmonics $Y_{2-2}$ and $Y_{2-1}$ is not provided by the 6-node Lebedev grid.

Now, we can use the insights from the best-case scenario spherical PC model based on the Lebedev quadrature (Fig. \ref{fig:models}A) to understand the traditional PC fitting problem with atom-centered charges and the probe grid that follows the solvent-accessible surface (vdW grid, Fig. \ref{fig:models}B). From the point of view of the inverse electrostatic model, the atom-centered PC fitting  corresponds to a numerical solution using an irregular and suboptimal integration grid to represent the source charge distribution. This problem can be treated by SVD of the LS matrix $\mathbf{K}$:
\begin{equation} \label{eq:numerical2}
    \mathbf{q} =\sum_{i=1}^{r} \frac{ \mathbf{\Phi} \cdot \mathbf{v}_i } {\mu_i} \mathbf{u}_i,
\end{equation}
where the maximum value of rank $r$ is the number of atoms in the molecule, i.e. $r=6$ in the case of methanol.

We can see that even in this case the singular vector $\mathbf{u}_1$ with the largest singular value $\mu_1$ corresponds the total charge (Fig. \ref{fig:bases}), which is reproduced with only a slight slight numerical deviation ($<0.01$), a consequence of the molecular charge density spillover beyond the solvent-accessible surface defining the vdW grid.\cite{Ivanov2015,Laio2002}

Although the other singular vectors do not exactly match the corresponding spherical harmonics, the $\mathbf{u}_2$--$\mathbf{u}_4$ vectors can be roughly related to the three components of the dipole moment (Fig. \ref{fig:bases}), and the corresponding singular values are commensurate with the  singular value $\mu_l$ ($l=1$) obtained for the Lebedev grid model (Fig. \ref{fig:singulars}).
The remaining singular values $\mu_5$ and $\mu_6$ are significantly distorted from the singular value $\mu_l$ ($l=2$), so the components of the quadrupole moment are not reproduced as precisely as the the dipole moment components (Table \ref{tab:multipoles}).

Among all singular vectors $\{ \mathbf{u}_i \}_{i=1}^6$, the singular vector $\mathbf{u}_6$ with the lowest singular value $\mu_6$, which is 100 times smaller than $\mu_1$, is dominated by the contribution from the methyl carbon atom (Fig. \ref{fig:bases}).
Since such small singular values cause numerical instabilities of the LS solution, once can use a regularization technique such as truncated SVD (tSVD) that reduces the rank $r$ by removing the lowest-$\mu_i$ vector(s) from the SVD expansion.\cite{Hansen2013}
Removal of $\mathbf{u}_6$ that decreases the rank to $r=5$ leads to dramatic change in the methyl group charges---the carbon atom charge in particular, which drops from $0.22$ to $-0.06$. Yet, these changes lead only to marginal changes in the the multipole moment and RMSD values, a typical example of the buried atom effect (Tables \ref{tab:multipoles}, \ref{tab:atomic_charges}).
This suggests a natural way to impose a restraint on the buried atom charges without introducing a restraining function into the LS sum $\chi^2$, an addition that can negatively affect the electrostatic properties of the PC model.\cite{Sigfridsson1997, Verstraelen2014}

Further removal of the singular vectors $\mathbf{u}_5$ and $\mathbf{u}_6$ (i.e. $r=4$) leads to severe deterioration of the LS solution,
as the corresponding multipole moment strongly deviate from the reference values and the RMSD significantly increases (Tables \ref{tab:multipoles}, \ref{tab:atomic_charges}).
Thus, it appears that the tSVD approach should be applied only to the singular vectors that strongly depend on the buried atoms, an important point that will be discussed in detail elsewhere.

\begin{figure*}[!htbp]
    \includegraphics{fig4
\caption{\label{fig:bases}
The orthonormal bases of the right singular vectors: basis of spherical harmonics $\widetilde{\mathbf{Y}}_{S_a}$ (A), basis from the numerical SVD of the LS matrix in two-sphere PC model (B), and atom-centered model (C).
}
\end{figure*}

\section{The total-charge constraint revisited}

Commonly used PC fitting approaches also modify the LS sum (eq. \ref{eq:chi2}) by adding a Lagrange multiplier $\lambda$ in order to constrain the total charge to the correct value:\cite{Besler1990, Francl1995, Cieplak1995}
\begin{equation}
    \chi^2 \left( \mathbf{q} \right) = |\mathbf{\Phi} - \mathbf{K} \mathbf{q}|^2 + \lambda (\mathbf{1}^\intercal \cdot \mathbf{q} - Q_{0}),
\end{equation}
which increases the dimension of the Hessian matrix
$\mathbf{H} = \mathbf{K}^\intercal \mathbf{K}$ in the normal equation (eq. \ref{eq:normal}):
\begin{equation} \label{eq:normal2}
    \begin{bmatrix}
        \mathbf{H} & \mathbf{1} \\
        \mathbf{1}^\intercal & 0 \\
    \end{bmatrix}
    \begin{bmatrix}
        \mathbf{q} \\ \lambda\\
    \end{bmatrix}
    =
    \begin{bmatrix}
        \mathbf{K}^\intercal \mathbf{\Phi} \\ Q_0 \\
    \end{bmatrix},
\end{equation}
where $\mathbf{1}$ is an all-ones column-vector.

However, as we have seen, both in the case of the idealized Lebedev grid and the less-than-ideal atom-centered PC models, the Hessian eigenvector with the largest curvature corresponds to the total charge (Fig. \ref{fig:bases} and also Ref. \cite{Ivanov2015}).
Thus, in the case of the two-sphere PC fitting, the total charge is reproduced exactly ($Q_0<10^{-5}$), while in the atom-centered PC model the total charge only slightly deviates from the exact value due the close proximity of the vdW grid and slight distortion of the total-charge vector $\mathbf{u}_1$ from its analytical analog $\widetilde{\mathbf{Y}}_0$ ($Q_0=0.003$ for methanol, Table \ref{tab:multipoles}).

\begin{table*}[!htbp]
\caption
{\label{tab:multipoles}
Effect of the rank $r$ (eq. \ref{eq:numerical}), degree $n$ (eq. \ref{eq:analytical}) and type of the charge constraint on the methanol multipole moments and the RMSD (kcal/mol) within the Lebedev grid PC model ($a=2$ au, $n=1,2$ and $N=6,14$) with probe sphere $S_R$ ($R=8$ au, $T=194$) and atom-centered PC model with vdW-type grid.
}

\begin{ruledtabular}
\begin{tabular}{ccccccccccccc}
\multirow{1}{*}{PC model, probe grid}
& Details & $Q_{0}$ & $Q_{10}$ & $Q_{11}$ & $Q_{1-1}$ & $Q_{20}$ & $Q_{21}$ & $Q_{2-1}$ & $Q_{22}$ & $Q_{2-2}$ & RMSD \\
\hline
\multirow{3}{*}{$S_a(N=6)$, $S_R (T=194)$}
& eq. \ref{eq:analytical}, $n=1$& 0.000 & 0.000 & -0.325 & -0.565 &  0.000 & 0.000 & 0.000 &  0.000 & 0.000 & 1.762 \\
& SVD, $r=4$ & 0.000 & 0.000 & -0.323 & -0.562 &  0.000 & 0.000 & 0.000 &  0.000 & 0.000 & 1.762 \\
& SVD, $r=6$ & 0.000 & 0.000 & -0.323 & -0.562 & -0.881 & 0.000 & 0.000 & -0.497 & 0.000 & 1.668 \\
\hline
\multirow{3}{*}{$S_a(N=14)$, $S_R (T=194)$}
& eq. \ref{eq:analytical}, $n=2$& 0.000 & 0.000 & -0.325 & -0.565 & -0.887 & 0.000 & 0.000 & -0.503 & 2.882 & 0.559 \\
& SVD, $r=9$ & 0.000 & 0.000 & -0.325 & -0.566 & -0.881 & 0.000 & 0.000 & -0.497 & 2.865 & 0.559 \\
& SVD, $r=14$& 0.000 & 0.000 & -0.325 & -0.566 & -0.881 & 0.000 & 0.000 & -0.497 & 2.865 & 0.425 \\
\hline
\multirow{5}{*}{atom-centered, vdW}
& SVD, $r=6$ & 0.007  & 0.000 & -0.335 & -0.546 & -1.020 & 0.001 & 0.001 & 0.110 & 2.761 & 2.457   \\
& tSVD, $r=5$ & 0.009  & 0.000 & -0.338 & -0.543 & -1.060 & -0.002 & 0.000 & 0.219 & 2.688 & 2.625 \\
& tSVD, $r=4$ & 0.003 & 0.001 & -0.270 & -0.377 & 0.574 & -0.003 & -0.003 & 0.225 & -1.317 & 8.729 \\
& Lagrange, $Q_0 = 0$
        & 0.000 & 0.000 & -0.332 & -0.543 & -0.989 & 0.001 & 0.001 & 0.073 & 2.739 & 2.587 \\
& Elimination, $Q_0 = 0$
        & 0.000 & 0.000 & -0.332 & -0.543 & -0.989 & 0.001 & 0.001 & 0.073 & 2.739 & 2.587 \\
& SVD, $Q_0 = 0$
        & 0.000 & 0.000 & -0.331 & -0.544 & -1.010 & 0.001 & 0.001 & 0.097 & 2.756 & 2.597 \\
& Trivial, $Q_0 = 0$
        & 0.000 & 0.000 & -0.331 & -0.544 & -1.010 & 0.001 & 0.001 & 0.097 & 2.755 & 2.597 \\

\hline
\multicolumn{2}{c}{Reference}
& 0.000 & 0.000 & -0.325 & -0.565 & -0.887 & 0.000 & 0.000 & -0.503 & 2.882 & --   \\
\end{tabular}
\end{ruledtabular}
\end{table*}

\begin{table*}[!htbp]
\caption
{\label{tab:atomic_charges}
Effect of the numerical rank $r$ (SVD in eq. \ref{eq:numerical2})
and the total-charge constraint
on the values of atom-centered PCs of methanol and the RMSD (kcal/mol).
}
\begin{ruledtabular}
\begin{tabular}{ccccccccccccc}
& & $q_C$ & $q_{H_g}$ & $q_{H_t}$ & $q_O$ & $q_H$ & RMSD & \\
\hline
& SVD, $r=6$ &  0.215 & -0.018 &  0.048 & -0.592 &  0.371 & 2.457 \\
& tSVD, $r=5$ & -0.058 &  0.056 &  0.118 & -0.532 &  0.370 & 2.625 \\
& tSVD, $r=4$ &  0.007 &  0.089 & -0.101 & -0.070 & -0.010 & 8.729  \\
& Lagrange, $Q_0=0$
        &  0.276 & -0.035 &  0.030 & -0.603 &  0.367 & 2.587 \\
& Elimination, $Q_0=0$
        &  0.276 & -0.035 &  0.030 & -0.603 &  0.367 & 2.587 \\
& SVD, $Q_0=0$
        &  0.214 & -0.019 &  0.047 & -0.593 &  0.370 & 2.597 \\
& Trivial, $Q_0=0$
        & 0.214 &  -0.019 &  0.047 & -0.593 &  0.370 & 2.597 \\
\end{tabular}
\end{ruledtabular}
\end{table*}

Addition of the Lagrange multiplier leads to an extra eigenvector $\mathbf{u}_7$ that appears in the eigenbasis of the Hessian matrix (Table VIII in the supplementary material\cite{SI}).
The curvature along this vector is the smallest in the magnitude ($\kappa_7=-0.009$) and the vector itself primarily depends on the Lagrange multiplier $\lambda$, with only marginal contribution from the PC values.
At the same time, remaining eigenvectors $\{\mathbf{u}\}_{i=1}^6$ preserve the structure of the original eigenbasis, with negligible contribution from the Lagrange multiplier $\lambda$ (Table VIII in the supplementary material\cite{SI}).
Thus, application of the the total charge constraint in addition to already strong restraint (imposed by the eigenvector $\mathbf{u}_1$) appears to be redundant. Moreover, addition of the Lagrange multiplier aggravates the rank deficiency of already ill-conditioned LS problem.\cite{Francl1995, Sigfridsson1997}

Alternatively, the total charge can be constrained by incorporating condition on the proper total charge directly into the LS sum,\cite{Cox1981, Hinsen1997, Jakobsen2014} by eliminating one of the charges and setting it to:
\begin{equation}
	q_n = Q_0^{mol} - \sum_i^{N-1} q_i,
\end{equation}
where $n$ is the index of the eliminated charge. This reduces the dimension of the LS problem by one:
\begin{equation}
	\chi^2 \left( \mathbf{q} \right) = \sum_j^T \left[ \Phi_j - \frac{Q_0^{mol}}{r_{nj}} - \sum_i^{N-1} \left( \frac{1}{r_{ij}} - \frac{1}{r_{nj}} \right) q_i  \right]^2
\end{equation}
and modifies the elements of the Hessian matrix:
\begin{equation}
	H_{km} = \sum_j^{T} \left( \frac{1}{r_{kj}} - \frac{1}{r_{nj}} \right) \left( \frac{1}{r_{mj}} - \frac{1}{r_{nj}} \right).
\end{equation}

Although the solution obtained with this approach is numerically equivalent to the solution with Lagrange multiplier, regardless which atom has been eliminated (Elimination, $Q_0=0$ in Tables \ref{tab:multipoles} and \ref{tab:atomic_charges}), the structure of the right singular vectors becomes disrupted, (Fig. 8 in the supplementary material\cite{SI}) which prevents the application of the truncated SVD to improve the numerical stability of the solution.

Given that even for the atom-centered PC/vdW probe model the total charge value deviates only very slightly from the reference value, it should be possible to correct for this deviation without exacerbating the numerical instabilities of the LS problem, e.g. using the total charge vector $\mathbf{u}_1$. To do that, we convert the SVD solution (eq. \ref{eq:numerical2}) to a system of linear equations:
\begin{equation}\label{eq:svdQ0corr}
    \mathbf{q} = \sum_{i=0}^r \underbrace{\frac{\Phi \cdot \mathbf{v}_i} {\mu_i}}_{c_i} \mathbf{u}_i =
    \mathbf{U} \mathbf{c},
\end{equation}
\begin{equation} \label{ref:sys}
    \mathbf{U}^\intercal \mathbf{q} = \mathbf{c}.
\end{equation}
Then, we replace $\mathbf{u}_1$ in $\mathbf{U}^\intercal$ by an all-ones vector $\mathbf{1}$, and set the corresponding coefficient $c_1$ in $\mathbf{c}$ to the exact value of the molecular total charge $Q_0^{mol}$:
\begin{equation}\label{eq:svdQ0corr1}
    \mathbf{U}_{Q_0}^\intercal \mathbf{q} = \mathbf{c}_{Q_0},
\end{equation}
where
\begin{equation}
    \mathbf{U}_{Q_0}^\intercal = \begin{bmatrix} \mathbf{1} & \mathbf{u}_2 & \cdots & \mathbf{u}_N \end{bmatrix}^\intercal,
\end{equation}
\begin{equation}\label{eq:svdQ0corr3}
    \mathbf{c}_{Q_0} = \begin{bmatrix} Q_0^{mol} & c_2 & \cdots & c_N \end{bmatrix}^\intercal.
\end{equation}

This approach does not introduce any redundant constraints, preserves the electrostatic properties of the unconstrained solution, and results only in to minor changes in the PC values (SVD, $Q_0=0$ in Tables \ref{tab:multipoles} and \ref{tab:atomic_charges}) and is compatible with truncated SVD.
Also, the error in the total charge value is small enough and can be corrected by simply distributing the $Q_0$ error correction across the atomic charges; this trivial total charge correction gives result nearly identical to eq. \ref{eq:svdQ0corr1} (Trivial, $Q_0=0$ in Tables \ref{tab:multipoles} and \ref{tab:atomic_charges}).

\section{Conclusions}

To understand the origins of the ill-conditioning of the least-squares (LS) point charge (PC) fitting problem, we revisited the PC representation of the molecular electrostatic potential (MEP) from the first principles, as an example of the inverse problem.

Based on the properties of the Coulomb potential that can be expanded in terms of spherical harmonics, we introduce a model where the MEP of a molecule is exactly reproduced by a charged sphere that has the same multipole moments $Q_{lm}$ as the molecule. Using Lebedev quadrature this continuous model is converted into a discrete PC model, where the PC values are evaluated analytically from the multipole moments $Q_{lm}$ up to the maximum value determined by the quadrature.

In this context, the traditional atom-centered PC model can be viewed as an irregular numerical quadrature, poorly suited to reproduce the multipolar expansion of the MEP. As such, this quadrature only allows integration of the monopole and, approximately, dipole terms. The corresponding large-curvature---or `stiff'\cite{Machta2013,Transtrum2015}---Hessian eigenvectors $\mathbf{u}_i$ can still be related to the corresponding multipoles $Q_{lm}$. This explains previously observed correspondence between the highest-curvature Hessian eigenvectors and the total charge and the dipole moment components;\cite{Ivanov2015,Laio2002} this correspondence quickly breaks down for the higher multipole moments.

This consideration then reveals the origins of the ill-conditioning of the PC fitting due to the presence of low-curvature---or `sloppy'\cite{Machta2013,Transtrum2015}---vectors  $\mathbf{u}_i$. The intrinsic ill-conditioning arises even in the case of the ideal spherical model: since the higher-rank multipole moments $Q_{lm}$ have smaller contribution to the MEP, the singular values $\mu_l$  decay as $l$ increases. The ill-conditioning is further exacerbated in the numerical treatment of the Lebedev grid model because the number of PCs does not match the dimension of the basis formed by Lebedev quadrature. The remaining singular values/curvatures are even lower in magnitude and do not correspond to particular multipole moments $Q_{lm}$. The same rank-deficiency problems apply to the atom-centered PC grids. However, in that case most of the eigenvectors do not have a direct correspondence to the multipole moments, which leads to even wider spread-out of the singular values/curvatures.

These insights can suggest several ways to alleviate the ill-conditioning of the problem. For instance, the buried atom problem can be addressed by truncating the sloppy singular vectors with dominant contribution from these atom, instead of introducing additional restraining functions\cite{Bayly1993, Dupradeau2010, Zeng2013,Burger2013} that can negatively affect the overall electrostatic properties of the molecule.\cite{Verstraelen2014,Sigfridsson1997} Also, slight deviations of the total charge of the fitted PC solution can be fixed by adjusting the stiff total-charge vector $\mathbf{u}_1$ and the corresponding coordinate $Q_0^{mol}$, rather than introducing a Lagrange multiplier that increases the rank-deficiency of the Hessian matrix.\cite{Francl1995, Sigfridsson1997}

The results presented here can help further application of the PC model in biomolecular simulations. Although the force fields using point charges may not be as accurate as the force fields that explicitly include multipoles and/or polarization effects, the simplicity and computational efficiency of the PC model has ensured its continued survival.\cite{Cerutti2013} In fact, representation of multipoles using the Lebedev grid PC model can provide an alternative to the multipole moment expansion;\cite{Rogers2015} it also can be used to extend recently proposed Distributed Charge Model.\cite{Devereux2014,Gao2014}

\section{Computational details}
MEP and multipole moments were calculated at the B3LYP/aug-cc-pVDZ level\cite{Lee1988, Becke1993, Stephens1994, Dunning1989} as implemented in Q-Chem package.\cite{QCHEM4}
For atom-centered PC fitting the reference MEP was generated as the cubic grid with linear density 2.8 points/\AA,
followed by the removal of the points outside of 1.0-2.0 van der Waals radii range around each atom (vdW grid).
For the two-sphere PC model the Lebedev quadrature rules were used as implemented in PyQuante package.
\cite{pyquante, pyquante2}
Charge fitting procedures were implemented in the in-house developed \textit{fftoolbox} Python library.\cite{fftoolbox}
SVD was performed using \textit{numpy} library.\cite{numpy}
Spherical harmonics were accessed from \textit{scipy} library.\cite{scipy}

\begin{acknowledgments}
Q.K.T. is a recipient of the National Science Foundation (NSF) CAREER award CHE-1255641, and Marquette University Way-Klinger Young Scholar Award.
M.V.I. is a recipient of Bournique Memorial Fellowship (Marquette University).
This work used the high-performance computing cluster Père at Marquette University funded by NSF awards OCI-0923037 and CBET-0521602, and the Extreme Science and Engineering Discovery Environment (XSEDE) supported by NSF grant ACI-1053575.
\end{acknowledgments}