\section{Introduction}

Shared digital whiteboards are becoming increasingly popular in educational and workplace settings as a natural mechanism for collaboration and communication~\cite{jensen2018remediating,fang2019augmented,davidson2021drawn,perlich2018cooperative,10.1145/3059454.3059477,10.1145/1054972.1055046}. The sharing aspect offers tremendous scope for interaction and a richer session experience. Unfortunately, shared whiteboards sometimes present situations where malicious participants draw controversial content~\cite{gameemblem}. Such activities tend to impair the collective experience of participants. Therefore, it is important to have scalable mechanisms for efficiently identifying and responding to such activities.

The popular social sketching game of Pictionary\textsuperscript{\texttrademark}, which we employ as a use case in this paper, also presents scenarios involving atypical sketched content. Pictionary is a wonderful example of cooperative game play to achieve a shared goal in communication-restricted settings~\cite{10.1145/3313831.3376316,10.1145/1378704.1378719,10.1145/3173574.3173767}. The game consists of a time-limited episode involving two players - a Drawer and a Guesser. The Drawer is tasked with conveying a given target phrase to a counterpart Guesser by sketching on a whiteboard~\cite{fay2017deconstructing}. The larger the number of target phrases correctly identified and the earlier the phrases are identified from the drawn sketch, greater the number of points accrued for the participating players. The rules of Pictionary forbid the Drawer from writing text on the whiteboard. This is usually not an issue when players are physically co-located. In the anonymized, web-based version of the game, however, the Drawer may cheat by writing text related to the target word on the digitally shared whiteboard, thus violating the rules. Intervention is possible by physically monitoring game sessions. However, such manual intervention is impractical and not scalable to an online setting involving a large number of multiple concurrent game sessions. Providing user interface options for player-triggered flagging of rule violation is another possibility. But such mechanisms are not completely reliable since the Guesser benefits from the content written on the canvas and does not have real incentive to use the flagging mechanism.

Apart from malicious game play, atypical sketch content can also exist in non-malicious, benign scenarios. For instance, the Drawer may choose to draw arrows and other such icons to attract the Guesser's attention and provide indirect hints regarding the target word (see Fig.~\ref{fig:ano-ex}). Accurately localizing such activities can aid statistical learning approaches which associate sketch-based representations with corresponding target words~\cite{8509167}. Considering both malicious and benign scenarios, the broad requirement is for a framework which can respond to a variety of atypical whiteboard sketch content in a reliable, comprehensive and timely manner. To this end, we make the following contributions:

\begin{itemize}
    \item We introduce \textsc{AtyPict} - the first ever dataset of atypical whiteboard content.
    \item We introduce \textsc{DrawMon}, a distributed system for sketch content-based alert generation (Sec.~\ref{sec:drawmon}). We analyze sessions with \textsc{DrawMon} deployed for  Pictionary setting and demonstrate its effectiveness (Sec.~\ref{sec:drawmonuserstudy}).
\end{itemize}

Although presented in a Pictionary game context, our contributions serve as a design guide for developing  response frameworks involving shared and interactive whiteboards. For code, models and additional details, visit the project page \url{https://drawm0n.github.io}

\section{Related Work}

\noindent \textbf{Detecting and Flagging Anomalous Gameplay:}
Some approaches employ a diverse mix of techniques for detecting cheating in online games~\cite{10.1007/978-3-642-38577-3_53,10.1145/1852611.1852643}. Dinh et al.~\cite{7725645} use hand-crafted game features and unsupervised machine learning approaches for offline detection of anomalous behavior. In our work, we introduce automatic deep learning based detection and flagging of anomalous gameplay in Pictionary. However, our system is also designed to detect secondary non-anomalous canvas entities which can potentially aid statistical understanding of canvas contents.

\noindent \textbf{Sketch datasets:} Existing sketch datasets (e.g. TU-Berlin~\cite{eitz2012hdhso}, Sketchy~\cite{sketchy2016}, QuickDraw~\cite{jongejan2016quick}) have been created primarily in the context of sketch \textit{object} recognition problem -- assign a categorical label to a hand-drawn sketch. The category labels correspond to objects (nouns). Therefore, these datasets lack abstract sketches which tend to be drawn when words from other parts of speech (verbs, adjectives) are provided as targets. Existing datasets are also unnatural because they do not include canvas actions such as erase strokes or location emphasis. Also, no intermediate guess words are associated with sketched content. For a similar reason, these datasets do not contain atypical activities unlike the dataset we introduce. Sarvadevabhatla et al.~\cite{8509167} explore neural network based generation of human-like guesses, but for pre-drawn object sketches. However, they do not accommodate interactivity and non-sketch drawing canvas activities (e.g. erase, pointing emphasis). The Kondate dataset~\cite{6981047} contains on-line handwritten patterns of text, figures, tables, maps, diagrams etc. The OHFCD dataset~\cite{Awal2011FirstEO} pertains to online handwritten flowcharts.  Although challenging in their own way, these datasets are considerably more structured than our setting. Additionally, they share the sketch datasets' shortcoming of being too cleanly curated because actions such as erase are absent. As a unique aspect, our combination of a game setting and a time limit unleashes greater diversity and creativity, causing sketches in our dataset to be more spontaneous and less homogeneous compared to existing datasets.

\noindent \textbf{Detecting canvas items:} Recognizing atypical activities can be thought of as a stroke segmentation problem wherein each sketch stroke is labelled as either belonging to an atypical class or the default class (drawing). Stroke segmentation has been employed for labelling parts in object sketches either from stroke sequence information~\cite{yang2020sketchgcn,8784880,kaiyrbekov2019strokebased,qi2019sketchsegnet+} or within an image canvas~\cite{wang2020efficient,li2018fast}.
\begin{figure*}[t!]
  \centering
  \includegraphics[width=\textwidth]{anomaly_examples3.png}
  \caption{Some examples of atypical sketch content in Pictionary game sessions are shown as canvas screenshots. The content instances span \textcolor{red}{text}, \textcolor{cyan}{numbers}, \textcolor{ForestGreen}{question marks}, \textcolor{blue}{arrows}, \textcolor{Maroon}{circles} and other \textcolor{orange}{icons} (e.g. tick marks, addition symbol) categories - refer to Sec~\ref{sec:atypical-data} for details.}
  \label{fig:ano-ex}
\end{figure*}

Recognizing atypical sketch content can also be posed as an object detection problem. In this case, the objective is to obtain 2-D spatial bounding boxes enclosing sketch strokes corresponding to the atypical content. We adopt this approach because it is faster and more amenable to near real-time operation compared to segmentation. Handwritten text is the most common atypical sketch content class in Pictionary. Hence, it is reasonable to consider approaches solely designed for text detection in domains such as outdoor scenes and documents~\cite{jaderberg2014reading,Deng_2019,liu2019omnidirectional,Liu_2016,liao2018textboxes++,zhou2017east,Baek_2019_CVPR}. Similarly, detection-based approaches have been proposed for mixed graphic structures~\cite{8395187,schafer2019arrow,Elyan2020DeepLF}. However, graphic elements in these scenarios are  more structured compared to our Pictionary setting.

\noindent \textbf{Pictionary-like guessing games:} Borrowing terminology from the seminal work of von Ahn and Dabbish~\cite{10.1145/1378704.1378719}, Pictionary can be considered an `inversion game' with full transparency. Riberio and Igarashi~\cite{10.1145/2380116.2380154} employ a sketching-based interactive guessing game to progressively learn visual models of objects. A review of Pictionary-like word guessing games involving drawing can be found in the work by Sarvadevabhatla et al.~\cite{8509167}.  In general, most of the existing works are confined to idealized toy settings~\cite{huang2020scones}, with some not even containing any sketching aspect~\cite{clark2021iconary,El-Nouby_2019_ICCV}. Unlike what we propose in this paper, they do not discuss the possibility of atypical content.

\section{Data Collection}
\label{sec:data-collection}

\subsection{Game Sessions}
\label{sec:game-browser}

Our browser-based game portal (Figures \ref{fig:pictguess},\ref{fig:pictguess2}) is compatible with mouse and touch inputs, scalable and can handle up to $50$ multiple concurrent Pictionary sessions. Consent is obtained and game instructions are provided when a player accesses the system for the first time. Players are assigned random names and paired randomly as Drawers and Guessers. The targets provided to the Drawers are sampled from a dictionary of $200$ guess phrases. We re-emphasize that the target phrases can be nouns (e.g. airplane, bee, chair), verbs (e.g. catch, call, hang) or adjectives (e.g. happy, lazy, scary). To ensure uniform coverage across the dictionary, the probability of a guess phrase being selected for a session is inversely proportional to the number of times it has been selected for elapsed sessions. The game has a time limit of $120$ seconds. The game ends when the Guesser enters a word deemed `correct' by the Drawer or when the time limit is reached.

For the Guesser, a text box is provided for entering guess phrases. For the Drawer, the interface provides a canvas with tools to draw, erase and highlight locations (via a time-decaying spatial animation `ping') for emphasis (see Fig.~\ref{fig:pictguess}). In addition, \faThumbsOUp $\:$ and \faThumbsODown $\:$ buttons enable Drawer to provide `hot/cold' feedback on guesses. A question (\faQuestionCircleO) button is provided to the Guesser for conveying that the canvas contents are not informative and confusing. The canvas strokes are timestamped and stored in Scalable Vector Graphic (SVG) format for efficiency. In addition to canvas strokes (drawing and erasure related), secondary feedback activities mentioned previously (\faThumbsODown$\:$,\faThumbsOUp $\:$, \faQuestionCircleO, highlight) are also recorded with timestamps as part of the game session.

Via our portal, we successfully gathered $3220$ timestamped episodes of diverse, realistic game play involving a total of $479$ participants in a large age range ($14$ years to $60$ years) and educational demographics (middle and high school students, graduate and undergraduate university students and working professionals). Please refer to project page for sample videos of game sessions, architectural overview of the game play system and plots with additional game session statistics.

\begin{table}[!t]
    \resizebox*{\linewidth}{!}
    {
    \begin{tabular}{llccc}
    \toprule
    \multicolumn{2}{c}{Sketch Content Type}& Number of & Number of & Number of \\
    \multicolumn{2}{c}{class} & occurrences & sessions & target phrases \\
    &&&containing&containing\\
    \hline
    \multicolumn{2}{l}{\textbf{\textit{Text}}}&\textbf{2419} & \textbf{478} & \textbf{180} \\
    \cline{3-5}
    &Individual letter&2244& 460& 178\\
    &Running hand&175& 103& 81\\
    \hline
    \multicolumn{2}{l}{\textbf{\textit{Numbers}}}&\textbf{331}& \textbf{73}& \textbf{28}\\
    \hline
    \multicolumn{2}{l}{\textbf{\textit{Circles}} }& \textbf{110}& \textbf{90}&\textbf{ 67}\\
    \hline
    \multicolumn{2}{l}{\textbf{\textit{Iconic}}} & \textbf{750} &\textbf{377} &\textbf{ 147}\\
    \cline{3-5}
    &Arrow & 497& 292& 129\\
    &Question mark & 158& 116& 78\\
    &Miscellaneous&95&54& 37\\
  \bottomrule
\end{tabular}
}
\caption{Statistics of atypical sketch content categories in game sessions.}
\label{tab:ano}
\end{table}

\subsection{Atypical Data}
\label{sec:atypical-data}
An atypical sketch content instance can be thought of as a subsequence of sketch curves relative to the larger sequence of curves that comprise the game session. We first describe the categories of atypical content usually encountered in Pictionary sessions:
 \begin{itemize}
     \item \textit{Text}: Drawer directly writes the target word or hints related to the target word on the canvas.
     \item \textit{Numerical}: Drawer writes numbers on canvas.
     \item \textit{Circles}: Drawers often circle a portion of the canvas to emphasize relevant or important content.
     \item \textit{Iconic}: Other items used for emphasizing content and abstract compositional structures include drawing a question mark, arrow and other miscellaneous structures (e.g. double-headed arrow, tick marks, addition symbol, cross) and striking out the sketch (which usually implies negation of the sketched item).
 \end{itemize}

 Examples can be viewed in Fig.~\ref{fig:ano-ex}. It is important to remember that we consider only \textit{Text} writing as a rule violation in Pictionary. Other categories mentioned above are atypical but their presence is not considered a violation of game rules.

 \begin{figure}[!t]
   \includegraphics[width=\linewidth]{labelling-tool.pdf}
   \caption{An illustration of annotation using the \textsc{CanvasDash} interface.
   }
   \label{fig:Dashboard}
 \end{figure}
To annotate atypical content, we use our custom-designed, browser-based annotation and visualization tool dubbed \textsc{CanvasDash} (see Fig.~\ref{fig:Dashboard})  - please refer to project page for details.

Using the described annotation procedure, we obtain our atypical Pictionary sketch dataset \textsc{AtyPict}. The occurrence statistics of atypical sketch categories across game sessions can be viewed in Table~\ref{tab:ano}. Representative visual examples can be viewed in Fig.~\ref{fig:ano-ex}.  Although we had earlier defined atypical content in terms of curve subsequences, the illustrations in Fig.~\ref{fig:ano-ex} show that the content can have a defined 2-D spatial extent and context relative to the canvas. For ease of processing, we consider this latter interpretation. In other words, we consider atypical content instances to be category-labelled 2-D spatial patterns sketched over time on the canvas. In terms of Fig.~\ref{fig:ano-ex}, detecting such content therefore corresponds to accurately detecting and categorizing 2-D spatial extents of entities shown color-coded on the canvas.

\begin{figure*}[t!]
  \includegraphics[width=0.8\textwidth]{CanvasNet.pdf}
  \caption{Architecture of \textsc{CanvasNet} deep neural network. Refer to Sec.~\ref{sec:canvasnetdeepnet} for details.}
  \label{fig:networkdiagram}
\end{figure*}

\section{CanvasNet}
\label{sec:canvasnet}

An effective approach for detecting atypical sketch instances needs to tackle the diversity in scale and appearance of various categories - a glance at Fig.~\ref{fig:ano-ex} makes this amply clear. In addition, the approach needs to utilize spatial context and be robust to the presence of similar looking yet semantically distinct regular (sketch) canvas content. To meet these requirements, we cast the problem as image-based object detection. In our case, the drawing canvas containing accumulated sketch strokes is the image. Any atypical content instances (e.g. \textit{Text}) present are considered spatially localized 2-D objects to be detected. For the object detection, we design a novel deep neural network which we dub \textsc{CanvasNet}. Before delving into details of \textsc{CanvasNet}, we first describe the data setup employed for its training and evaluation.

\subsection{Data Preparation}
\label{sec:strokedatapreparation}

As mentioned previously, the canvas elements for a given game session are represented as timestamped SVG curve elements. Each SVG element is either a  drawing stroke or an erasure stroke. We group drawing strokes into subsequences which are separated by erase strokes. The curves are converted to a 2-D point sequence representation and adaptively downsampled into line-based strokes using Ramer–Douglas–Peucker~\cite{RAMER1972244} algorithm  ($\epsilon=2$). The strokes are rendered on a $512 \times 512$ 2-D canvas with a stroke thickness of $4$ for the purpose of data annotation and representation. Note that a drawing stroke either belongs to one of the atypical classes (Sec.~\ref{sec:atypical-data}) or is a normal sketch stroke.
The spatial extents of labelled stroke subsequences are used to automatically generate ground-truth data for training and evaluation of \textsc{CanvasNet} (Sec.~\ref{sec:canvasnet}). Examples of ground-truth bounding boxes for atypical categories can be seen as solid (non-dashed/non-dotted) rectangles in Fig.~\ref{fig:CanvasNetDetections}.

\noindent \textbf{Data Augmentation:} The number of game sessions containing atypical instances are considerably smaller compared to the total number of game sessions. To increase the amount of data available for deep network training in a realistic manner, we first isolate atypical instance stroke subsequences. We sample from this set and add the resulting subsequences to other sessions which share the same target phrase, but do not contain any atypical content. The atypical content subsequences are also randomly rotated and localized carefully. This ensures they are spatially disjoint from strokes of the reference game sessions (which originally lack such atypical entities) -- examples can be viewed in project page.

\subsection{\textsc{CanvasNet} Deep Network}
\label{sec:canvasnetdeepnet}

Inspired by the success of deep networks which attempt to detect text in photos~\cite{liu2016ssd,shen2017dsod,liao2018textboxes++}, we adopt a similar efficient approach for our \textsc{CanvasNet} deep network to detect atypical sketch instances on a drawing canvas. \textsc{CanvasNet} consists of two stages.

\textbf{Feature Extractor:} The first stage consists of a stem block containing three $3 \times 3$ unit separable convolution layers~\cite{chollet2017xception}. Our choice of seperable convolution layers is motivated by the reduction in number of parameters and operations involved. The first layer in stem block uses a stride of $2$ to downsample the input features. The stem block is then followed by a repeating three segment structure consisting of (i) a $2 \times 2$ max pooling layer (ii) a $6$ layered dense block~\cite{huang2017densely} with a growth rate of $48$. Each of the $6$ layers consists of a $1 \times 1$ separable convolution followed by a $3 \times 3$ separable convolution (iii) a $1 \times 1$ separable convolution layer denoted as SC Transition layer. This three segment structure is then repeated three times with similar parameters - see project page for additional details.

\textbf{Multi-scale Predictor:} Atypical content (e.g. handwritten text, arrows) can occupy varying amounts of drawing canvas area, to detect them we use multi-scale predictor for prediction on multiple scales of feature maps~\cite{liao2018textboxes++, shen2017dsod}. For our setting, we use a customized multi-scale predictor for both handwritten text and non-text classes. The second stage sub-network is responsible for generating multi-scale feature maps and generating predictions over each of the feature maps. The output of third segment structure of the Feature Extractor is considered the first scale of the multi-scale feature map. The other feature map scales are obtained as outputs of successive downsampling blocks applied to the Feature Extractor's output (Fig.~\ref{fig:networkdiagram}) - see project page for additional details. The feature maps are passed individually through a prediction block. The multi-scale prediction features are concatenated and non-maximal suppression is applied to generate the final bounding box predictions.

\textit{Prediction block:} This consists of a $3 \times 5$ separable convolution layer, followed by a fully connected layer comprising the prediction. The rectangular filter dimensions ($3 \times 5$) used in the block ensure that elongated objects can be detected reliably.

\textit{Anchor boxes with vertical offsets:} Among the atypical object categories, words have larger aspect ratios and range of box orientation. Therefore, we set anchor aspect ratios to $1, 2, 3, 4, 5, 1/2, 1/3, 1/5$ with $\pm0.25$ as the vertical offset.

\noindent \textbf{Optimization:} We formulate the loss function for \textsc{CanvasNet} as a combination of a classification loss $L_{cls}$ and a bounding-box localization loss $L_{loc}$:

\begin{equation}
    L = \alpha \left( \frac{1}{N} \sum_{i=1}^{N} L_{cls}(P_i, G_i) \right) + \left( \frac{1}{M} \sum_{i=1}^{N} \sum_{j=1}^{c} G_{ij}*L_{loc}(B_{i}, B^{gt}_{i}) \right)
\end{equation}

where $G$ is the ground truth label matrix ($G_{ij} = 1$ if $i$-th anchor box belongs to atypical category $j$, else $G_{ij} = 0$), $P$ is the predicted confidence score matrix ($P_{ij}$ indicates the confidence score that $i$-th anchor box belongs to category $j$), this means $P_i$ is a row of confidence scores for $i$-th box to belong to various atypical categories, similarly $G_i$ is a row where values for all atypical category is $0$ except for one(to which $i$-th box belongs).
In the above loss function formulation, $B^{gt}_i$ and $B_i$ denotes the ground truth and predicted offsets for the $i$-th anchor box. $N$ is the total number of anchor boxes, $M$ is the total number of ground truth anchor boxes not belonging to background class, $c$ represents the number of atypical categories. We use focal loss\cite{lin2017focal} for $L_{cls}$ which prioritizes a sparse set of hard examples and prevents large number of negatives from overwhelming the detector. For localisation task, we adopt Distance-IoU Loss~\cite{zheng2020distance}:

\begin{equation*}
  \mathcal{L}_{D I o U}=1-I o U+\frac{\rho^{2}\left(\mathbf{b}, \mathbf{b}^{g t}\right)}{c^{2}}
\end{equation*}

where $b$ and $b^{gt}$ denote the central points of predicted and ground truth bounding box, $\rho^{2}\left(\mathbf{b}, \mathbf{b}^{g t}\right)$ gives the square of the Euclidean distance between them, and $c$ is the length of the diagonal of the smallest enclosing box covering the two bounding boxes.

\section{DrawMon}
\label{sec:drawmon}

\begin{figure}[t!]
  \centering
    \includegraphics[width=0.4\textwidth]{PictGuess_architecture.png}
    \caption{System architecture for Pictionary game setup and \textsc{DrawMon} (Sec.~\ref{sec:drawmon}).}
    \label{fig:pictguess2}
\end{figure}

\begin{figure*}[!t]
  \includegraphics[width=0.8\textwidth]{DrawMon.pdf}
  \caption{A pictorial overview of \textsc{DrawMon} - our distributed atypical sketch content response system (Sec.~\ref{sec:drawmon}). Also see Fig.~\ref{fig:pictguess2} for additional architectural details.}
  \label{fig:judge}
\end{figure*}

\begin{table*}[!t]
    \caption{Results for atypical content detection}
    \begin{subtable}[t]{.36\textwidth}
        \caption{\textsc{CanvasNet} performance for atypical content classes. IoU=0.5 refers to detection threshold. Refer to Sec.~\ref{sec:evalprotocol} for details.}
     \label{tab:multiclass}
        \raggedright
        \resizebox{\linewidth}{!}
{
           \begin{tabular}{c|cc}
\toprule
 \multirow{2}{*}{Atypical Content Category}& \multicolumn{2}{c|}{IoU=0.5} \\
 \cline{2-3}
      &  mAP & mAR \\% & Precision & Recall \\
\hline
      Text & 0.58 & 0.80 \\
      Number & 0.44 & 0.61 \\
      Icon & 0.55 & 0.68 \\
      Circle & 0.72 & 0.85 \\
      \bottomrule
    \end{tabular}
    }
    \end{subtable}%
    \hfill
   \begin{subtable}[t]{.6\textwidth}
        \raggedleft
        \caption{Performance comparison with baselines. mAP = mean Average Precision, mAR = mean Average Recall, \#Parameters = the number of trainable weights in the corresponding deep network in millions, ADT = average detection time per image in milliseconds.}
  \label{tab:baselines}
        \resizebox{\linewidth}{!}
{
        \begin{tabular}{l|cc|cc|c|c}
    \toprule
    \multicolumn{1}{c|}{\multirow{2}{*}{Method}} & \multicolumn{2}{c|}{Text only} & \multicolumn{2}{c|}{Multiclass} & \# Parameters & ADT\\
    \cline{2-5}
    & mAP & mAR & mAP & mAR & M=million & (m.sec)\\
    \hline
    \textbf{\textsc{CanvasNet}} & \textbf{0.78} & \textbf{0.90 } & \textbf{0.41} & \textbf{0.53} & 1.90 M & 35\\
    \hline
    BiLSTM+CRF~\cite{darvishzadeh2019cnn} & 0.06 & 0.04 & 0.02 & 0.03 & 0.01 M &85\\

    SketchsegNet+\cite{qi2019sketchsegnet+} & 0.56 & 0.32 & 0.04 & 0.11 & 3.90 M &21\\

    Tiny-YOLOv4\cite{wang2020scaled} & 0.26 & 0.65 & 0.31 & 0.51 & 5.88 M & 40\\

    TextBoxes++\cite{liao2018textboxes++} & 0.41 & 0.65 & 0.25 & 0.39 & 29.31 M & 51\\

    DSOD\cite{shen2017dsod} & 0.43 & 0.66 & 0.25 & 0.40 & 17.49 M & 52\\

    CRAFT \cite{Baek_2019_CVPR} & 0.47 & 0.69 & 0.17 & 0.30 & 1.18 M & 34 \\

    \bottomrule
 \end{tabular}
 }
    \end{subtable}
\end{table*}

Consider a scenario with multiple online Pictionary game sessions in progress. We require a  framework for automatic and concurrent monitoring of these game sessions for any atypical activities (e.g. a rule violation such as writing text on canvas). Such a framework needs to be reliable, scalable and time-efficient. To meet these requirements, we propose \textsc{DrawMon} - a distributed alert generation system (see Fig.~\ref{fig:judge}). Each game session is managed by a central Session Manager which assigns a unique session id (Fig.~\ref{fig:pictguess2}). For a given session, whenever a sketch stroke is drawn, the accumulated canvas content (i.e. strokes rendered so far) is tagged with session id and relayed to a shared Session Canvas Queue. For efficiency, the canvas content is represented as a lightweight Scalable Vector Graphic (SVG) object. The contents of the Session Canvas Queue are dequeued and rendered into corresponding $512 \times 512$ binary images by Distributed Rendering Module  in a distributed and parallel fashion. The rendered binary images tagged with session id are placed in the Rendered Image Queue. The contents of Rendered Image Queue are dequeued and processed by Distributed Detection Module. Each Detection module consists of our custom-designed deep neural network \textsc{CanvasNet} which processes the rendered image as input. \textsc{CanvasNet} outputs a list of atypical activities (if any) along with associated meta-information (atypical content category, 2-D spatial location).

The outputs from multiple distributed \textsc{CanvasNet} instances within the Distributed Detection Module are routed to the Alert Generator Module. An activity Record Table within this module records information related to ongoing game sessions and atypical content instances. This table is analyzed with respect to a Rule Base sub-module which generates appropriate alerts and relays them to the
appropriate game sessions. Since rule violations are of predominant interest, other atypical content  alerts can be filtered out. Incoming alerts are finally displayed on the game session user interface (UI)  - see Fig.~\ref{fig:pictguess}.

Also note that two manually-controlled mechanisms related to alert generation exist within the game UI. The Guesser player can press a button labelled `Drawer is violating rule!'. This  simply generates the alert (but does not highlight the canvas location where violation occurs). On the Drawer's side, the button `False Alarm' can be used to dismiss false positive alerts (see Fig.~\ref{fig:pictguess}). In the current deployment, we utilize the \textit{Text} detection variant of \textsc{CanvasNet} to detect text writing event on canvas.

\section{Experiments and Results}
\label{sec:experiments}

We first describe the experiments and results for atypical content detection. Following standard machine learning protocols, we divide data into training, validation and test splits. For each target phrase, the sessions containing atypical content are randomly split in the ratio $70$ (train) : $15$ (validation) : $15$ (test). Since we perform data augmentation on atypical content-free sessions, we divide such sessions in the aforementioned ratio as well. The respective data splits are combined to obtain the final groups.

\subsection{Atypical content detection}
\label{sec:evalprotocol}

\noindent \textbf{Baselines:} All along, our approach for detecting atypical activities treats the canvas as a 2-D image. In effect, the game session is considered to be a video-like frame sequence of 2-D canvas images. We also consider alternate approaches wherein the game session is processed as a sequence of curves. Each curve is labelled either as a regular sketch stroke or one associated with an atypical content category. Briefly, the baselines we consider are the following:
\textit{BiLSTM+CRF~\cite{darvishzadeh2019cnn} - } This classifies each stroke in a sketch sequence as one of the atypical classes using a bidirectional Long Short Term Memory (BiLSTM) neural network~\cite{hochreiter1997lstm} and Conditional Random Field (CRF)~\cite{sutton2006introduction}. The input to this model is a set of hand-crafted features of the strokes. \textit{SketchSegNet+~\cite{qi2019sketchsegnet+} - } This classifies each point in a sketch sequence as one of the atypical classes using bidirectional LSTM. For image based models, we train appropriately modified versions of two state-of-the-art generic object detectors -- \textit{DSOD~\cite{shen2017dsod}} and \textit{Tiny-YOLOv4}~\cite{wang2020scaled}. We also train modified versions of two popular scene text detection models -- \textit{TextBoxes++~\cite{liao2018textboxes++}} and \textit{CRAFT~\cite{Baek_2019_CVPR}}. Please see project page for architectural details of baselines.

\noindent \textbf{Evaluation Protocol:} We conduct evaluation using two protocols. In the first protocol, we consider models trained to detect all atypical content classes. To score performance, we use the standard object detection measures --  mean-average-precision (mAP) and  mean-average-recall (mAR)~\cite{zhou2019iou}. These measures are typically reported on a $[0,1]$ scale -- larger the better. mAP and mAR are reported at an Intersection-over-Union (IoU) threshold of $0.5$. In other words, an overlap of 50\% or greater between predicted bounding box and ground-truth bounding box is deemed correct (assuming predicted category label also matches).

\noindent \textbf{Training:} For training \textsc{CanvasNet}, we employ Adam optimizer~\cite{kingma2014adam} with a mini-batch size of $8$, the exponential decay rate for $1st$ and $2nd$ moment estimates set to $0.9$ and $0.999$ respectively. We stop training after $50$ epochs. Classification loss weight ($\alpha$) is set to $1000$ for quick convergence. The training takes approximately $4.5$ hours on two GTX 1080Ti 11GB GPUs. We use \textit{mish}~\cite{misra2020mish} activation function since it provides fast convergence and improved performance compared to the standard \textit{relu} activation. For training the text-only model, we used \textit{hard mining training regime} wherein false-positive samples having IoU overlap value with ground truth in the range $[0.45, 0.55]$ are chosen for mining. We observed significant increase in $mAP$, $mAR$ compared to regular training regime due to this regime. For training the multiclass model, we used mini-batch re-sampling with class balancing to counteract the presence of imbalanced per-class sample distribution.

\noindent \textbf{Detection Results:} \textsc{CanvasNet}'s performance for all the atypical categories can be viewed in Table~\ref{tab:multiclass}. Fig.~\ref{fig:CanvasNetDetections} depicts examples of \textsc{CanvasNet} detections, including some failure cases. We conducted ablation experiments to determine the relative importance of our architectural and optimization choices. Details of the ablation configurations and results can be viewed in the project page. The comparative evaluation of \textsc{CanvasNet} with baseline approaches is summarized in  Table~\ref{tab:baselines}. Note that the comparison also includes compute-related aspects (number of trainable parameters in the approaches and average detection time).

Table~\ref{tab:multiclass} shows \textsc{CanvasNet}'s performance for various atypical content categories. The consistent depictions of \textit{Circle} enables good performance for the category. Detecting isolated numbers is slightly more challenging. Empirically, we observed that sketched content resembling letters (e.g. a mountain sketch) or numerals (e.g. vertical bar groups) accounted for most of the false positive detections. Text spanning a significantly large extent of the canvas and unusually oriented numbers accounted for majority of the missed detections (false negatives). Performance scores for ablative variants of \textsc{CanvasNet} are included in project page.

From Table~\ref{tab:baselines}, we see that \textsc{CanvasNet} clearly outperforms a variety of   baseline approaches (Sec.~\ref{sec:evalprotocol}). This is predominantly due to the carefully considered architectural and optimization choices in designing \textsc{CanvasNet}. The results also illustrate the superiority of image-based approaches compared to the sketch stroke processing approaches (\textit{BiLSTM
+CRF, SketchSegNet+}). Keeping the rule-violation detection scenario in mind, we also trained variants designed to detect the single class \textit{Text}. As the `Text only' column in Table~\ref{tab:baselines} shows, \textsc{CanvasNet} remains the best performer. Consequently, we utilize this model variant as part of \textsc{DrawMon} in our game deployment scenario (Sec.~\ref{sec:drawmon}). From the table (column named `Parameters'), we also note that \textsc{CanvasNet} achieves its superior performance despite containing a smaller number of parameters relative to most of the baselines.

\begin{figure*}[!t]
    \centering
    \includegraphics[width=0.9\textwidth]{Multiclass_Samples.jpeg}
    \caption{Examples of atypical content detection by \textsc{CanvasNet}. False negatives are shown as dashed rectangles and false positives as dotted rectangles. Color codes are: \textcolor{red}{text}, \textcolor{cyan}{numbers}, \textcolor{ForestGreen}{question marks}, \textcolor{blue}{arrows}, \textcolor{Maroon}{circles} and other \textcolor{orange}{icons} (e.g. tick marks, addition symbol). }
    \label{fig:CanvasNetDetections}
\end{figure*}

\subsection{DrawMon User Study Experiments}
\label{sec:drawmonuserstudy}

To quantify the efficacy of \textsc{DrawMon}, we analyzed game session data with \textsc{DrawMon} deployed to detect text. The canvas contents are relayed to \textsc{DrawMon} every $1$ second. We deployed $4$ \textsc{CanvasNet} instances within the Distributed Detection Module on two 2080Ti GPUs alongside $16$ worker processes for svg to image conversion. The combined peak usage of GPU RAM was $20$ GB while peak CPU RAM usage was $15$ GB. $23$ participants ($11$ male, $12$ female) in the age group $19-25$ (mean=$20.8$, std.=$3.1$), recruited using social media and from the institution's student pool, participated in the study. Each session had an average duration of $47.5$ sec (std.= $37.4$) with the maximum being $120$ seconds. Over the study period, the maximum number of concurrent game sessions managed by \textsc{DrawMon} was $4$. From the resulting set of $145$ game sessions, $69$ contained atypical text activities. During the sessions, we recorded timestamped alerts from \textsc{DrawMon}, false alarm notifications by the Drawer player and rule violation notifications from the Guesser player. The results from the study are summarized in Table~\ref{tab:userstudy}. To determine \textsc{DrawMon}'s throughput, we measured two quantities. The first, processing time (\texttt{p-time}), is the average elapsed time between the canvas representation being sent to \textsc{DrawMon} and receiving an alert. In case no alerts were generated, the timestamp corresponding to end of \textsc{CanvasNet} processing was considered. From our data, \texttt{p-time} was $0.4$s. The other measurement was the maximum number of concurrently active sessions (\texttt{n-sess}) -- this was $4$. Defining the effective throughput rate to be \texttt{tpr} = \texttt{p-time}/\texttt{n-sess}, we obtain an average processing rate of $10$ items per second.

\begin{table}[!t]
\captionof{table}{User study statistics with \textsc{DrawMon} deployed.}
\label{tab:userstudy}
\resizebox{\linewidth}{!}
{
\centering
    \begin{tabular}{lc}
    \toprule
    Game Event Type & Count \\
    \midrule
    (True Positive) \textsc{DrawMon} generates `Rule Violation' alert. Drawer doesn't press `False Alarm' button. & $62$ \\
    (False Positive) \textsc{DrawMon} generates `Rule Violation' alert. Drawer presses `False Alarm' button.   & $32$ \\
    (False Negative) No `Rule Violation' alert. Guesser presses `Drawer is violating rule' button. & $6$ \\
    \bottomrule
    \end{tabular}
}
\end{table}

 \noindent \textbf{Results and Analysis:} The results from \textsc{Drawmon} deployment user study are summarized in Table~\ref{tab:userstudy}. From the table, we see that a significant fraction of \textsc{DrawMon} generated alerts are valid (see `True Positives'). From the results, \textsc{DrawMon}'s precision is $0.66$ while recall is $0.91$. Post the user study, we conducted a brief survey with Likert-type questions on a $1$ to $5$ scale with $5$ being the best. `Q: How responsive was \textsc{DrawMon} to valid rule violations?' : The average score was $3.63$ (s.d.=$0.74$), indicating reasonably high system throughput despite multiple concurrent sessions. This is also supported by the recorded throughput rate (\texttt{tpr}) mentioned previously in this section. `Q: How was the overall game experience?' : The score was $3.91$ (s.d.=$0.60$), suggesting a positive session experience and  satisfaction with rule violation detection and response mechanisms.

User study plots and sample videos of game sessions with \textsc{DrawMon} in action can be viewed in project page.

\subsection{Application Scenarios}
\label{sec:limitationsotherscenarios}

\noindent \textbf{Application Scenarios:} Although we have used Pictionary as a use case scenario, we expect \textit{DrawMon} to be suitable for other shared and interactive whiteboard scenarios. For instance, in a writing related setting, the notion of atypical categories can be the exact opposite of Pictionary scenario: text on canvas would be routine while drawings might be considered abnormal. This can be tackled by appropriate data labelling, for e.g. using our \textsc{CanvasDash} annotation tool, and subsequently retraining \textsc{CanvasNet} deep network. In another scenario, consider participants grouped into teams for a collaborative scene drawing task~\cite{davidson2021drawn,10.1145/3313831.3376517}. \textsc{DrawMon}, using a \textsc{CanvasNet} configured for sketched scene recognition~\cite{zou2018sketchyscene}, can alert the instructor on progress and task completion. For this task, a \textsc{CanvasNet} instance trained to recognize individual objects and iconic components from our dataset (e.g. arrows, `addition mark') could also be included as additional detection component for expanding the detection capability.

\section{Conclusion and Future Work}

\textsc{DrawMon} is a distributed framework for monitoring multiple shared interactive whiteboards for detecting atypical content.  We use a Pictionary-like sketching game as the use case scenario.  \textsc{DrawMon} is enabled by a number of equally important lateral contributions - (i) \textsc{CanvasDash} - an intuitive dashboard UI for annotation and visualization (ii) \textsc{AtyPict} - a first of its kind dataset for atypical sketch content (iii) \textsc{CanvasNet} - a deep neural network for atypical content detection. Together, these reusable contributions create the possibility of developing similar frameworks for other shared and interactive whiteboard scenarios. Apart from atypical content detection, we expect our game session dataset to be a valuable resource in itself for analyzing player characteristics and strategies in communication restricted non-adversarial games~\cite{10.1145/3313831.3376316}. In addition, we plan to develop practical AI agents which can mimic human Pictionary players in a more interactive, realistic manner compared to existing non-interactive works~\cite{bhunia2020sketch,huang2020scones}.

\bibliographystyle{ACM-Reference-Format}


\section{Introduction}
The main challenge in perturbative Quantum Field Theory at colliders is the computation of multiloop scattering amplitudes. In order to face this difficult task, a novel methodology has been developed, the loop-tree duality (LTD)~\cite{Catani:2008xa,Rodrigo:2008fp,Bierenbaum:2010cy,Bierenbaum:2012th,Buchta:2014dfa,Tomboulis:2017rvd,Runkel:2019yrs,Capatti:2019ypt} which opens any Feynman loop diagram into a sum of connected trees. The effort to deepen in the LTD is reflected in their appealing features~\cite{Buchta:2015wna,Hernandez-Pinto:2015ysa,Sborlini:2016gbr,Sborlini:2016hat,Jurado:2017xut,Driencourt-Mangin:2017gop,Driencourt-Mangin:2019yhu,Plenter:2020lop,Aguilera-Verdugo:2019kbz,Prisco:2020kyb}, furthermore, the latest progress in this framework includes a clever reformulation which was first presented in Ref.~\cite{Verdugo:2020kzh}. This result allowed to exploit its most remarkable property, the existence of a manifestly causal representation, opening the path for a significant evolution~\cite{snowmass2020,Aguilera-Verdugo:2020kzc,Aguilera-Verdugo:2020nrp,Ramirez-Uribe:2020hes,Sborlini:2021owe,Bobadilla:2021pvr,TorresBobadilla:2021ivx,Aguilera-Verdugo:2021nrn,Ramirez-Uribe:2021ubp,RamirezUribe:2021nlj,Uribe:2021cdb,Ramirez-Uribe:2022gxz}.

The possibility to work in a causal LTD scenario, in which noncausal singularities are absent, gives the advantage to work with more numerically stable integrands~\cite{Aguilera-Verdugo:2020kzc,Ramirez-Uribe:2020hes}.
To achieve a causal LTD representation it is required to identify the causal configurations of the multiloop topology of interest.
In this presentation, we review the solution presented in Ref.\cite{Ramirez-Uribe:2021ubp,RamirezUribe:2021nlj,Uribe:2021cdb} for the unfolding of those configurations fulfilling causal conditions by the application of a modified Grover's quantum algorithm~\cite{Grover:1997fa} for querying multiple solutions over unstructured databases~\cite{Boyer:1996zf}. A variational quantum eigensolver approach has recently been presented in Ref.~\cite{Clemente:2022nll}.

\section{LTD framework}
The LTD representation of any multiloop scattering amplitude is computed by the iterative evaluation of the Cauchy's residue theorem in the terms of nested residues~\cite{Verdugo:2020kzh,Ramirez-Uribe:2020hes}.
Regarding the causal LTD expression, it is obtained by summing over all the nested residues and followed by an ingenious rearrangement~\citep{Aguilera-Verdugo:2020kzc,Ramirez-Uribe:2020hes} that allows to achieve the following causal dual form,
\begin{equation}
{\cal A}_D^{(L)} = \int_{\vec \ell_1 \ldots \vec \ell_L}
\frac{1}{x_n} \sum_{\sigma  \in \Sigma} \frac{{\cal N}_{\sigma(i_1, \ldots, i_{n-L})}}{\lambda_{\sigma(i_1)}^{h_{\sigma(i_1)}} \cdots \lambda_{\sigma(i_{n-L})}^{h_{\sigma(i_{n-L})}}}
+ (\lambda_p^+ \leftrightarrow \lambda_p^-)~,
\label{eq:AD}
\end{equation}
where $x_n = \prod_n 2\qon{i}$.
This expression is written in terms of causal propagators, $1/\lambda_p^\pm$, with $\lambda_p^\pm = \sum_{i\in p} \qon{i} \pm k_{p,0}$ , on-shell energies $\qon{i}=\sqrt{\mathbf{q}_i^2+m_i^2-\imath 0}$ and $k_{p,0}$ a linear combination of the energy components of external momenta.
After every propagator in the set $p$ is set on shell, and according the sign of $k_{p,0}$, either $\lambda_p^-$ or $\lambda_p^+$ becomes singular.
The collection of entangled causal propagators is encoded $\Sigma$, representing all combinations of compatible causal thresholds.

To work in a general context, we recall the concepts of edges and eloops~\cite{TorresBobadilla:2021ivx,Sborlini:2021owe}. Edges are defined as the union of an arbitrary number of propagators connecting two interaction vertices.
The algorithm and the implementation is presented in terms of eloops, or loops made of edges, given that in a causality scenario the causal singular configurations only consider the setting of having all the momentum flow of all the propagators in an edge aligned in the same direction.

\section{ Feynman loop integrals in a quantum computer}

The feasibility of associating quantum computing and Feynman loop integrals is due to the fact that the two on-shell states of a Feynman propagator can be stored in a qubit. The initial momentum flow for a specific topology is represented with the state $\ket{1}$, whereas $\ket{0}$ encodes those states with inverse flow orientation.

\subsection{Modified Grover's quantum algorithm}
\begin{figure}[t]
\raisebox{15pt}{\includegraphics[scale=0.2]{2eloop-2edges.pdf}} \qquad
\includegraphics[scale=0.30]{2eloop_qcx.pdf}
\centering
\caption{Six-edge two-eloop topology (left) and its corresponding quantum circuit (right) used to bootstrap its causal configurations.
\label{fig:two-eloop}}
\end{figure}

The global structure of Grover's quantum algorithm relies in three main components: the uniform superposition of all possible states, an oracle and a diffuser operator. Regarding Feynman loop integrals, the states refer to the associated configuration of the diagram. The total number of possible configurations is $N=2^n$, with $n$ the number of edges. Additionally, the causal configurations represent the winning states and the noncausal configurations the orthogonal states.

The uniform superposition of all the $N$ states is written in terms of the uniform superposition of the causal $\left(~\ket{w}~\right)$ and noncausal $\left(~\ket{q_\perp}~\right)$ states as: $\ket{q} = \cos \theta \, \ket{q_\perp} + \sin\theta \, \ket{w}$. The mixing angle between the causal and noncausal states is denoted by $\theta = \arcsin \sqrt{r/N}$, with $r$ the number of causal states.
The oracle operator, $U_w$, flips the state $\ket{x}$ if $x\in w$, and does nothing otherwise; the diffuser operator, $U_q$, amplifies the probability of these elements by performing a reflection around the initial state $\ket{q}$.
The iterative application of both operators $t$ times gives $(U_q U_w )^t \ket{q} = \cos \theta_t \, \ket{q_\perp } +  \sin \theta_t \, \ket{w}$, with $\theta_t = (2t +1) \, \theta$.

To define a proper number of iterations, $\theta_t$ has to fulfill that the noncausal state probabilities are much smaller than causal state probabilities.
Based on this condition, if the initial mixing angle is less or similar to $\pi /6$, Grover's quantum algorithm is considered a feasible framework.

Regarding the identification of causal singular configurations of multiloop Feynman diagrams, it is known from classical~\cite{Aguilera-Verdugo:2020kzc,Ramirez-Uribe:2020hes} and quantum~\cite{Ramirez-Uribe:2021ubp} computations that in most cases they do not fit in a favorable scenario for a direct application of Grover's algorithm.
A suitable adjustment is given by fixing one qubit to reduce the number of causal solutions, taking advantage that given one causal solution the mirror configuration, is also a causal solution~\cite{Ramirez-Uribe:2021ubp,RamirezUribe:2021nlj,Uribe:2021cdb,Ramirez-Uribe:2022gxz}; if the case requires it, also the total number of states can be increased by introducing an ancillary qubit in the $\ket{q}$ register~\cite{Nielsen2000}.

The modified Grover's quantum algorithm requires three registers and one extra qubit used as a marker in the oracle. The $n$ edges are encoded in $n$ qubits $q_i$;
the binary clauses compare the orientation of two adjacent edges through $c_{ij} \equiv (q_i = q_j)$ and $\bar c_{ij} \equiv (q_i \ne q_j)$ with $i,j \in\{0, \ldots, n-1\}$, which are stored in $\ket{c}$.
The loop clauses encoded in the register $\ket{a}$ combine the qubits from $\ket{c}$ to test if any subloop configurations generate a cyclic circuit.

The structure of the algorithm takes as a first step the initialization of all the registers involved.
The qubits encoding the edges are set in a uniform superposition, $\ket{q} = H^{\otimes n} \ket{0}$; the registers $\ket{a}$, $\ket{c}$ are set to the state $\ket{0}$; and the Grover's marker is set to, $\ket{out_0} = \left(\ket{0} - \ket{1}\right)/\sqrt{2}$.

The implementation of $\bar c_{ij}$ needs two CNOT gates acting between $q_i$, $q_j$ and a qubit in $\ket{c}$. In the case of $c_{ij}$, an extra XNOT gate is used to operate on the respective qubit in $\ket{c}$. The loop clauses are set by combining binary clauses and implemented through multicontrolled Toffoli gates.

To test the causal conditions we define $f(a,q)$. If the causal conditions are satisfied then $f(a,q) = 1$, if not $f(a,q) = 0$.
The implementation of the oracle, to identify and to mark, is given by
$U_w \ket{q} \ket{c} \ket{a} \ket{out_0} = \ket{q} \ket{c} \ket{a} \ket{out_0 \otimes f(a,q)}$, where $\ket{out_0 \otimes f(a,q)}$ is $-\ket{out_0}$ if $q \in w$, and $\ket{out_0}$ if $q \not\in w$.
The process continues by applying the oracle operations in the inverse order.
As a final step before measuring, the amplification of the probabilities is applied by implementing $U_q$ to $\ket{q}$. The function of the diffuser operator is taken from IBM Qiskit (\texttt{https://qiskit.org/}).

\subsection{Two-eloop topology with six edges}

 \begin{figure}[t]
\includegraphics[width=\textwidth]{2eloopHist.pdf}
\centering
\caption{Probabilities of causal and noncausal configurations of a two-eloop topology with six edges.
\label{fig:output}}
\end{figure}
The implementation of the algorithm described in the previous section is illustrated with the two-eloop topology depicted in Fig.~\ref{fig:two-eloop} (left), consisting of six edges.
From a classical computation it is found that $\theta \thickapprox \pi/3$, therefore, a suitable modification is required.
Halving the number of solution leads to $\theta \thickapprox \pi/5$, which allows to implement Grover's quantum algorithm, furthermore, under this arrangement only one iteration is needed to achieve a desire probability amplification.

The total number of qubits required for the implementation of the algorithm is sixteen, six for encoding the edges,  six to store the binary clauses, three to test the causal conditions and the one standing as the Grover's oracle marker.
The six-edge two-eloop diagram (Fig.~\ref{fig:two-eloop}) requires to test three subloops:
$a_0 = \neg \left( c_{01} \wedge c_{14}\wedge c_{45}\right)~, \quad a_1 = \neg \left( c_{01} \wedge \bar c_{13} \wedge c_{23} \right)~$ and $~a_2 = \neg \left( c_{23} \wedge c_{34} \wedge c_{45} \right)$.

The Boolean condition marking the causal configurations is $f^{(2)}(q,a)=(a_0\wedge a_1 \wedge a_2)\wedge q_0$, where $q_0$ is fixed.
The quantum circuit associated to the algorithm is depicted in Fig.~\ref{fig:two-eloop} (right), successfully implemented in the IBM Qiskit quantum simulator {\tt{qasm}}. The output shown in Fig.~\ref{fig:output} successfully identifies the expected twenty three causal configurations, corresponding to the forty six causal states when including the mirror states obtained by inverting the momentum flows.

\section{Conclusions}
A modified Grover’s quantum algorithm, for the identification of causal singular configuration of multiloop Feynman integrals was described.
The algorithm was illustrated with the six-edge two-eloop topology, implemented through the quantum simulator provided by IBM Qiskit, and successfully determined all the causal singular configurations.

The development of this algorithm represents a relevant resource to the LTD framework, providing an efficient procedure to search the causal singular configurations required to bootstrap the causal LTD representation of mutiloop scattering amplitudes.

\section*{Acknowledgments}
I want to thank G. Rodrigo for the guidance through the development of this work, and A. Renter\'ia for the companionship and support.
Support for this work has been received in part by MCIN/AEI/10.13039/501100011033, Grant No. PID2020-114473GB-I00, Consejo Nacional de Ciencia y Tecnología and Universidad Aut\'onoma de Sinaloa.


\section{Introduction}
\label{sec:introduction}
In recent times, sporadic infectious disease outbreaks have created massive disruptions not only in public health functioning but also in socio-economic developments. Incomplete knowledge of disease transmission mechanisms, and unavailability of effective medical treatment cause difficulties in controlling disease spread. Advances in modern transportation technology further accelerate the spatial spread of infection and complicate the management of disease spread in affected regions~\citep{brockmann2009human, findlater2018human}. Mobility can occur on different spatial scales within a country. The complex interconnections among multiple scales limit the predictability of spatio-temporal patterns of disease spread~\citep{alessandretti2020scales, barbosa2018human,schlapfer2021universal}. Alternatively, inter-country migration often acts to introduce new pathogens in non-affected countries and within-country mobility can then quickly spread the disease throughout the country~\citep{wilder2008geographic}. Inter-country migration occurs primarily in two ways: international air-travel and daily commuting between countries sharing a border. The availability of digital technology makes it possible to keep track of the fluxes of global air-travel passengers. Moreover, during an outbreak, it helps to implement measures like airport screening, quarantining upon arrival, and contact tracing, which reduce the risk of disease import due to air travel. On the other hand, for daily cross-country mobility, due to a lack of convenient real-time database management~\citep{cavallaro2019cross}, it is often difficult to estimate the mobility flow, and, in turn, to quantify the component of infection risk because of human mobility.

Mathematical models are important tools for assessing the risk of disease import from a neighboring country and investigating the consequences of different control strategies. Reliable information on cross-border mobility flow not only hints at the effectiveness of intervention strategies like lockdowns or border-control measures but is also a key input parameter for multi-national epidemiological models in evaluating the performance of spatially-explicit intervention strategies. In very few cases, mobility flux information can be obtained from technology companies, public transit data, census data, and survey data. However, these data are very expensive to obtain, and often not publicly available due to privacy laws. Therefore, estimation of cross-country mobility flux remains a challenging task in the context of infectious disease modelling.

In spatial epidemiological analyses, mobility flux is usually accounted for via the gravity model, the radiation model, or related approaches~\citep{balcan2009multiscale,jandarov2014emulating, xia2004measles, senapati2019impact, mari2012modelling,yang2015transmission, mertel2022fine}. However, despite their wide use in the literature, these models have certain limitations as they rely on tuneable parameters which can vary with spatial location~\citep{simini2012universal}. Similarly, several theoretical frameworks have been proposed to reconstruct the underlying connection topology from the observed dynamics of the focal system~\citep{timme2007revealing, casadiego2017model,shandilya2011inferring}. These methods have been mainly applied to those systems that exhibit oscillatory dynamics, such as coupled chaotic oscillators~\citep{gao2022autonomous, ching2015reconstructing,shi2021inferring}, EEG time series data~\citep{alves2022eeg}, or neuronal systems~\citep{sysoeva2021reconstruction}. Recently, \cite{hempel2020estimating} proposed a method to estimate the coupling strength between two regions from invasion time, which is the time taken to reach the focal region after an infection is seeded in a neighboring region. In practice, however, the invasion time will not likely be a reliable enough basis for extracting mobility flow. This is because disease monitoring systems typically face the greatest difficulties in the early phase of a disease outbreak due to incompleteness or delays in reporting.  Consequently, all existing methods for estimating spatial connectivity are either data-hungry or depend heavily on precise information related to data reporting. A need therefore exists for a practical method that overcomes these limitations and can estimate cross-border mobility from imperfect monitoring data.

To this end, we propose an intuitive approach based on maximum-likelihood estimation to quantify the cross-border mobility flux for a pair of regions situated in two different countries. Our method assumes the underlying disease dynamics is governed by a simple Susceptible-Infected-Recovered (SIR) model and estimates the mobility flux from the observed difference in the timing of peak of infection between two regions. We first theoretically investigate, under realistic epidemiological assumptions, how the observed difference in peak-timing can effectively recover the underlying mobility between two regions and how the estimate of mobility varies depending on the disease transmission rate. Next, as a case study, we use our method to estimate the cross-country mobility flux from COVID-19 incidence data for the pair of regions which are connected directly through multiple modes of transportation and are located on either side of the Poland-Germany border.

\section{Methods}\label{sec:methods}
\subsection{Deterministic model}
We consider a standard deterministic two-patch SIR model incorporating short-term or commuting-type migration. Each spatial unit in our study, such as a district or even a state which shares border, can be considered a patch. In each patch $i$, the total population is divided into three classes based on health status: susceptible ($S_{i}$), infected ($I_{i}$), and recovered ($R_{i}$). Since the patches are connected through migration, the disease dynamics in each patch has two components: (i) disease dynamics within the patch, and (ii) disease dynamics between the patches. The mathematical equations describing the mechanism of disease dynamics in this two-patch scenario is given as follows:
\begin{eqnarray}
\begin{array}{lll}
\displaystyle \frac{d S_{i}}{dt} &=& -\displaystyle  \sum_{j=1}^{2}\beta_{j}m_{ij}S_{i}\Bigg(\frac{\displaystyle \sum_{k=1}^{2}m_{kj}I_{k}}{\displaystyle \sum_{k=1}^{2}m_{kj}N_{k} } \Bigg),\\

\displaystyle \frac{d I_{i}}{dt} &=& \displaystyle  \displaystyle  \sum_{j=1}^{2}\beta_{j}m_{ij}S_{i}\Bigg(\frac{\displaystyle \sum_{k=1}^{2}m_{kj}I_{k}}{\displaystyle \sum_{k=1}^{2}m_{kj}N_{k} } \Bigg) - \gamma_{i} I_{i},\\

\displaystyle \frac{d R_{i}}{dt} &=& \displaystyle  \gamma_{i} I_{i}, \hspace*{1cm}~i~=~1,~2.
\end{array}
\label{EQ:Model}
\end{eqnarray}
The rate at which the susceptible population in patch $i$ gets infected in patch $j$ is given by the product of three terms:  the risk of infection in patch $j$ ($\beta_{j}$), the number of susceptibles from patch $i$ who are currently in patch $j$ ($m_{ij}S_{i}$), and the proportion of the population that is infected in patch $j$ ($\frac{\sum_{k=1}^{2}m_{kj}I_{k}}{ \sum_{k=1}^{2}m_{kj}N_{k} }$). Here, $N_{j}$ denotes the total population in patch $j$ and $\sum_{k=1}^{2}m_{kj}N_{k}$ gives the total population currently present in patch $j$ is called effective population of patch $j$~\citep{bichara2015sis}. The model,~\eqref{EQ:Model}, has two epidemiological parameters: transmission rate ($\beta_{i}$), and infectious period ($\frac{1}{\gamma_{i}}$). The matrix $m_{ij}$ is called the residence-time matrix, where the element $m_{ij}$ can be interpreted as the fraction of people from patch $i$ that visited patch $j$~\citep{bichara2015sis}. Since $m_{ij}$ denotes the fraction of mobility flux, we have the constraint $\sum_{j=1}^{2}m_{ij}=1$, for $i=1,2$.

For a fixed set of model parameters and initial conditions, we define $t_{d}$ as the difference between the peak timing of these two regions, which can be expressed:
\begin{equation*}
    t_{d}= |t_{1}^{\max} - t_{2}^{\max}|,
\end{equation*}
    where, $t_{1}^{\max}$ and $t_{2}^{\max}$ denote the time points at which the trajectories of the infected populations $I_{1}(t)$ and $I_{2}(t)$ reach their respective maxima.

\subsection{Stochastic model}
\label{subsec:stochastic_model}
Here we are interested in estimating the cross-border mobility flux from the observed difference in peak timing, $t_d$. It should be noted that in a deterministic setup, for given model parameters and initial conditions, we will always end up with exactly one $t_{d}$. Therefore, we introduce randomness into the model as a form of demographic stochasticity. Here demographic stochasticity refers to the fluctuations in the disease propagation process arising from the random nature of disease transmission and recovery events at the individual level. We consider a standard event-driven approach called the \textit{tau-leaping} method~\citep{gillespie2001approximate}. This method is a modified version of the Gillespie algorithm where the inter-event duration is kept fixed instead of following an exponential distribution, which increases computational efficiency~\citep{gillespie2001approximate,gillespie2003improved}.

Let us denote the number of transmission and recovery events at time $t$ for each patch $i$ (i=1, 2) by $\mathcal{E}_{\rm Trans}^{i}(t)$ and $\mathcal{E}_{\rm Rec}^{i}(t)$ respectively. The transition probabilities for the events $\mathcal{E}_{\rm Trans}^{i}(t)$ and $\mathcal{E}_{\rm Rec}^{i}(t)$ that occur during the small but fixed time interval $\delta t$ can be written as

\begin{eqnarray}
\begin{array}{lll}
\mathbb{P}(\delta \mathcal{E}_{\rm Trans}^{i}=1|S_{i}, I_{i}, R_{i})=\displaystyle  \displaystyle  \sum_{j=1}^{2}\beta_{j}m_{ij}S_{i}\Bigg(\frac{\displaystyle \sum_{k=1}^{2}m_{kj}I_{k}}{\displaystyle \sum_{k=1}^{2}m_{kj}N_{k} } \Bigg) \delta t + o(\delta t),\\

\mathbb{P}(\delta \mathcal{E}_{\rm Rec}^{i}=1|I_{i})= \gamma_{i} I_{i} \delta t + o(\delta t),
\end{array}
\end{eqnarray}
where, $\delta \mathcal{E}_{\rm Trans}^{i}=\mathcal{E}_{\rm Trans}^{i}(t+\delta t) - \mathcal{E}_{\rm Trans}^{i}(t),$ and
 $\delta \mathcal{E}_{\rm Rec}^{i}=\mathcal{E}_{\rm Rec}^{i}(t+\delta t) - \mathcal{E}_{\rm Rec}^{i}(t)$ for $i=1,2$.

For small $\delta t$, the increments $\delta \mathcal{E}_{\rm Trans}^{i}$ and $\delta \mathcal{E}_{\rm Rec}^{i}$ are approximated as Poisson distributed random variables:
\begin{eqnarray}
\begin{array}{lll}
\delta \mathcal{E}_{\rm Trans}^{i}\approx Poisson(\displaystyle  \sum_{j=1}^{2}\beta_{j}m_{ij}S_{i}\Bigg(\frac{\displaystyle \sum_{k=1}^{2}m_{kj}I_{k}}{\displaystyle \sum_{k=1}^{2}m_{kj}N_{k} } \Bigg) \delta t ),\\

\delta \mathcal{E}_{\rm Rec}^{i}\approx Poisson(\gamma_{i}I_{i} \delta t).
\end{array}
\end{eqnarray}

The population sizes in different compartments are then updated as,
\begin{eqnarray}
\begin{array}{lll}
S_{i}(t+\delta t)= S_{i}(t) -\delta \mathcal{E}_{\rm Trans}^{i} ,\\
I_{i}(t+\delta t)= I_{i}(t) +\delta \mathcal{E}_{\rm Trans}^{i} - \delta \mathcal{E}_{\rm Rec}^{i},\\
R_{i}(t+\delta t)=R_{i}(t) +\delta \mathcal{E}_{\rm Rec}^{i}.
\end{array}
\end{eqnarray}

Given model parameters $\beta_{i}, \gamma_{i}, N_{i}, I_{i}(0)$, for a particular residence-time matrix $m_{ij}$, we generate multiple stochastic realizations. In each realization, we calculate the difference in peak timing ($t_{d}$) and consequently the distribution of it. We show an example in Fig.~\ref{fig:fig_schematic} (A)-(B), where we consider the parameters as $\beta_1=0.30$, $\beta_2=0.20$, $\gamma_1=\gamma_2 = 0.07$, $N_1=10000$, $N_2=20000$, $I_{1}(0)=10$, $I_{2}(0)=1$ and generate multiple stochastic realizations. In each of the realizations, we calculate $t_d$ (see Fig.~\ref{fig:fig_schematic} (A)) and plot its distribution as a histogram ( see Fig.~\ref{fig:fig_schematic} (B)).

To produce the likelihoods, we first consider a range of values for the mobility parameter and divide it into 100 linearly-spaced values. Then for each value of mobility parameter, we produce 10000 stochastic simulations and for each simulation, we calculate $t_d$. The full range of $t_d$ values is then divided into 100 bins. Then we calculate the number of simulations for which the corresponding bin contains the observed $t_d$ and consequently, this gives the likelihood of mobility parameter given $t_d$. Finally, we take the logarithm of these likelihoods (i.e. log-likelihoods (LL)) and calculate the maximum likelihood estimate (MLE) of mobility parameter given $t_d$.

\begin{figure}[H]
    	\begin{center}
    	\includegraphics[width=1\textwidth]{timeseries_histd_td.png}
    	\end{center}
    	\caption{ (A) Example of stochastic realizations of the model for a fixed combination of model parameters. The red and blue trajectories refer to $I_{1}(t)$, and $I_{2}(t)$ respectively. The dashed red and blue line indicate the points at which $I_{1}(t)$ and $I_{2}(t)$ attain their maxima respectively. The difference between these time points is shown as $t_{d}$. (B) Histogram of $t_{d}$ for 10,000 stochastic realizations.    }
    	\label{fig:fig_schematic}
\end{figure}
\section{Results}\label{sec:results}
We first investigate the relationship between difference in peak timing ($t_{d}$) and cross-country mobility ($m_{ij}$). For the sake of simplicity we assume the mobility in both  directions is the same, i.e, $m_{12}$=$m_{21}$=$m$. We also take the recovery rates to be the same for both the regions, $\gamma_{1}=\gamma_{2}=\gamma$. We fix $\gamma=\frac{1}{14}~day^{-1}$ throughout our study unless specified otherwise. For the transmission rate, we consider two cases: (i) $\beta_{1}$ and $\beta_{2}$ are equal, i.e $\beta_1 =\beta_2=\beta$, and (ii) $\beta_1$ and $\beta_2$ are not equal. We consider the total population $N_{1}=10000$ and $N_{2}=20000$, which is fixed over time, and we assume initially there is no recovered individual in either region, i.e $R_{1}(0)=0$, and $R_{2}(0)=0$.

Following the method described in Section~\ref{subsec:stochastic_model}, for each value of the mobility parameter ($m$), we take $10,000$ realizations and calculate $t_d$ for different combinations of initial infected population ($I_{1}(0)$, and $I_{2}(0)$) and disease transmission rate ($\beta$). We see that the difference in peak timing ($t_{d}$) gradually decreases with increasing mobility flux $m$ (see Fig.~\ref{fig:fig_m_vs_td}(A)-(F)).

\begin{figure}[H]
    	\begin{center}
    	\includegraphics[width=1\textwidth]{m_vs_td.png}
    	\end{center}
    	\caption{(A)-(F)Behaviour of difference in peak-timing ($t_{d}$) with respect to coupling strength ($m$) for different combinations of transmission rate and initial infection. Here, the solid line represents the mean of the stochastic realizations and the shaded region represents $95\%$ confidence interval around the mean.}
    	\label{fig:fig_m_vs_td}
\end{figure}
We see that for a particular value of $m$, the difference in peak timing is higher for low disease transmission rate, and if we increase the transmission rate it starts to decrease. For instance, if we fix initial infecteds as $I_{1}(0)=1$, and $I_{2}(0)=10$, we see that for lower transmission rate $\beta=0.15$, and mobility parameter $m=0.025$, the difference in peak timing varies in the range [0,18] days (see Fig.~\ref{fig:fig_distri_td}(A)). If we slightly increase the transmission rate, i.e $\beta=0.2$, we observe that the range of $t_{d}$ shrinks to [0,10] days (see Fig.~\ref{fig:fig_distri_td}(B)). For higher transmission rate ($\beta=0.3$), the range of $t_d$, in this case, reduces to [0,6] days (see Fig.~\ref{fig:fig_distri_td})(C). Similar patterns are observed even if we change the number of initially infected individuals. For $I_{1}(0)=10$, and $I_{2}(0)=1$, we see, for $\beta=0.15$, the range of $t_{d}$ is [0,20] (see Fig.~\ref{fig:fig_distri_td}(D)). It shrinks to [0,14] for $\beta=0.2$ and finally to [1,9] for $\beta=0.2$ (see Fig.~\ref{fig:fig_distri_td}(E)-(F)). We also note that for low transmission rate values, the distribution of $t_d$ is flatter than that for higher transmission rates (see Fig.~\ref{fig:fig_distri_td}). This indicates that estimation of the mobility parameter from observed differences in peak timing ($t_d$) would be associated with greater uncertainty in the scenario where the transmissibility of the infection is comparatively lower.

\begin{figure}[H]
    	\begin{center}
    	\includegraphics[width=1\textwidth]{td_vs_density.png}
    	\end{center}
    	\caption{(A)-(F) Distribution of $t_d$ corresponding to the mobility parameter $m=0.025$, for different parameter combinations as mentioned in Fig~\ref{fig:fig_m_vs_td}.   }
    	\label{fig:fig_distri_td}
\end{figure}

Now we explore the feasibility of estimating the mobility flux ($m$) from observed differences in peak timing ($t_d$). We calculate the maximum likelihood estimate (MLE) of the parameter $m$, given observed value of $t_d$, and other model parameters ($\beta, \gamma,N_{1},N_{2}, I_{1}(0), I_{2}(0)$). Here, observed value of $t_d$ in each case is obtained by integrating the deterministic model~\eqref{EQ:Model} with same set of fixed parameters. Since it is evident from the above analysis that the transmission rate ($\beta$) is influential in calculating $t_{d}$, for an observed value $t_d$, we consider two different values of $\beta$ and estimate the MLE of the parameter $m$.

We calculate the log-likelihood (LL) of the parameter $m$ for a range of possible values using the stochastic simulation described in Section~\ref{subsec:stochastic_model}. Then using cubic-splines, we smooth the likelihood profile to obtain the maximum likelihood estimate and the corresponding confidence interval using a likelihood ratio test.
For example, we consider a case where the observed value of $t_{d}$ is $5$ days. If we fix $\beta=0.2$, then MLE of $m$ is $0.031$ and if $\beta$ is taken as $0.3$, the MLE of $m$ is reduced to $0.021$ (see Fig.~\ref{fig:mle_theory}(A), (C)). In another example, we consider higher observed value of $t_{d}$, i.e 10 days. In this case, for $\beta=0.2$ the MLE of $m$ is $0.018$ and for $\beta=0.3$ it becomes $0.01$ (see Fig.~\ref{fig:mle_theory}(B),(D)). From Fig.~\ref{fig:mle_theory}(A)-(D), in all four cases, we can see that MLEs of $m$ (vertical solid lines) are quite close to the true value of $m$ (dashed vertical line), the value of mobility parameter obtained from~\eqref{EQ:Model} for the corresponding value of $t_{d}$.

Note also that for a given $t_{d}$, the uncertainty in estimating the mobility parameter decreases with increasing disease transmission rate $\beta$. For both the lower and higher observed $t_{d}$ scenarios, the confidence intervals on $m$ shrink when we increase $\beta$ (see Fig.~\ref{fig:mle_theory}(A), (C) and (B), (D)). For a fixed value of $\beta$, if the observed value of $t_d$ is higher, then we also observe reduced uncertainty. For example, if $\beta$ is fixed as $0.2$, then the confidence interval on $m$ is wider for $t_d=5$ days than $t_d=10$ (see Fig.~\ref{fig:mle_theory}(A), (B)). A similar trend is also observed for higher transmission rates, e.g. $\beta=0.3$ (see Fig.~\ref{fig:mle_theory} (C), (D)).

\begin{figure}[H]
    	\begin{center}
    	\includegraphics[width=1\textwidth]{theory_mle.png}
    	\end{center}
    	\caption{ (A)-(D) Maximum likelihood estimates of the mobility parameter ($m$) given the difference in peak timing ($t_d$). The solid circles are the log-likelihood obtained from stochastic simulations and the solid curve is the smoothed likelihood profile. The solid vertical line is the maximum likelihood estimate of $m$ and the dashed vertical line is the true value of $m$, for which the deterministic model~\eqref{EQ:Model} gives the corresponding observed value of $t_{d}$. The shaded region depicts the $75\%$ confidence interval. The remaining parameters are fixed as $\gamma=\frac{1}{14}$, $N_{1}=10000$, $N_{2}=20000$, $I_{1}(0)=10$, $I_{2}(0)=1$.   }
    	\label{fig:mle_theory}
\end{figure}

Now we consider the case when $\beta_{1} \neq \beta_{2}$. We will explore how the estimate of the mobility parameter changes with the difference in transmission rates. To illustrate this, we consider a case where the observed value of $t_{d}$ is $5$ days. If we fix $\beta_1=0.25$, and $\beta_2 =0.20$ then MLE of $m$ is $0.076$ and if $\beta_1$, and $\beta_2$ are taken as $0.30$, and $0.20$ respectively, the MLE of $m$ is increased to $0.108$ (see Fig.~\ref{fig:mle_theory_diff_beta}(A), (C)). In another example, we consider a larger observed value of $t_{d}$, i.e 10 days. In this case, for $\beta_1=0.25$, and $\beta_2=0.20$, the MLE of $m$ is $0.038$ and for $\beta_1=0.30$, and $\beta_2 = 0.20$, it becomes $0.051$ (see Fig.~\ref{fig:mle_theory_diff_beta}(B),(D)). Similar to the case when $\beta_1=\beta_2$, from Fig.~\ref{fig:mle_theory_diff_beta}(A)-(D), we can see that MLEs of $m$ (vertical solid lines) are quite close to the true value of $m$ (dashed vertical line), the value of mobility parameter obtained from~\eqref{EQ:Model} for the corresponding value of $t_{d}$.

\begin{figure}[H]
    	\begin{center}
    	\includegraphics[width=1\textwidth]{theory_different_beta_mle.png}
    	\end{center}
    	\caption{ (A)-(D) Maximum likelihood estimate of the mobility parameter ($m$) given the difference in peak timing ($t_d$), when $\beta_{1}\neq \beta_2$. The solid circles are the log-likelihood obtained from stochastic simulation and the solid curve is the smoothed likelihood profile. Solid vertical line is the maximum likelihood estimate of $m$ and the dashed vertical line is the true value of $m$, for which the deterministic model~\eqref{EQ:Model} gives the corresponding observed value of $t_{d}$. The shaded region depicts the $75\%$ confidence interval. The remaining parameters are fixed as $\gamma=\frac{1}{14}$, $N_{1}=10000$, $N_{2}=20000$, $I_{1}(0)=10$, $I_{2}(0)=1$.   }
    	\label{fig:mle_theory_diff_beta}
\end{figure}

\subsection*{Case study: COVID-19 incidence in Poland-Germany border region}
We apply our proposed method to COVID-19 incidence data during the period $20^{\rm th}$ February 2021 to $20^{\rm th}$ June 2021, for the pair of regions which are located on either side of the Poland-Germany border and are directly connected through different modes of transportation e.g., bus-transportation, train-transportation (see~Fig.\ref{fig:mle_data} (A)). For each pair of regions, we obtain the corresponding information on total population, and initial number of infections (see Table~\ref{table:table_case_study}). All these datasets are obtained from the database pipeline by~\cite{abdussalam2022scalable}. We fix the value of average infectious period ($\frac{1}{\gamma}$) at $14$, following ~\cite{byrne2020inferred}.

We have previously demonstrated that the disease transmission rate plays an influential role in the estimation of mobility flux given the difference in peak-timing. Therefore, in the case of real data, instead of fixing the disease transmission rate arbitrarily, we follow a simpler approach to estimate these value from the initial phase of the epidemic. Since, early in the epidemic the daily number of new cases grows exponentially over time, for a SIR-type model the number of infection at time $t$ can be approximated as~\citep{ma2020estimating}:
\begin{equation*}
    I(t)\approx I(0)\exp{(\lambda t)}.
\end{equation*}
Here, $\lambda$ is the growth rate and is given by $\lambda=\beta - \gamma$. Since cumulative number of new cases per day ($C(t)$) and new cases per day ($I(t)$) are linearly related~\citep{favier2006early, pinho2010modelling}, i.e $I(t)$ $\sim$ $\lambda$ $C(t)$, we estimate the growth rate $\lambda$ using a linear fit. Therefore, once the parameter $\gamma$ is kept fixed, the disease transmission rate $\beta$ can be approximated as $\beta \approx \lambda + \gamma$. Following this approach, for each region we estimate the disease transmission rate from the time series data of daily new cases (see Table~\ref{table:table_case_study}).

 \begin{table}[H]
  \tiny

  	\centering
  	\begin{tabular}{|p{3.1cm}|p{1.5cm}|p{1.5cm}|p{1.3cm}|p{1.3cm}|p{1.9cm}|}

  		\hline
  		\centering{\textbf{Pair (Powiat, Kreise)}} &\centering{\textbf{Transmission rate ($\beta_{1}$, $\beta_2$)} } &\centering{\textbf{Total Population ($N_{1}$, $N_{2}$)}} & \textbf{Initial infection $(I_{1}(0), I_{2}(0))$} & \textbf{Observed difference in peak-timing ($t_{d}$)} & \textbf{Estimated mobility flux ($m[75\% CI]$)} \\
  		[2ex]
  		\hline
  		\centering (Gryfi\'{n}ski, Uckermark) & (0.10,0.075) & (82951,119552)  & (11,8)  & 20 days & 0.038[0.022,0.064]\\[5ex]
  		\centering (S\l{}ubicki, Frankfurt(Oder)) & (0.13,0.14)  & (47068,57873) & (8,1) & 10 days & 0.032[0.004,0.052]\\[5ex]
  		\centering(Zgorzelecki,G\"{o}rlitz) & (0.10,0.11) & (90584,254894) & (9,23) & 20 days  & 0.012[0.002,0.018]\\[4ex]

  		\hline
  	\end{tabular}
  	\caption{Parameters and observed difference in peak-timing ($t_{d}$) for each pair of region used for estimating the mobility flux $m$. The estimated value of $m$ is presented in the last column of the table. }
  	\label{table:table_case_study}
  \end{table}

We first consider the regions: Gryfi\'{n}ski, and Uckermark, which are located at the upper part of the Poland-Germany border. The observed difference in peak-timing obtained from the real data is $20$ days (see Fig.~\ref{fig:mle_data}(B)). Given all the other model parameters, we estimate the mobility flux as $0.038 [0.022,0.064]$ (see Fig.~\ref{fig:mle_data}(C)). Next, we move to the middle part of the border and choose the regions: S\l{}ubicki, and Frankfurt(Oder). In this case, the observed difference in peak-timing is $10$ days (see Fig.~\ref{fig:mle_data} (D)). The maximum likelihood estimate of $m$ is $0.032[0.004,0.052]$ (see Fig.~\ref{fig:mle_data} (E)). Now for the regions: Zgorzelecki, and G\"{o}rlitz, the observed difference in peak timing is $20$ days (see Fig.~\ref{fig:mle_data} (F)). The estimated value of the parameter $m$ in this case is $0.012[0.002,0.018]$ (see Fig.~\ref{fig:mle_data} (G)).

Comparing the first case and third case, we see that there is a substantial difference between the estimates of mobility parameters even though the observed differences in peak-time ($t_{d}$) for both cases are same. Instead, the difference in estimated transmission rates (i.e, $\beta_{1}$ and $\beta_2$) is higher in the first case (Gryfi\'{n}ski, and Uckermark) than that of the third case (Zgorzelecki, and G\"{o}rlitz) (see Table~\ref{table:table_case_study}). This marked difference in transmission rates likely explains the different estimates of the coupling parameter.

\begin{figure}[H]
    	\begin{center}
    	\includegraphics[width=1\textwidth]{mle_border_data.png}
    	\end{center}
    	\caption{(A) The map shows the three pairs of regions that share a border between Poland and Germany. The Polish regions (called \textit{Powiat}) are filled with dark blue color, and the German regions (called \textit{Kreise}) are shaded dark green. The name of the regions are: top-(Gryfi\'{n}ski, Uckermark); middle- (S\l{}ubicki, Frankfurt(Oder)); and bottom-(Zgorzelecki,G\"{o}rlitz). (B), (D), (F) The time series of number of infecteds in the selected regions. The blue lines represent the time series for the Polish regions and green lines represent the time series for the German regions. Here, time zero refers to the date $20^{\rm th}$ February 2021 and time 120 refers to the date $20^{\rm th}$ June 2021. (C), (E), (G) The likelihood profiles of the cross-border mobility $m$. The gray circles are the log-likelihoods obtained from the stochastic simulation and the solid curve is the smoothed likelihood profile. The dashed vertical line is the maximum likelihood estimate of $m$ and the shaded region denotes the $75\%$ confidence interval.  }
    	\label{fig:mle_data}
\end{figure}

\section{Discussion}\label{sec:discussion}
In general, understanding the interplay between human mobility and spatio-temporal propagation of infectious diseases is a long-standing problem in the field of infectious disease modelling. The unavailability of reliable information on human mobility poses challenges to designing effective public health strategies. In this study, we proposed a method to address the issue of estimating cross-country mobility flux during a disease outbreak, which is important yet less explored in the literature. We showed that from the observed  difference in peak-timing between any two regions which share a border, the underlying mobility flow between them can be retrieved effectively. Since, the information on peak timing is less affected than the invasion time by the delay or incompleteness in reporting infection, therefore the difference in peak-timing can be considered as a more robust quantity than invasion time. Consequently, our approach would be more broadly applicable for estimating mobility flux in real-world case studies than existing alternatives that rely on the observed disease invasion time~\citep{hempel2020estimating}.

 We investigated how different model parameters influence the difference in peak-timing in a simple scenario, where the disease transmission rates in both regions are the same. We see that the disease transmission rate plays key role in determining the distribution of difference in peak-timing. The difference in peak-timing decreases with the disease transmission rate (Fig.~\ref{fig:fig_m_vs_td}). This is because of the fact that higher disease transmission rate causes earlier occurence of the peak in infection trajectory and thus the difference in peak-timing between the two regions also decreases. It has also been revealed that, for lower transmission rates, the distribution of differences in peak-timing becomes flatter and can have a wider range of possible values (Fig.~\ref{fig:fig_distri_td}(A), (D)). On the other hand, with increasing transmission rates the distribution progressively narrows (Fig.~\ref{fig:fig_distri_td} (B), (C), (E), (F)).

Based on a simple SIR-type model, with the help of stochastic simulations, we provided maximum likelihood estimate of the cross-border mobility under different epidemiologically relevant scenarios. From our numerical investigation, it is evident that the estimation method based on the observed difference in peak-timing can effectively recover the underlying mobility rate for both scenarios: (i) when the disease transmission rates in both regions are assumed to be the same (Fig.~\ref{fig:mle_theory}(A)-(D)), and (ii) when the transmission rates are assumed to be different (Fig.~\ref{fig:mle_theory_diff_beta}(A)-(D)). In the first scenario, for a given difference in peak-timing, the uncertainty of the estimate decreases with the transmission rate (Fig.~\ref{fig:mle_theory}(A), (C) and Fig.~\ref{fig:mle_theory}(B), (D)). Also for a fixed transmission rate, the uncertainty in estimation decreases with increasing values of the observed difference in peak timing (Fig.~\ref{fig:mle_theory}(A), (B) and Fig.~\ref{fig:mle_theory}(C), (D)). In the second scenario, for a given observed difference in peak-timing, the estimate of mobility parameter increases with the difference between the disease transmission rates (Fig.~\ref{fig:mle_theory_diff_beta}(A), (C) and Fig.~\ref{fig:mle_theory_diff_beta}(B), (D)). Similar to the first scenario, for fixed disease transmission rates, the uncertainty in the estimate decreases with the observed difference in peak-timing (Fig.~\ref{fig:mle_theory_diff_beta}(A), (B) and Fig.~\ref{fig:mle_theory_diff_beta}(C), (D)).

To demonstrate its real-world utility, we applied our proposed method to COVID-19 incidence data obtained for the regions located along the border of Poland and Germany. We chose only three pairs of regions located in different parts of the Polish-German border which are well connected through different modes of transportation. However, the method can also be applied to any other pair of regions that shares border. Using least squares techniques, we first estimated the disease transmission rate for each of the regions from the early phase of epidemics. We then calculated the maximum likelihood estimate of the mobility parameter for three pairs of regions located in different areas along the border. It is noteworthy that, if we compare the pair of regions Gryfi\'{n}ski, Uckermark and Zgorzelecki, G\"{o}rlitz, we see that, even though the differences in peak-timing are the same, the estimated values for the cross-border mobility are very different. Since the difference in transmission rates between the regions Gryfi\'{n}ski and Uckermark is higher than that of Zgorzelecki and G\"{o}rlitz, we see higher mobility flow for the regions Gryfi\'{n}ski and Uckermark. This result is in good agreement with our theoretical findings (Fig.~\ref{fig:mle_theory_diff_beta}).

Even though our study addressed the important issue of estimation of cross-border mobility and our proposed method worked well in the real case study, it has some limitations. We assumed that the mobility between two regions is symmetric, which is not true in general, specifically if we consider the mobility between a big city and a village. We considered the regions on either side of the border as isolated and ignored the contribution of the adjacent regions which is not realistic. We leave these modifications and extensions for future study.

Summing up, this work provides a simple and intuitive method for estimating mobility flow between the regions connected through human migration. The methodology and findings can serve as a basis for designing intervention strategies in situations where direct information on human mobility is inaccessible.

\section*{Acknowledgements}
This work was partially funded by the Center of Advanced Systems Understanding (CASUS), which is financed by Germany’s Federal Ministry of Education and Research (BMBF) and by the Saxon Ministry for Science, Culture, and Tourism (SMWK) with tax funds on the basis of the budget approved by the Saxon State Parliament.

\bibliographystyle{elsarticle-harv}