\section{Introduction}
Vehicle-to-vehicle (V2V) communication is seen as one of enabling technologies for intelligent transportation systems and a key enabler for many smart road and traffic management systems \cite{Geo}, as it allows critical information dissemination. Moreover, spurred by the availability of in-vehicle infotainment (IVI) systems disseminating entertainment and information content to passengers \cite{YWu}, there is a strong need for a high-speed and stable delivery of large multimedia files, such as videos, photos and songs, to the various cars within a vehicular network. For effective dissemination of such diverse content across vehicular networks, there is a need for effective content placement strategies to maximize the throughput of the system \cite{DR}. Moreover, to reap the benefits of V2V content dissemination, we must address different challenges including optimal cache placement and resource allocation \cite{MM}.\setlength{\parskip}{0pt}

Cache placement in V2V networks has recently attracted significant attention such as in \cite{Haowu}-\cite{Shosny2}.
In such works, the popular contents are offloaded to the storage of a number of well-chosen cars and devices at off-peak hours in order to serve requests during peak traffic hours. In these scenarios, cars and devices that do not have the cached content will not have to download the content from wireless base stations (BSs). Instead, they can request the content directly from other cars having cached data which can eventually lead to a reduction in the traffic load at the BSs. Meanwhile, caching using local storage can reduce latency due to a shorter communication distance. However, the benefits of caching are highly dependent on the set of cars and devices chosen for caching the popular contents \cite{Smitha}.

In \cite{Haowu}, the authors use vehicle mobility data for content dissemination and combine the idea of opportunistic forwarding, trajectory based forwarding and geographical forwarding to develop a mobility-centric algorithm to place content in vehicular networks.
Meanwhile, the work in \cite{Guan} applies available users mobility patterns to develop a polynomial-time solution to maximize the saved cost by caching contents in local storage.
Moreover, in \cite{Smitha}, the authors use location information and subscription-based information to divide vehicles into different groups and, then, design a spatio-temporal multicast routing protocol to construct an optimized dissemination mesh network. The work in \cite{Shosny2} presents an optimal caching policy by using both mobility information and users' demands and proposes a greedy caching algorithm with polynomial order complexity to obtain bounds of the caching policy whose complexity grows exponentially with the number of users.
However, existing works, including \cite{Haowu}-\cite{Shosny2}, do not take into account the temporal dynamics of V2V communication networks, such as the frequency of occurrence of different V2V links, which can be a key metric for content placement. For example, using only location information to select the set of cars that will cache the content can lead to choosing cars which are unable to communicate with each other or with other cars, effectively limiting the benefits of caching.\setlength{\parskip}{0pt}

The main contribution of this paper is a novel framework for spatio-temporal caching in vehicular networks that is cognizant of intrinsic spatial and temporal patterns in V2V communications. In particular,
using the tools of \emph{temporal network analysis} \cite{UAlon}, we identify temporal motifs in the V2V network as key communication patterns observed among vehicles that appear more frequently compared with what is expected in a baseline, randomized reference system. After identifying the spatial-temporal motifs, the proposed approach then finds the best candidate cars for content placement. To our best knowledge, this is \emph{the first work that exploits spatial-temporal motifs for optimizing content dissemination in vehicular networks}. Simulation results using real traces from the city of Cologne, Germany, show that the proposed motif-based approach yields significant performance gains in terms of the average data rate, compared to a conventional location-based scheme.

The rest of the paper is organized as follows. Section \uppercase\expandafter{\romannumeral2} presents the system model and problem formulation. In Section  \uppercase\expandafter{\romannumeral3}, we present the proposed motif-based approach. Section \uppercase\expandafter{\romannumeral4} provides the simulation results. Conclusions are drawn in Section \uppercase\expandafter{\romannumeral5}.\vspace{-0.1cm}

\section{System Model and Problem Formulation}\vspace{-0.1cm}
\begin{figure}[!t]
\centering
\includegraphics[width=3.5in,height=1.8in]{System_model.pdf}\vspace{-0.4cm}
\DeclareGraphicsExtensions.
\caption{V2V communication networks and its macroscopic communication graph. \vspace{-0.6cm}}\vspace{-0.4cm}
\label{systemmodel}
\end{figure}

Consider a vehicular network in an urban environment, composed of a set $\mathcal{N}$ of $N$ vehicles. In particular, as shown
in Fig.\ref{systemmodel}, we consider that a nearby BS $b$ provides wireless coverage to the cars in $\mathcal{N}$. The BS frequently uses beacon signals to keep track of the vehicles' positions by using the received signal-strength measurements in \cite{Martin}, and collects the V2V communication data when exchanging information with the communication facility within the vehicles as proposed in \cite{Hubaux}.

In our system, the BS seeks to seed multimedia files and traffic data in the storage unit of a set $\mathcal{C}\!\!\subseteq\!\!\mathcal{N}$ of $C$ cars to provide the passengers in neighboring cars with streaming services or to disseminate delay-sensitive information, such as upcoming road incidents to help drivers decide on their routes. Using this mechanism, the seeded cars, refereed to hereinafter as \emph{serving cars}, can disseminate the cached content to nearby cars via V2V links. Consequently, the BS can reduce its traffic load, as it no longer needs to transmit the same content to multiple cars. Fig. \ref{systemmodel} shows a case in which the BS caches content at car $c$. Here, an arbitrary car $a \!\!\in \!\!\mathcal{N}\!\setminus\!\! \{c\} $ sends a request for content $m$ to the BS $b$.
If the content is already cached at car $c$, the BS would inform the requesting car $a$ to use V2V communication to obtain the content from $c$.
Otherwise, the BS has to directly transmit the requested content to $a$, which increases its load. Hereinafter, we refer to cars that do not have the cached content as \emph{non-serving cars}. To increase the spectrum efficiency, we assume that V2V links reuse the spectral resources of the cellular network as in \cite{Yu}.\vspace{-0.1cm}

\subsection{Problem formulation}
The signal-to-interference-plus-noise ratio (SINR) of a V2V communication link between a serving car $c$ and a non-serving car $a$ will be:
\begin{align}\label{1}
\gamma_{c,a}=\frac{P_{c,a}g_{c,a}}{I_{b,a}+I_{a}^{\prime}+\sigma^{2}},
\end{align}
where $P_{c,a}$ is the transmission power from $c$ to $a$, and $\sigma^{2}$ is the variance of the Gaussian noise at the receiver. In addition, $g_{c,a}$ represents the channel gain between $c$ and $a$ and can expressed as $g_{c,a} \!= \! \eta d_{c,a}^{-\alpha}$, where $\eta$ is the fading gain, $d_{c,a}$ is the distance, and $\alpha$ is the path loss exponent.
We consider that the cellular links and the V2V links will experience a Rayleigh fading as done in \cite{Fang}.
$I_{b,a}$ and $I_{a}^{\prime}$ capture, respectively, the interference generated by the links between the BS and other vehicles and by other V2V communication links.
These interference terms are given by,
\begin{align}\label{2}
I_{b,a}\!\!= \!\!\sum_{a^{\prime}\in \mathcal{N}\setminus a}\!\!(L_{b,a^{\prime}}) P_{b,a^{\prime}}g_{b,a},
I_{a}^{\prime}\!\!=\!\!\sum_{c_{1}, a_{1}\in \mathcal{N}\setminus a}\!\! (L_{c_{1},a_{1}})  P_{c_{1},a_{1}}g_{c_{1},a}\,
\end{align}
where $L_{i,j}$ is a binary variable that captures the feasibility of  communication link between $i$ and $j$. In fact, $L_{i,j}=1$ if the SINR at the receiving car $j$ from the transmitting car $i$ exceeds a target threshold for V2V communication,
otherwise, $L_{i,j}=0$.

Accordingly, the achievable data rate for the V2V link between the serving car $c$ and non-serving car $a$ is
\begin{align}\label{5}
R_{c,a}=\omega\log_{2}(1+\gamma_{c,a}),
\end{align}
where $\omega$ is the bandwidth. The corresponding SINR $\gamma_{b,a}$ and the achievable rate $R_{b,a}$ between the BS and non-serving car $a$, are:
\begin{align}\label{51}
\gamma_{b,a}= \frac{P_{b,a}g_{b,a}}{I_{a}^{\prime}+\sigma^{2}},
\vspace{1in}R_{b,a}=\omega\log_{2}(1+\gamma_{b,a}).
\end{align}

We consider a set $\mathcal{M}$ of $M$ popular contents with the same size, sorted from high to low popularity. Due to the limited storage, the serving nodes will choose to cache a limited number of contents. Therefore, to increase the possibility of meeting the requirements from non-serving nodes, the serving cars will cache $F\!\! \leq \!\!M$ most popular files, under the capacity constraints.
Therefore, when a non-serving node requests one file out of the $F$ files, the node can acquire it via a V2V link as long as there exists one serving node within its communication range.
To leverage V2V communications for disseminating cached content, the network must determine which vehicles act as serving nodes so as to maximize the average data rate achieved by non-serving cars. To this end, we define the binary variable $x_i$ for each car $i$, where $x_{i}\!\!=\!\!1$ if car $i$ is selected as serving node, otherwise $x_{i}\!\!=\!\!0$, and set $\mathcal{C}$ contains $C\!\!=\!\!\sum_{i}^{N}\!\!x_{i}$ serving nodes. Therefore, we can formulate the problem as follows:
\begin{align} \label{11}
 &\max_{\mathcal{C} \subseteq \mathcal{N}}\!\!\frac{1}{(N\!\!-\!\! C)}\! \!\sum_{a \in \mathcal{N}\setminus \mathcal{C}}\! \! \left(\sum_{m=1}^{F}\text{Pr}_{a}(m) \sum_{c\in \mathcal{C}}\beta_{c,a}R_{c,a}\!+\!\! \!\! \!\sum_{m=F+1}^{M}\! \!\! \!\text{Pr}_{a}(m)R_{b,a}\right)      \\ \vspace{-0.4cm}
&\text{s.t.}\hspace{0.1in}
P_{c,a}\leq P_{\text{max}}, \label{size} \\
&\hspace{0.3in}
\bar{\gamma}\leq\gamma_{c,a}, \label{1111}
\hspace{0.4in}
\end{align}
where $\text{Pr}_{a}(m)$ is the probability mass function (pmf) of the request for file $m$ by car $a$.  This distribution can be modeled by the Zipf distribution with
pmf $\text{Pr}_{a}(m)\!\! =\!\! \frac{1}{m^{\theta_{r}}}/\sum_{x=1}^{M}\frac{1}{x^{\theta_{r}}}$, where $\theta_{r}$ is the Zipf exponent that determines the skewness of the
distribution \cite{Malak}. The indicator variable $\beta_{c,a}$ is such that $\beta_{c,a}\!\!=\!\!1$ if the serving car $c$ is the nearest car to the requester $a$ in the serving set $\mathcal{C}$; $\beta_{c,a}\!\! =\!\!0$, otherwise. $\beta_{c,a}$ ensures that a non-serving car will always choose the closest serving car $c$ with the requested file.
Constraint (\ref{size}) guarantees that the transmission power of vehicles will not surpass the maximum power level $P_{\text{max}}$, and
constraint (\ref{1111}) ensures that SINR of the V2V links is above a threshold, $\bar{\gamma}$. \vspace{-0.1cm}

\subsection{V2V macroscopic communication graphs}

As long as we find the optimal set of serving vehicles, we can assign non-serving nodes to serving nodes according to the spatial distance, and solve the problem given by (\ref{11})-(\ref{1111}). However, finding the optimal set $\mathcal{C}$
is a 0-1 integer programming where determining whether each individual node in the set $\mathcal{N}$ should be considered as either a serving or a non-serving node. In fact, the problem is one of Karp’s 21 NP-complete problems \cite{Karp}, which is hard to solve directly.

Alternatively, we can find a sub-optimal solution to the optimization problem by using the information within the vehicular network. In particular, in addition to the the position and demand of each vehicle, we can also leverage the temporal domain information.
This is because V2V networks are naturally dynamic and exhibit the temporal features. For example, the number of communication links between two arbitrary cars may vary with time. Such temporal information on the frequency of communication is valuable to determine which cars are more likely to better disseminate the content. Therefore, time domain information is also important to choose the optimal set $\mathcal{C}$ for cache placement and solve the problem (\ref{11})-(\ref{1111}).

To capture the dynamics in the time domain, we propose to use collected V2V communications data and model the system as a directed \emph{temporal graph} $G \!= \!(\mathcal{N}, \mathcal{E})$, whose vertices are the cars in $\mathcal{N}$ and whose temporal edges, in set $\mathcal{E}$, denote the time-stamped communication events among vehicles. We represent a wireless communication link between two different cars
as a \emph{3-tuple} edge labeled as $i$, $e_{i}\!=\!<\!c_{i},a_{i},t_{i}\!>\!$, where $c_{i}$ and $a_{i}$ denote the serving car and the non-serving car, respectively. $t_{i}$ is the time of initiating the transmission from $c_{i}$ to $a_{i}$. Here, we assume that, at any time, a car cannot communicate with more than one car simultaneously. In this case, there are no two edges with the common vehicle element initiated at the same time.

To obtain the set of cars that are more active and more likely to participate in the V2V communication in a period of time, we introduce the time constraint $T$ and devide the graph $G \!= \!(\mathcal{N}, \mathcal{E})$ into multiple \emph{macroscopic communication graphs}, where we can change the value of $T$ to filter outdated V2V links. That is, for any edge $e_{i}$ in these macroscopic graphs, there always exists at least one other edge $e_{j}, i\neq j$, in the same graph, that meets the following requirements:
\begin{itemize}
  \item Two edges share at least one node, $(c_{i},a_{i}) \bigcap (c_{j},a_{j}) \neq \emptyset$;
  \item If a wireless connection $i$ occurs before another connection $j$, $0 \leq t_{j}-t_{i} \leq T$;
  \item If a wireless connection $j$ occurs before another connection $i$, $0 \leq t_{i}-t_{j} \leq T$.
\end{itemize}

Given these definitions, using the graph-theoretic framework of \emph{temporal networks analysis}, next, we present a novel graph-theoretic sub-optimal solution to the original problem in (\ref{11})-(\ref{1111}).\vspace{-0.1cm}

\section{Proposed Strategy Based on Spatio-Temporal Motifs}\vspace{-0.1cm}
\begin{algorithm}
\caption{ Microscopic Subgraphs Searching Algorithm}\label{euclid}
\textbf{Input:} A V2V macroscopic communication graph $G_{1}=(\mathcal{N}_{1},\mathcal{E}_{1})$, and target subgraph size $k$. \\
\textbf{Output:} The set of microscopic subgraph $\mathcal{V}_{1}$ with size $k$.
\begin{algorithmic}[1]\footnotesize
\State  $\mathcal{V}_{1}=\varnothing$
\State \textbf{for} {arbitrary edge $e_{i} \in \mathcal{E}$ } \textbf{do}
\State \hspace{10pt} $\mathcal{V}_{1}\leftarrow \{e_{i}\}$
\State \hspace{10pt} \textbf{for} edge $e_{j}\in \mathcal{E}-e_{i}$  \textbf{do}
\State \hspace{20pt} \textbf{if} $\{c_{i},a_{i}\}\cap\{c_{j},a_{j}\}\neq \varnothing$ \textbf{and} $j>i$
\State \hspace{30pt} $\mathcal{V}_{2}\leftarrow \{e_{j}\}$
\State \hspace{20pt} \textbf{end}
\State \hspace{10pt} \textbf{end}
\State \hspace{10pt} \textbf{call} \textproc{EdgeExtension}($\mathcal{V}_{1}, \mathcal{V}_{2}, k$)
\State \textbf{return} $\mathcal{V}_{1}$
\State \textbf{function} \textproc{EdgeExtension}($\mathcal{E}_{\text{sub}}$, $\mathcal{V}_{3}$, $k$)
\State  \textbf{If} $|\mathcal{E}_{\text{sub}}|= k$ \textbf{then output} $G[\mathcal{E}_{\text{sub}}]$ and \textbf{return}
\State \hspace{0pt} \textbf{while} $\mathcal{V}_{3}\neq \emptyset$ \textbf{do}
\State \hspace{10pt} \textbf{for} arbitrary edge $e_{x} \in \mathcal{V}_{3}$  \textbf{do}
\State \hspace{20pt} Remove $e_{x}$ from $\mathcal{V}_{3}$
\State \hspace{20pt} \textbf{for} edge $e_{y}\in \mathcal{V}_{3}-e_{x}$  \textbf{do}
\State \hspace{30pt} \textbf{if} $\{c_{x},a_{x}\}\cap\{c_{y},a_{y}\}\neq \varnothing$ \textbf{and} $y>x$
\State \hspace{40pt} $\mathcal{V}_{4}\leftarrow \{e_{y}\}$
\State \hspace{30pt} \textbf{end}
\State \hspace{20pt}\ \textbf{end}
\State \hspace{20pt} $\mathcal{V}_{3}^{\prime}=\mathcal{V}_{3} \cup \mathcal{V}_{4}$
\State \hspace{20pt} \textbf{call} \textproc{EdgeExtension}($\mathcal{E}_{\text{sub}}\cup \{e_{x}\}$, $\mathcal{V}^{\prime}_{3}$, $k$)
\State \hspace{10pt}\ \textbf{end}
\State  \textbf{return}
\newline
\vspace{-0.4cm}
\end{algorithmic}
\end{algorithm}
\setlength{\parskip}{0pt}

To find a sub-optimal solution to (\ref{11})-(\ref{1111}), we need to develop an algorithm to detect motifs in V2V communication macroscopic graphs. To this end, we follow three steps. The first step is searching \emph{microscopic subgraphs} defined as basic communication units with typical sizes (number of edges) in macroscopic communication graphs. The second step is collecting microscopic subgraphs with similar isomorphic structure.
The third step is determining the frequently occurring motifs after comparing each subgraph's frequency $f$ of occurrence with its counterpart $f^{\prime}$ in a baseline randomized V2V communication network. In particular, we use the notion of \emph{Z-score}, expressed as \vspace{-0.1cm}
\begin{align}\label{zscore}
Z=\frac{f-f^{\prime}}{\sigma_{m}},
\end{align}
where $\sigma_{m}$ captures the standard deviation of the corresponding subgraph in the reference system. If $Z\!\! >\!\! Z_{th}$, where $Z_{th}$ is a given threshold, we can classify the subgraph as a motif \cite{Temporal}. Finally,
given the detected motifs, our spatio-temporal caching strategy can select $\mathcal{C}$, $\mathcal{C}\!\! \subseteq \!\!\mathcal{N}$, to solve the optimization problem in (\ref{11})-(\ref{1111}). \vspace{-0.4cm}
\begin{figure*}[!t]
\centering
\includegraphics[width=5.5in]{ValidSubgraph.pdf}\vspace{-0.3cm}
\DeclareGraphicsExtensions.
\caption{Considering a graph with $4$ nodes and $6$ edges in the corner, we can obtain all microscopic subgraphs with size as $3$ after following Algorithm 1.}
\vspace{-0.7cm}
\label{systemmodel111a}
\end{figure*}
\vspace{-0.2cm}

\subsection{Searching for V2V microscopic subgraphs}
To effectively detect existing motifs in V2V communication networks, we first decompose V2V macroscopic communication graphs into microscopic subgraphs with the same size.

The algorithm used for finding V2V microscopic subgraphs with target size $k$ in a macroscopic graph $G_{1} \!= \!(\mathcal{N}_{1}, \mathcal{E}_{1})$ is shown in Algorithm \ref{euclid}. Given a set of labeled cars and labeled connection edges, there are two ways to obtain microscopic subgraphs from the macroscopic graph. One approach is to consider the nodes connected by edges, and another approach is to consider the set of edges in which an arbitrary edge can find another edge sharing the common node. In contrast to the work in \cite{SW}, the proposed algorithm is based on the second approach. As shown in Fig. \ref{systemmodel111a}, from the macroscopic graph, we can obtain a pair of edge sets, $\mathcal{V}_{1}$ and $\mathcal{V}_{2}$, where the first set contains only one edge $e_{i}$, labeled as $i$, and the second set has $e_{i}$'s all neighboring edges with greater labels.
The mechanism will first add an arbitrary edge $e_{j}$, labeled as $j$, from the second set to the first set and update the first set. By calling \textproc{EdgeExtension($\mathcal{E}_{\text{sub}}$, $\mathcal{V}_{3}$, $k$)}, the second set could be extended by first merging the set of edges, which are neighbor to the newly added edges $e_{j}$ and with greater labels in the macroscopic graph, and then removing edge $e_{j}$ and other edges with smaller labels compared with $j$.
Next, we repeat the aforementioned steps for the first set with the added edge and the new second set until we obtain the microscopic subgraphs meeting the size requirement, i.e., $k$.

\textbf{Example:} When the two sets of edges from the macroscopic graph are \{1\} and \{2,3,4,6\}, as shown in Fig. \ref{systemmodel111a}, the algorithm first adds edge 2 to the first set \{1\} and obtains the first updated set \{1,2\}. Then, it updates the second set as \{3,4,5,6\} by first merging the neighboring set with greater labels, i.e., \{3,5,6\} with \{2,3,4,6\} and then deleting the edge 2. Repeatedly, the algorithm follows the same processes for the updated set \{1,2\} and the edge set \{3,4,5,6\}, and we can obtain the microscopic subgraphs \{1,2,3\} and \{1,2,6\}, if $k=3$.

As shown in Fig. \ref{systemmodel111a}, by using this algorithm for any one edge set and the set of its surrounding edges with greater labels, we can finally collect all V2V microscopic subgraphs with the required size in a macroscopic graph. To obtain the whole set of microscopic subgraphs existing in the V2V network, we can use this algorithm for different macroscopic graphs. Although the complexity of the algorithm will increase with the increment of the vehicles, we can apply the algorithm in scenarios with capacity limitations, like intersections,  parking lots, and parts of the highway, or a subset of the network, thus reducing the complexity and the processing time. \vspace{-0.1cm}

\subsection{Classifying V2V microscopic subgraphs}

\begin{figure}[!t]
\centering
\includegraphics[width=2.8in,height=1.8in]{CanonicalLabeling.pdf}\vspace{-0.4cm}
\DeclareGraphicsExtensions.
\caption{Examples on the canonical labeling of a microscopic subgraph and other microscopic subgraphs sharing the same canonical labeling.}
\label{systemmodel2}\vspace{-0.6cm}
\end{figure}

To sort V2V microscopic subgraphs with the same isomorphic structure, i.e., to find microscopic subgraphs which contain the same number of cars connected in the same way, we exploit the notion of \emph{canonical labeling}. By permuting the elements in the adjacency matrices obtained from a microscopic subgraph, we could construct many lists of integers. By viewing each list of integers as a string of 1s and 0s, we sort them based on lexicographic ordering to obtain a \emph{canonical labeling} defined as the string with the minimum value. Due to the fact that the canonical labeling of two subgraphs will be identical as long as they have the similar isomorphic structure \cite{R.C}, the problem of
determining isomorphic structures among subgraphs is equivalent to deciding whether given microscopic subgraphs
have the same canonical labeling or not.

For example, as shown in Fig. \ref{systemmodel2}, by concatenating rows or columns one after the other in the permuted adjacency matrices of the first microscopic extracted from Fig. \ref{systemmodel111a}, we can find two strings of 1s and 0s.
According to lexicographic ordering, we notice the string corresponding to the second matrix is smaller than its counterpart in the first matrix, i.e., ``000010101''<``0001011'', and ``0001011'' can be chosen as the canonical labeling for the first microscopic subgraph compared with other permutations. Similarly, we can observe the second microscopic subgraph shares the same canonical labeling with the first subgraph, and, thus, these two subgraphs have the same isomorphic structure.

Accordingly, we are capable of completing structure classification for all the microscopic subgraphs in the V2V communication network, and then, we can calculate the occurrence frequency $f$ of microscopic subgraphs having a similar isomorphic structure.
Moreover, we also repeat the same steps for a randomized V2V communication network so as to obtain the mean frequency $f^{\prime}$ and the standard deviation $\sigma_{m}$ of the corresponding subgraph. We can use (\ref{zscore}) to determine whether the subgraph is motif or not.  \vspace{-0.25cm}

\subsection{Proposed spatio-temporal caching strategy}
\begin{figure}[!t]
\centering
\includegraphics[width=3in,height=2.6in]{final2.pdf}
\DeclareGraphicsExtensions.\vspace{-0.2cm}
\caption{Flow chat of the proposed spatio-temporal caching strategy.}
\label{Pro2}\vspace{-0.5cm}
\end{figure}
To obtain the sub-optimal solution for the optimization problem (\ref{11})-(\ref{1111}), as shown in Fig. \ref{Pro2}, the BS will first select the serving cars in the set $\mathcal{C}$ based on the motifs detected from the V2V communication graph. In particular, we assume there are $W$ motifs and the Z-score of a motif $i$ is $Z_{i}$. Note that for any $i\in \{1,...,W\}$, $Z_{i}>Z_{th}$. By observing the motif structure, we can determine the outdegree of each car, which is defined as the number of outgoing edges emanating from the car. For example, for the first microscopic subgraph in Fig. \ref{systemmodel2}, the outdegrees of $a, b, c$ are 1, 2, and 0. Then, we choose the influential car as the node with the maximum outdegree in the corresponding motif, since the connections that originate from this car can reach more recipients compared with other cars, leading to a more effective content dissemination. Next, we can statistically acquire the frequency of car $j$'s being the influential car in the $i-$th motif as $f_{ij}, j \in \mathcal{N}$. Based on that, as shown in Fig. \ref{Pro2}, we can obtain the sum frequency of being the influential car in the motifs for car $j$ as,\vspace{-0.2cm}
\begin{equation}\label{mott}
\setlength{\abovedisplayskip}{3 pt}
\setlength{\abovedisplayskip}{3 pt}
f(j)=\sum_{i=1}^{W} w_{i}f_{ij},
\end{equation}
where $w_{i}$ captures the weight value of a motif $i$, expressed as $w_{i}=\frac{Z_{i}}{\sum_{j=0}^{W}Z_{j}}$.
After calculating the sum frequency of each car, we sort nodes from the car with the highest frequency to the one with the lowest frequency. Finally, the choice of best $C$ candidates to cache will be the first $C$ elements in the array.

After choosing the serving cars, the next step is to assign each receiving car to its optimal serving car. In particular, if the BS receives the request from one receiving car, the BS would inform the car the nearest serving car with required content based on the collected location information. To solve (\ref{11})-(\ref{1111}), we decompose it into two problems. One is determining the optimal set of serving cars, which is solved by exploiting the temporal motifs. The other is completing the best assignments between the serving cars and the non-serving cars using the spatial knowledge.\vspace{-0.1cm}

\section{Simulation results}\vspace{-0.1cm}

\begin{table}[!t]
    \large
	\begin{center}

		\caption{\small Simulation parameters.}
		\vspace{-0.2cm}
		\label{table_example}
		\resizebox{9cm}{!}{
			\begin{tabular}{|c|c|c|}
				\hline
				\textbf{Parameter} & \textbf{Meaning} & \textbf{Value} \\ \hline \hline
				$ w_{l} $  & Width of each lane & $3.5$~m  \\ \hline
                $P_{b,n}$ &   Transmission power of base station     &  $20$~W     \\ \hline
				$P_{\text{max}}$ & Transmission power of V2V links  & $20$~dBm  \\ \hline
                $\alpha$  & Path loss exponent & $3$ \\ \hline
                $\bar{\gamma}$  & SINR threshold  & $10$~dB  \\ \hline
                $\sigma^{2}$ & Power of noise  & $-94$~dBm  \\ \hline
                $\omega$  &Bandwidth of the system  & $75$~MHz  \\ \hline
                $\theta_{r}$ &Zipf exponent & $2$ \\ \hline
                $M$  & Total number of files in the network   & $10$ \\ \hline
                $F$  & The maximum number of files the car can cache & $3$ \\ \hline
                $d_{b}$ &Approximate distance from BS to freeway & $10.0$~km \\ \hline
                $T$ & Time constraint & 100~s \\ \hline
			\end{tabular}}
		\end{center}\vspace{-0.7cm}
	\end{table}

For our simulations, we use the vehicular mobility dataset within the city of Cologne, Germany, which is collected by the TAPASCologne project \cite{TAPA}.
The dataset contains information about roads and vehicles as well as the trip information for each individual car in one day.
In particular, we take into account a $5$ km-length freeway (Autobahn 4 in
Cologne) with three lanes in each direction and the nearest BS is colonius fernsehturm. We collect the location information of a network of $53$ vehicles that coexist on the freeway from the data.
All simulation parameters are summarized in Table. \ref{table_example}.

Since the Poisson distribution can be used to capture the number of events occurring within a fixed interval \cite{Poisson}, we assume that the number of wireless communication links between two arbitrary cars $i$ and $j$ that in proximity of one another, in a given period of time, follows a Poisson distribution with parameter $\lambda_{i,j}$. To better simulate real-time data, we assume that $\lambda_{i,j}$ is inversely proportional to the distance between $i$ and $j$. This is because a closer distance leads to a better communication environment, resulting in a higher probability to build  communication links. Then, we randomly assign a time stamp to each V2V communication links.
In a baseline randomized V2V communication network, the communication instances are randomly generated and given a time stamp.
Once the temporal graph data is generated, the motifs can be detected based on the approach in Section \uppercase\expandafter{\romannumeral3}. In our analysis, a structure is identified as a motif, when the Z-core of the structure is at least $2$.
For comparison, we use a location-based caching strategy. In this strategy, the basic principle is that the BS would always select a set of cars that can realize the least summation of distance between remaining cars and the corresponding closest cars the in the selected set.

\begin{figure}[!t]
\centering
\includegraphics[width=2.0in,height=0.25in]{Motif_large.pdf}
\DeclareGraphicsExtensions.\vspace{-0.2cm}
\caption{Motifs found in the first scenario.}
\label{Average11}\vspace{-0.5cm}
\end{figure}
\begin{figure}[!t]
\centering
\subfigure[$31$ cars]{
\label{Motif1}
\includegraphics[width=1.33in,height=0.25in]{Motif_31.pdf}}\vspace{-0.1cm}
\subfigure[$33$, $35$, and $37$ cars]{
\label{Motif2}
\includegraphics[width=1.65in,height=0.25in]{Motif_35.pdf}}\vspace{-0.0cm}
\subfigure[$39$ and $41$ cars]{
\label{Motif3}
\includegraphics[width=2.0in,height=0.25in]{Motif_39.pdf}}\vspace{-0.1cm}
\subfigure[$43$, $45$, $47$, $49$, $51$ and $53$ cars]{
\label{Motif4}
\includegraphics[width=1.65in,height=0.25in]{Motif_45.pdf}}\vspace{-0.1cm}
\caption{Motifs found in the second scenario with different number of cars.}
\label{fig:1}\vspace{-0.7cm}
\end{figure}

We consider two simulation scenarios.
In the first scenario, we consider total number of $53$ cars (serving and non-serving). As we change the number of serving cars, the number of non-serving cars will be modified accordingly.
In the second scenario, we first randomly choose several car sets with different total number of cars. Then, based on the proposed strategy and the location-based strategy, we choose a fixed number of non-serving cars chosen out of the selected sets. In particular, we select twelve car sets having a total number of cars ranging from $31$ to $53$ with a step of $2$, and we fix the number of receiving cars to $30$.
According to the proposed method, we can detect different motifs in both scenarios, sorted from the highest to the lowest in terms of Z-score, from the real trace data and generated wireless communication data, as shown in Fig. \ref{Average11} and Fig. \ref{fig:1}.
Based on these structures, we can observe the outdegree of each node. Then, using the proposed algorithm in
Section \uppercase\expandafter{\romannumeral3}, we can determine the set of cars used for caching.

Fig. \ref{Average} shows the average transmission rate achieved by non-serving cars under the location-based caching strategy and the proposed approach, with the total number of cars fixed at $53$ (first scenario). From Fig. \ref{Average},
we observe that the spatio-temporal caching strategy yields a better performance compared with the location-based cache strategy in terms of average date rate per non-serving car. In particular, the performance advantage reaches up to $45\%$ when the number of serving car is $20$. Furthermore, as the number of serving cars increases, the number of non-serving cars will decrease. In particular, when the number of serving cars goes to $50$, there are only $3$ cars requesting for content, leading to a reduced interference and an increase in the average data rate, as also seen in Fig. \ref{Average}.

\begin{figure}[!t]
\centering
\includegraphics[width=2.89in]{large2.pdf}
\DeclareGraphicsExtensions.\vspace{-0.2cm}
\caption{Average data rate for non-serving cars where the number of cars is $53$.}
\label{Average}\vspace{-0.38cm}
\end{figure}
\begin{figure}[!t]
\centering
\includegraphics[width=2.89in]{main3.pdf}
\DeclareGraphicsExtensions.\vspace{-0.2cm}
\caption{Average data rate for non-serving cars where the number of non-serving cars is $30$.}
\label{Pro}\vspace{-0.38cm}
\end{figure}

\begin{figure}[!t]
\centering
\includegraphics[width=2.89in]{CDFICNC.pdf}
\DeclareGraphicsExtensions.\vspace{-0.2cm}
\caption{CDFs for the data rate of non-serving cars where the number of serving cars are $5$, $11$, $15$, $21$.}
\label{Pro11}\vspace{-0.6cm}
\end{figure}

Fig. \ref{Pro} shows the average data rate for $30$ non-serving cars as we vary the number of serving cars. We can observe that the proposed strategy outperforms location-based caching strategy
up to $15\%$ when there are $11$ cars acting as serving nodes.
Further, the average data rate for both strategies will not increase all the time. This is due to the fact that an increase in the number of serving cars will raise the interference over the V2V links. When the impact of interference cannot be compensated by the gain from V2V communication, the average data rate for V2V links will decrease, as seen in Fig. \ref{Pro}. Moreover, Fig. \ref{Pro11} shows the cumulative distribution functions (CDFs) of the data rate of non-serving cars for $5$, $11$, $15$, and $21$ serving cars. Compared with the data rate resulting from location-based strategy, the non-serving cars are more likely to achieve a higher data rate when employing the proposed caching strategy. Fig. \ref{Pro11} also shows under both caching strategies, the probability of achieving a higher data rate for $15$ serving cars is greater than the counterparts for $5$, $11$, and $21$ serving cars. In particular, when the number of serving cars is $15$ and the probability is  $0.8$, the proposed caching strategy improves the data rate of about $20\%$ compared to the location-based strategy.\vspace{-0.3cm}

\section{Conclusions}
In this paper, we have proposed a novel spatio-temporal caching policy in vehicular networks. In contrast to traditional location-based caching strategies, we have leveraged temporal graph motifs, which represent subgraphs with high frequency of occurrence in the V2V communication graph, to determine car candidates to cache popular content. We have developed an approach to detect the motifs and, then, have used the results to determine the preferred set of cars for popular content placement. Simulation results using real car location traces have shown that the proposed spatio-temporal caching strategy can yield significant gains in terms of the average data rate per car for two practical scenarios.
\vspace{-0.5cm}
\def.88{.88}


\section{Introduction} The  relation between the average order of the Goldbach generating function and the Riemann Hypothesis was introduced by Granville \cite{Gra07}.
In an attempt to generalise such relations to sums of primes in arithmetic progressions it was observed that  asymptotic results on the average number of restricted
representations of integers   assuming a conjecture on distinct zeros of Dirichlet L-functions
did not obviously yield an equivalence with the generalised Riemann Hypothesis since the case of the real part of non-trivial zeros being exactly 1 could not
be eliminated. In this note we succeed in doing so for the Riemann zeta function in a simple way and establish an equivalence of average orders of the original Goldbach function with what we call a quasi-Riemann hypothesis, that the real parts of the non-trivial zeros of of the Riemann zeta function
are    strictly less than $1$. This completes the proof of Theorem 1A in \cite{Gra07} i.e. a good asymptotic order
of the Goldbach function is indeed equivalent to the Riemann Hypothesis. The same idea has since been applied  to the restricted sums of two primes both in the same congruence class (Section 8 \cite{BHMS}).

The average order of the function \[ G(n) = \sum_{k_1+k_2=n} \Lambda(k_1) \Lambda(k_2) , \] where $\Lambda$ is the von Mangoldt function
was studied by Fujii \cite{Fuj91} and Granville (op.cit) who showed that under the Riemann Hypothesis  the summatory function can be expressed as
\begin{equation}
 S(x) = \sum_{n\leq x} G(n)
        = \frac{x^2}{2} + O(x^{3/2}) . \end{equation}
The error term was  eventually improved to its conjectured value $O(x^{1+\varepsilon})$ \cite{B_SP} but here it is enough for us to assume any error term less than $x^2$ and from the square root of the associated power series  estimate well the  Chebyshev function,
$\psi(x)=\sum_{n\leq x}\Lambda(n)$,   using a  Dirichlet kernel.

\section {Quasi-Riemann Hypothesis}
  \begin{Th}
  If the asymptotic relation
  \begin{equation}\label{ass}
        S(x) = \frac{x^2}{2} + O(x^{2-\delta}), \ \delta>0
        \end{equation}
holds, then  there exists $0< \delta'  <1  $ such that
 for  the non-trivial zeros $\rho$ of the Riemann zeta function $ \Re \rho < 1- \delta'$.
\end{Th}

\begin{proof}
We consider the power series
 \[ F(z) = \sum _{n=1}^\infty \Lambda(n) z^n, \ |z| <1  \]
whose square is written in terms of the Goldbach function as
  \[ F(z)^2 = \sum_{n=1}^\infty G(n)z^n .\]
 Writing the above as a telescoping series with $G(n)=S(n)-S(n-1)$ we get $$F(z)^2= (1-z) \sum_{n=1}^\infty S(n)z^n.$$
  From the assumption \eq{ass} on the Goldbach asymptotic order we infer

   \begin{equation}\label{asym_sum}
   \sum_{n=1}^\infty  S(n)z^n = \sum_n \left(\frac{n^2}{2} + O(n^{2-\delta}) \right)z^n
   =\sum_n\frac{n^2}{2}z^n+O\left(\sum_n n^{2-\delta}|z|^n\right).
   \end{equation}

   We have
    \[ \sum_n \frac{n^2}{2} z^n =  \frac{1}{(1-z)^3} - \frac{3}{2(1-z)^2} +\frac{1}{2(1-z)} ,\]
 while the error term of (\ref{asym_sum} ) is estimated as
$$\sum_{n=1}^\infty n^{2-\delta}|z|^n \ll (1-|z|)^{\delta-3}  $$
by comparison with the power series expansion of the right hand side
and thus
 $$\sum_n S(n)z^n= \frac{1}{(1-z)^3} +O\left(\frac{1}{(1-|z|)^{3-\delta}  }\right) $$
which implies
 $$F(z)^2= (1-z) \sum_{n=1}^\infty S(n)z^n =  \frac{1}{(1-z)^2}
    + O\zj{ \frac{|1-z|}{(1- |z| )^{3-\delta} }}.$$

To estimate $\psi(N)$ we consider the circle $|z| =R= 1-1/N$ for a large positive integer $N$ and rewrite the last formula as
\begin{equation}\label{square}
F(z)^2 = \frac{1}{(1-z)^2} + O\bigl( |1-z| N^{3-\delta} \bigr).
\end{equation}

The error term is less than the absolute value of the  main term for
\begin{equation}\label{major}  |1-z| <  N^ {\frac{\delta}{3}-1}, \end{equation}
and in this range we take the complex square root of \eq{square} to obtain
\begin{equation}\label{root}
F(z) = \frac{1}{1-z} + O\bigl( |1-z|^2 N^{3-\delta} \bigr).
\end{equation}

\underline {Kernel.} We introduce the function
$$K(z)=z^{-N-1}\frac{1-z^N}{1-z}.$$
We have $K(z)\ll |1-z|^{-1}$ and  we  use it to estimate the partial sums
 \begin{equation}\label{kernel}
 \psi(N) = \frac{1}{2\pi i} \int_{|z|=R}  F(z)K(z) dz \\
= N + \frac{1}{2\pi i} \int_{|z|=R}\left(F(z)-\frac{1}{1-z}\right)K(z) dz .
  \end{equation}

\underline{Major Arc.} On the part of the circle \eq{major} the integrand is
\[  O\bigl( |1-z| N^{3-\delta} \bigr) = O\bigl(  N^{2-2\delta/3} \bigr)\]
so the contribution to the integral above is $ O\bigl(  N^{1-\delta/3} \bigr)$
from the power series estimate obtained in (\ref{root}).

\smallskip
\underline {Minor Arc.} On the rest of the circle, i.e. the arc $|1-z|>N^{\delta/3-1}$ we use the Cauchy-Schwarz inequality.
For $F$ we apply
the  estimate for the full circle

\begin{equation*}
\begin{split}
\int_{|z|=R} \ab{F(z)-\frac{1}{1-z}}^2 |dz|
&=\int_0^{2\pi}\sum_n(\Lambda(n)-1)R^ne^{int}\sum_m(\Lambda(n)-1)R^ne^{-imt}\\
&=2\pi \sum_{n=1}^\infty (\Lambda(n)-1)^2 \zj{1-\frac{1}{N}}^{2n} = O(N \log N).
\end{split}
\end{equation*}
We estimate the {\it square integral of the kernel} as follows.
\[
I= \int_{\substack{|z|=R \\ |1-z|>N^{\delta/3-1}}} |K(z)|^{2} |dz| \ll
   \int_{\substack{|z|=R \\ |1-z|>  N^{\delta/3-1} } }  \frac{ |dz| }{|1-z|^{2}}.
\]
With the parametrization $ z=Re^{it}$ the last integral becomes
 \[ 2\int_{t_0}^\pi \frac{1}{ |1-Re^{it} |^2} dt, \]
where $t_0$ is defined by $ |1-Re^{it_0} | =  N^{\delta/3-1} $.

We have
 \[  |1-Re^{it} |^2 = 1 + R^2 - 2R \cos t = (1-R)^2 + 4R^2 \sin^2 \frac{t}{2} \
\begin{cases} < t^2 + N^{-2} \\ > t^2/3, \end{cases}  \]
whence $t_0 \gg  N^{\delta/3-1} $ and the integral satisfies
 \[ I \ll \int_{t_0}^\pi t^{-2} \, dt < t_0^{-1} \ll N^{1-\delta/3} . \]

Now the Cauchy Schwarz inequality gives the bound for the
 minor arc as
\[
  \Big( \int_{|z|=R} \Big|F(z)-\frac{1}{(1-z)}\Big|^{2}dz
  \Big)^{1/2} I^{1/2}
  \ll  N^{1-\delta/6} ( \log N)^{1/2} .
\]
 The major and minor arcs together give
\[
  \psi(N)-N \ll  N^{1-\delta/6} ( \log N)^{1/2} ,
\]
which proves the theorem with $\delta'=\delta/6$.

   \end{proof}
\section{Remarks}
 1. Using a better kernel and a more careful calculation we can improve the bound to $\delta'=\delta/3$. It is of little importance in this context
as the methods of [1,2,4] give $\delta'=\delta$ as soon as we have $\delta'<1$, through an explicit expression of $S(x)$ using roots of the zeta-function.
 We are unable to get $\delta'=\delta$ directly by the simple method of this paper and we think it cannot be done. We used very little of the properties
 of the Liouville function and our method works equally for any function $f$ satisfying $ |f(n)| = O(n^\varepsilon)$ for all $\varepsilon>0$ and
  \[ \sum_{m+n\leq x} f(m)f(n) = cx^2 + O(x^{2-\delta}) .   \]
We doubt that
\[ \sum_{n\leq x} f(n) = c'x + O(x^{1-\delta+o(1)})  \]
follows in this generality.

2. In general we are yet unable to establish an equivalence between average orders of
restricted Goldbach functions and  appropriate quasi Riemann Hypotheses for $L$-functions of \cite{Gra07} or \cite{BHMS} though a particular case is handled
in [1] using the method of this note.


\section{Introduction}

Studies have shown that approximately 80\% of Americans would like to spend their final days at home if possible, but only 20\% do \cite{AMERICANSDIE}. In fact, up to 60\% of deaths happen in an acute care hospital, with patients receiving aggressive care in their final days. Access to palliative care services in the United States has been on the rise over the past decade. In 2008, 53\% of all hospitals with fifty or more beds reported having palliative care teams, rising to 67\% in 2015\cite{DumanovskySpecialReport}. However, despite increasing access, data from the National Palliative Care Registry estimates that less than half of the 7-8\% of all hospital admissions that need palliative care actually receive it\cite{REPORTSITE}. Though a significant reason for this gap comes from the palliative care workforce shortage \cite{Spetz2016FewRecommendations.}, and incentives for health systems to employ them, technology can still play a crucial role by efficiently identifying patients who may benefit most from palliative care, but might otherwise be overlooked under current care models.

We focus on two aspects of this problem.  First, physicians may not refer patients likely to benefit from palliative care for multiple reasons such as overoptimism, time pressures, or treatment inertia \cite{Christakis2000ExtentStudy}.  This may lead to patients failing to have their wishes carried out at end of life\cite{STANFMORT} and overuse of aggressive care.  Second, a shortage of palliative care professionals makes proactive identification of candidate patients via manual chart review an expensive and time-consuming process.

The criteria for deciding which patients benefit from palliative care can be hard to state explicitly.
Our approach uses deep learning \fixed{Here and for at least first half of paper machine learning is used instead of more specific deep learning in contrast to title} to screen patients admitted to the hospital to identify those who are most likely to have palliative care needs.
The algorithm addresses a proxy problem - to predict the mortality of a given patient within the next 12 months - and use that prediction for making recommendations for palliative care referral.
This frees the palliative care team from manual chart review of every admission and helps counter the potential biases of treating physicians by providing an objective recommendation based on the patient's EHR. Currently existing tools to identify such patients have limitations, and they are discussed in the next section.

\section{Related work}

Accurate prognostic information is valuable to patients, caregivers, and clinicians \cite{KutnerInformationIllness} \cite{Steinhauser2001PreparingProviders}. Several studies have shown that clinicians are generally over optimistic in their estimates of the prognoses of terminally ill patients\cite{Selby2011ClinicianCategories} \cite{Christakis2000ExtentStudy} \cite{ViganoTheCancer} \cite{Glare2008PredictingDisease}. It has also been shown that no subset of clinicians are better at late stage prognostication than than others  \cite{White2016AExperts} \cite{Macmillan1998AccuracyHospice} . However, clinician judgment remains the most common method of predicting survival in practice \cite{White2016AExperts}. Several solutions exist that attempt to make patient prognosis more objective and automated. Many of these solutions are models that produce a score based on the patient's clinical and biological parameters, and can be mapped to an expected survival rate.

\subsection*{Prognostic tools in Palliative Care}
The Palliative Performance Scale  \cite{Lau2006UsePrognostication} was developed as a modification of the Karnofsky Performance Status Scale (KPS) \cite{Karnofsky1948TheCarcinoma} to the Palliative care setting, and is calculated based on observable factors such as: degree of ambulation, ability to do activities, ability to do self-care, food and fluid intake, and state of consciousness. The Palliative Prognostic Score (PPS) was constructed for the Palliative Care setting as well, focusing on terminally ill cancer patients \cite{Pirovano1999ACare}.  The PPS is calculated with multiple regression analysis based on the following variables: Clinical Prediction of Survival (CPS), Karnofsky Performance Status (KPS), anorexia, dyspnea, total white blood count (WBC) and lymphocyte percentage. The Palliative Prognostic Index (PPI),  developed around the same time as PPS, also calculates a multiple regression analysis based score using Performance Status, oral intake, edema, dyspnea at rest, and delirium. These scores are difficult to implement at scale since they involve face-to-face clinical assessment and involve prediction of survival by the clinician. Furthermore, these scores were designed to be used within the palliative care setting, where the patient is already in an advanced stage of the disease --- as opposed to identifying them earlier.
\subsection*{Prognostic tools in the Intensive Care Unit}
There also are prognosis scoring models that are commonly used in the Intensive Care Unit. The APACHE-II (Acute Physiology, Age, Chronic Health Evaluation) Score predicts hospital mortality risk for critically ill hospitalized adults in the ICU \cite{Knaus1985APACHESystem.}. This model has been more recently refined with the APACHE-III Score, which uses factors such as major medical and surgical disease categories, acute physiologic abnormalities, age, preexisting functional limitations, major comorbidities, and treatment location immediately prior to ICU admission\cite{Knaus1991TheSystem}. Another commonly used scoring system in the ICU is the Simplified Acute Physiological Score, or SAPS II \cite{LeGall1993AStudy}, which is calculated based on the patient's physiological and underlying disease variables. While these score are useful for the treatment team when the patient is already in the ICU, they have limited use in terms of identifying patients who are at risk of longer term mortality, while they are still capable of having a meaningful discussion of their goals and values, so that they can be set on an alternative path of care.

\subsection*{Prognostic tools for Early Identification}
There have been a number of studies and tools developed that aim to identify terminally ill patients early enough for an end-of-life plan and care to be meaningful.

CriSTAL (Criteria for Screening and Triaging to Appropriate aLternative care) was developed to identify elderly patients nearing end of life, and quantifies the risk of death in the hospital or soon after discharge  \cite{Cardona-Morrell2015DevelopmentCriSTAL.}. CriSTAL provides a check list using eighteen predictors with the goal of identifying \textit{the dying patient}.

CARING is a tool that was developed to identify patients who could benefit from palliative care\cite{FischerACriteria}. The goal was to use six simple criteria in order to identify patients who were at risk of death within 1 year. PREDICT \cite{Richardson2015PREDICT:Directive} is a screening tool also based on six prognostic indicators, which were refined from CARING. The model was derived from 976 patients.

The Intermountain Mortality Risk score is an all-causes mortality prediction based on common laboratory tests \cite{Horne2009ExceptionalTests}. The model provides score for 30-day, 1-year and 5-year mortality risk. It was trained on a population of 71,921 and tested on 47,458.

\fixed{The paragraph structure here seems too explicitly formulaic}

Cowen, M et al \cite{Cowen2013MortalityActivities} proposed using a twenty-four factor based prediction rule at the time of hospital admission to identify patients with high risk of 30-day mortality, and to organize care activities using this prediction as a context. One of the their motivation was to have a rule from a single set of factors, and not be disease specific.  The model was derived from 56,003 patients.

Meffert, C et al \cite{Meffert2016IdentificationScore} proposed a scoring method based on logistic regression on six factors to identify hospitalized patients in need of palliative care. In this prospective study, they asked the treating physician at the time of discharge whether the patient had palliative care needs. The trained model was then used to identify such patients at the time of admission. The model was derived from 39,849 patients.

Ramachandran, K et al \cite{Ramchandran2013ARecord} developed a 30-day mortality prediction tool for hospitalized cancer patients. Their model used eight variables that were based on information from the first 24 hours of admission, and laboratory results and vitals. A logistic regression model was developed from these eight variables and used as a scoring function. The model was derived from 3,062 patients.

Amarasingham, R et al \cite{AmarasinghamAnData} built a tool to screen patients who were admitted with heart failure, and identify those who are at risk of 30-day readmission or death. Their regression model uses a combination of Tabak Morality Score \cite{Tabak2007UsingAdjustment}, markers of social, behavioral, and utilization activity that could be obtained electronically,  ICD-9 CM codes  specific to depression and anxiety, billing and administrative data. Though this study was not specifically focused on palliative care, the methodology of using EHR system data is relevant to our work. The model was derived from 1,372 patients.

Makar, M et al \cite{MakarShort-TermData} used only Medicare claims data on older population ($\ge$ 65 years) to predict mortality in six months. By limiting their model to use only administrative data, they hypothesized an easier deployment scenario thereby making automated prognostic models more prevalent. The model was derived separately on four cohorts (one per disease type) with 20,000 patients per cohort.

\subsection*{Prognosis in the age of Big-Data}
The proliferation of EHR systems in healthcare combined with advances in Machine Learning techniques on high dimensional data provides a unique opportunity to make contributions, especially in disease prognosis \cite{Sim2016TwoMedicine}\cite{Obermeyer2016PredictingMedicine.}. All the tools described above, and those we reviewed \cite{Rose2013MortalityLearning}  \cite{Wiemken2017PredictingApproaches.} \cite{MotwaniMachineAnalysis}  \cite{CooperAnMortality} \cite{VomlelMachine}, have at least one of the following limitations. They were either derived from small data sets (limited to specific studies or cohorts), or used too few variables (intentionally to make the model portable, or avoid overfitting), or the model was too simple to capture the complexities and subtleties of human health, or was limited  to certain sub-populations (based on disease type, age etc.) We address these limitations in our work.

\section{Methods}

We approach the problem of predicting mortality from the point of view of the palliative care team by being largely agnostic to disease type, disease stage, severity of admission (ICU vs non-ICU), age etc. We take a data driven approach and build a deep learning model that considers every patient in the EHR (with a sufficiently long history), without limiting our analysis to any specific sub-population or cohort. In order to make the problem of identifying patients with palliative care needs tractable, we use the following proxy problem statement instead:

\textit{Given a patient and a date, predict the mortality of that patient within 12 months from that date, using EHR data of that patient from the prior year.}

We treat this as a binary classification problem and build a supervised deep learning model to solve it.
Other than building a model that performs well on the above problem, we are also separately interested in the model performance on a sub-problem --- the ability to predict mortality of patients who are currently admitted. This is because it is much easier for the palliative care staff to intervene with admitted patients.

\subsection*{Data Source}

STRIDE (Stanford Translational Research Integrated Database Environment) \cite{LoweSTRIDEPlatform} is a clinical data warehouse supporting clinical and translational research at Stanford University. The snapshot of STRIDE used in our work comprises the EHR data of approximately 2 million adult and pediatric patients cared for at either the Stanford Hospital or the Lucile Packard Children's hospital between 1995 and 2014.

\subsection*{Constructing a Dataset for Supervised Learning}
Patients who have a recorded date of death are considered \textit{positive cases}; other patients are considered \textit{negative cases}. Further, we define the \textit{prediction date} of a patient to be the point in time that divides their health record timeline into virtual future and past events.  We use data from each patient's virtual past to make predictions about their survival 3-12 months in the future.  Note that we must take care when defining the \textit{prediction date} to not violate common sense constraints (described below) that could invalidate the labels.  We only include patients for whom it is possible to find a \textit{prediction date} that satisfies these constraints.

\begin{figure}
    \includegraphics[width=\columnwidth]{km}
    \caption{Right-censoring lengths shown as a survival plot. \zxie{Not sure whether Kaplan Meier plot would be good to include as well.}}
\label{fig:censor}
\end{figure}

\subsubsection*{Positive Cases}
The constraints for positive cases were decided based on the rationale that palliative care is most beneficial if the referral occurs 3-12 months prior to death. Predicting mortality within 3 months is considered too late due to the preparatory time required to start palliative care in general. On the other hand, a lead time longer than 12 months is problematic because making accurate predictions over such a long time horizon is difficult, and more importantly, palliative care interventions are a limited resource that are best focused on more immediate needs. The \textit{prediction date} for positive cases must meet all the following constraints:

\begin{itemize}
\item The \textit{prediction date} must be a recorded date of encounter.
\item The \textit{prediction date} must be at least 3 months prior to date of death (otherwise death is too near).
\item The \textit{prediction date} can be at most 12 months prior to date of death (otherwise death is too far).
\item The \textit{prediction date} must be at least 12 months after the date of first encounter (otherwise the patient lacks sufficient history on which to base a prediction).
\item In-patient admissions are preferred over other admission types for the \textit{prediction date}, as long as they meet the previous constraints (since it is easier to start the palliative care conversation with them).
\item The \textit{prediction date} must be the earliest among the possible candidate dates subject to previous constraints.
\end{itemize}

\subsubsection*{Negative Cases}
For negative cases (patients without a date of death), we require that the patient was alive for at least 12 months from the \textit{prediction date}.  We choose the \textit{prediction date} such that it satisfies all the following constraints:

\begin{itemize}
\item The \textit{prediction date} must be a recorded date of encounter.
\item The \textit{prediction date} must be at least 12 months prior to date of last encounter (to avoid ambiguity of death after date of EHR snapshot).
\item The \textit{prediction date} must be at least 12 months after the date of first encounter (otherwise insufficient history).
\item In-patient admissions are preferred over other encounter types for the  \textit{prediction date}, as long as they meet the previous constraints (to serve as controls for the admitted positive cases).
\item The \textit{prediction date} must be the latest among the possible candidate dates subject to previous constraints.
\end{itemize}

\subsubsection*{Admitted patients}
Those patients whose \textit{prediction date} corresponds to an in-patient admission are considered \textit{admitted patients}. Remaining patients are considered non-admitted (note that non-admitted patients could still have other recorded admissions in their history).  Further, for \textit{admitted patients}, their \textit{prediction date} it is re-adjusted by incrementing it to be the second day of admission. The rationale for doing this is that patient records are generally updated with the latest data (preliminary tests, diagnostics etc.) within 24 hours of admission, and the second day is better suited for making a more informed prediction. Note that the \textit{admitted patients} are a subset of the included patients (and NOT a separate set).

All the data after the \textit{prediction date} is censored in both the positive and negative cases. The KM-plot of censor lengths is shown in Fig \ref{fig:censor},  highlighting the separation between the two classes at 365 days.

\subsection*{Data Description}

\begin{table}
\begin{center}
\begin{tabular}{|l|r|r|r|}
\hline
       & Alive & Deceased  & Total\\ \hline
In EHR & 1,880,096 & 131,009 & 201,1105 \\ \hline
Selected & 205,571 & 15,713& 221,284 \\ \hline
Admitted & 9,648 & 1,131 & 10,779 \\ \hline
\end{tabular}
\newline
\newline
\caption{Breakdown of patient counts.}
\label{table:inclusion}
\end{center}
\end{table}

The inclusion criteria selected a total of 221,284 patients. Table ~\ref{table:inclusion} shows the breakdown of these patients based on inclusion and admission.

\begin{figure}
    \includegraphics[width=\columnwidth]{age}
    \caption{Age of patients at \textit{prediction time}. \fixed{If you have additional dimensions that could be plotted I would include some more.}}
    \label{fig:age}
\end{figure}

We observe that, unsurprisingly, the distribution of age at prediction time is not equal between the classes, and that the positive class (of deceased patients) is skewed towards older age (Fig ~\ref{fig:age}).

\begin{table}
\begin{center}
\begin{tabular}{|l|r|r|r|r|}
\hline
       & Training & Validation  & Testing & \\ \hline
Alive & 164,424 & 20,619 & 20,528 & 205,571\\ \hline
Deceased & 12,587 & 1,520 &  1,606 & 15,713\\ \hline
 & 177,011& 22,139 & 22,134 & 221,284\\
\hline
\end{tabular}
\newline
\newline
\caption{Data split for modeling.}
\label{tab:splits}
\end{center}
\end{table}

The included patients are randomly split in approximate ratio 8:1:1 into training, validation and test sets, as shown in Table \ref{tab:splits}.

The prevalence of death among the included patients is approximately 7\%.  Approximately 5\% were \textit{admitted patients} (i.e., \textit{prediction date} was the second day of an admission). Among the \textit{admitted patients}, the prevalence of death is about 11\%.

\subsection*{Feature Extraction}

For each patient, we consider the 12 months leading up to their \textit{prediction date} as their \textit{observation window}. Within the \textit{observation window} of each patient, we use ICD9 (International Classification of Diseases 9th rev) diagnostic and billing codes, CPT (Current Procedural Terminology) procedure codes, RxNorm prescription codes, and encounters found in that period to create features.

\begin{table}
\begin{center}
\begin{tabular}{|l|l|l|r|}
\hline
 & Start date & End date & Duration\\
\hline
\textit{observation window} & PD - 365 & PD & 365 \\
\hline
\textit{observation slice 1} & PD - 30 & PD & 30\\
\hline
\textit{observation slice 2} & PD - 90 & PD - 30 & 60 \\
\hline
\textit{observation slice 3} & PD - 180 & PD - 90 & 90 \\
\hline
\textit{observation slice 4} & PD - 365 & PD - 180 & 185 \\
\hline
\end{tabular}
\newline
\newline
\caption{Observation window and slices}
\label{tab:obswindow}
\end{center}
\end{table}

We create features as follows.  In order to capture the longitudinal nature of the data, we split the \textit{observation window} of each patient into four \textit{observation slices}, specified relative to the \textit{prediction date} (PD) as shown in Table \ref{tab:obswindow}

Thus, \textit{observation slice 1} is the most recent, and 4 is the oldest. The slice widths are intentionally uneven in order to give more emphasis to recent data. Within each \textit{observation slice}, we count the the number of occurrences of each code in each code category (prescription, billing, etc.) per patient. The count of every such code within the slice is considered a separate feature.

We also include the patient demographics (age, gender, race and ethnicity), and the following per-patient summary statistics in the \textit{observation window} for each code category:
\begin{itemize}
\item Count of unique codes in the category.
\item Count of total number of codes in the category.
\item Maximum number of codes assigned in any day.
\item Minimum number of codes (non-zero) assigned in any day.
\item Range of number of codes assigned in a day.
\item Mean of number of codes assigned in a day.
\item Variance in number of codes assigned in a day.
\end{itemize}

All these features (i.e, code counts in each of the four \textit{observation slices}, per category summary statistics over the \textit{observation window}, and demographics) were concatenated to form the candidate feature set. From this set, we pruned away those features which occur in 100 or fewer patients. This resulted in the final set of 13,654 features. Of the 13,654 features, each patient on average has 74 non-zero values (with a standard deviation of 62), and up to a maximum of 892 values. The overall feature matrix is approximately 99.5\% sparse.

\subsection*{Algorithm and Training}

Our model is a Deep Neural Network (DNN) \cite{LeCun2015DeepLearning} comprising an input layer (of 13,654 dimensions), 18 hidden layers (each 512 dimensions) and a scalar output layer. We employ the logistic loss function at the output layer and use the Scaled Exponential Linear Unit (SeLU) activation function \cite{Klambauer2017Self-NormalizingNetworks} at each layer. The model is optimized using the Adam optimizer \cite{Kingma2014Adam:Optimization}, with a mini-batch size of 128 examples. Intermediate model snapshots were taken every 250 mini-batch iterations, and the snapshot that performed best on the validation test was selected as the final model. Explicit regularization was not found necessary. The network configuration was reached by extensive hyperparameter search over various network depths (ranging from 2 to 32) and activation functions ($tanh$, $ReLU$ and $SeLU$).

\fixed{If you have GBM results may be nice to include as a baseline.}

The software was implemented using the Python programming language (version 2.7), PyTorch framework \cite{PyTorch}, and the scikit-learn library (version 0.17.1) \cite{scikit-learn}. The training was performed on an NVIDIA TitanX (12GB RAM) with CUDA version 8.0.

\subsection*{Evaluation}

Since the data is imbalanced (with 7\% prevalence), accuracy can be a poor evaluation metric \cite{HaiboHe2009LearningData}. The ROC curve can also be sometimes misleading on imbalanced problems\cite{DavisTheCurves} \cite{BoydUnachievableEvaluationb}. Therefore, we use the Average Precision (AP) score, also known as Area Under Precision-Recall Curve (AUPRC) for model selection \cite{Yuan2015Threshold-freeTests.}.

\section{ Results}

\begin{table*}[t]
\begin{center}

\begin{tabular}{|l|l|}
\hline
Patient MRN & XXXXXXX \\ \hline
Probability score & 0.946 \\ \hline
\end{tabular}
\newline
\begin{tabular}{|l|l|r|r|l|}
\hline
Factors & Code & Value & Influence & Description \\ \hline
Top Diagnostic factors              & V10.51    &         4   &   0.0051  &  Personal history of malignant neoplasm of bladder \\
              & V10.46     &        5    &  0.0019   & Personal history of malignant neoplasm of prostate \\
              & 518.5  &            1   &   0.0012   & Pulmonary insufficiency following trauma and surgery \\
              & 518.82   &          1    &  0.0008   & Other pulmonary insufficiency \\
               & 88.75    &         1    &  0.0006 &   Diagnostic ultrasound of urinary system \\ \hline
 Top Procedural factors             & 88331  &           1   &   0.0017  &  Pathology consultation during surgery with FS \\
               &75984       &      1    &  0.0014   & Transcatheter Diagnostic Radiology Procedure \\
               &72158       &      1   &   0.0013&    MRI and CT Scans of the Spine \\
     &Code\_Type\_Count    &        76   &   0.0011&   Summary statistic (count of all ICD/CPT codes) \\
              & 76005    &         1   &   0.0007&  Fluroscopic guidance and localization of needle or catheter tip for spine  \\      \hline
Top Medication factors& & & & \\ \hline

 Top Encounter factors           &Hx Scan    &        21   &   0.0012&  Number of scan encounters of all types \\
           &Inpatient     &       60     & 0.0004&  Number of days patient was admitted \\
   &Var\_Codes\_per\_Day    &         8    &  0.0002&  Summary statistic (variance in number of codes assigned per day)  \\
      &Code\_Day\_Count    &        88  &    0.0001&   Number of days any encounter code was assigned \\ \hline
   Top Demographic factors              & Age    &        81  &    0.0010& Age of patient in years at \textit{prediction time} \\ \hline
\end{tabular}
\newline
\newline
\caption{Prediction explanation generated on a random positive patient with high  probability score. Only factors that contributed to a drop in probability score are reported. \fixed{This table could use some formatting...tough to parse}}
\label{tab:explainpositive}
\end{center}
\end{table*}

\begin{table*}[t]
\centering
\begin{tabular}{|l|l|}
\hline
Patient MRN & YYYYYYY \\ \hline
Probability score & 0.909 \\ \hline
\end{tabular}
\newline
\begin{tabular}{|l|l|r|r|l|}
\hline
 Factors& Code & Value & Influence & Description \\ \hline

 Top Diagnostic factors              &197.7  &          16   &   0.1299 &  Malignant neoplasm of liver, secondary \\
               &154.1   &          3    &  0.1254&    Malignant neoplasm of rectum  \\
               &287.5    &         1  &    0.0194    &   Thrombocytopenia, unspecified \\
               &780.6    &         1   &   0.0171&   Fever and other physiologic disturbances of temperature regulation \\
              &733.90    &         1   &   0.0113&   Other and unspecified disorders of bone and cartilage \\ \hline
 Top Procedural factors            & 73560  &           1    &  0.0502&    Diagnostic Radiology (Diagnostic Imaging) Procedures of the Lower Extremities \\
    & Code\_Type\_Count &           20   &   0.0491&     Summary statistic (Number of unique ICD-9/CPT codes) \\
          &     74160    &         1    &  0.0381&     Diagnostic Radiology (Diagnostic Imaging) Procedures of the Abdomen \\
  & Max\_Codes\_per\_Day &            6  &    0.0234&    Summary statistic (Maximum number of codes in any day) \\
& Range\_Codes\_per\_Day  &           6  &    0.0233&     Summary statistic (Range of codes across days) \\ \hline
  Top Medication factors       &     283838  &           1   &   0.0619&     Darbepoetin Alfa  \\
          &     28889   &          1    &  0.0247&     Loratadine \\
 &Range\_Codes\_per\_Day   &          5  &    0.0023&     Summary statistic (Ranges of codes across days) \\
   &Max\_Codes\_per\_Day   &          5  &    0.0023&     Summary statistic (Maximum number of codes in any day) \\
   &  Code\_Type\_Count    &         6   &   0.0015&     Summary statistic (Number of unique medication codes) \\ \hline
 Top Encounter factors            &Hx Scan       &     19   &   0.2239&     Number of scan encounters of all types \\
     & Code\_Day\_Count    &        97  &    0.0284&     Number of days any encounter code was assigned \\
       &   Outpatient      &      22    &  0.0228&     Number of Outpatient encounters \\
  & Var\_Codes\_per\_Day    &         1  &    0.0074&     Summary statistic (variance in number of codes assigned per day) \\ \hline
Top Demographic factors&     & & & \\ \hline
\end{tabular}
\newline
\newline
\caption{Prediction explanation generated on a random  false positive patient with high probability score. Only factors that contributed to a drop in probability score are reported.}
\label{tab:explainfp}
\end{table*}
\begin{figure}
    \includegraphics[width=\columnwidth]{calibration}
    \caption{Reliability curve (calibration plot) of the model output probabilities on the test set data.}
\label{fig:calibration}
\end{figure}

\begin{figure}
    \includegraphics[width=\columnwidth]{pr}
    \caption{Interpolated Precision-Recall curve. The horizontal dotted line represents precision level of 0.9. The vertical dotted lines indicate the recall at which the curves achieve 0.9 precision.}
\label{fig:pr}
\end{figure}

\begin{figure}
    \includegraphics[width=\columnwidth]{roc}
    \caption{Receiver Operating Characteristic (ROC) of the model performance on the test set.}
\label{fig:roc}
\end{figure}

In this section we report technical evaluation results obtained on the test set using the model selected based on the best AP score on the validation set.

We observe that the model is reasonably calibrated (Fig \ref{fig:calibration}) with a \textbf{Brier score of 0.042}. In the high threshold regime, which is of interest to us, the model is a little conservative (under-confident) in its probability estimates, which should not hurt.

The interpolated Precision-Recall curve is shown in Fig \ref{fig:pr}. The model achieves an \textbf{AP score of 0.69} (0.65 on \textit{admitted patients}). Early recall is desirable, and therefore Recall at precision 0.9 is a metric of interest. The model achieves \textbf{recall of 0.34 at 0.9 precision} (0.32 on  \textit{admitted patients}). The Receiver Operating Characteristic curve is shown in Fig \ref{fig:roc}.  The model achieves an \textbf{AUROC of 0.93} (0.87 for \textit{admitted patients}). Both the ROC and Precision-Recall plots suggest that the model demonstrates strong early recall behavior. \fixed{Would be nice to have some qualitative statements about requirements of palliative care team and how the model exceeds these requirements and does well--baseline would also help for this purpose}

\subsection*{Qualitative Analysis}
\fixed{Would say more here or make a subsection / roll into another section}
It is worth recalling that predicting mortality was a proxy problem for identifying patients who could benefit from palliative care. In order to evaluate our performance on the original problem, we inspected false positives with high output probability.  Although such patients did not die within 12 months from their prediction dates, we noted that they were often diagnosed with terminal illness and/or are high utilizers of healthcare services. This can be seen in the positive and false positive examples shown in Section \ref{section:explain}.

Upon conducting a chart review of 50 randomly chosen patients in the top 0.9 precision bracket of the test set, the palliative care team found all were appropriate for a referral on their \textit{prediction date}, even if they survived more than a year. This suggests that mortality prediction was a reasonable (and tractable) choice of a proxy problem to solve.

\section{Explaining Predictions}
\label{section:explain}

Supervised machine learning techniques, and in particular Deep Learning techniques, have recently demonstrated tremendous success in predictive ability. However, better performance often requires larger, more complex models and thus sacrifice interpretability. It is worth drawing a distinction between interpreting a model, versus interpreting its decision \cite{RibeiroquotClassifierb} \cite{LiptonTheInterpretability}. While interpreting complex models (e.g very deep neural networks) may sometimes be infeasible, it is often the case that users only want an explanation for the prediction made by the model for a given example. It is important to establish the trust of the practitioner in the model's decisions for them to feel comfortable taking actions based on it. Providing explanations along with decisions help establish that trust.

We make the following observations to motivate our explanation technique.

\begin{itemize}
\item We can view the EHR data as a strictly growing log of events, and that new data is only added (nothing is modified or removed in general). This results in all our features being positive valued (as counts, means and variance of counts, etc).

\item We are most interested in explaining why a model assigns high probability to a patient. We are less interested in getting an explanation for why a healthy person was given a low probability (the reasons are also much less clear: the patient did not have brain cancer, did not have pneumonia, and so on).

\item Directly perturbing feature vectors (e.g sensitivity analysis or for techniques described in \cite{RibeiroquotClassifierb}) does not work well in our case . For example, perturbing the feature representing the ICD count for brain cancer from zero to non zero can increase the probability of death significantly, implying that it is an important factor in general.  However, that is not a very useful observation for a \emph{specific patient who does not have brain cancer}.

\end{itemize}

These observations motivate the following technique. For each ICD-9, CPT, RXNORM and Encounter code, we ablate \textit{all occurrences} of that code from the patient's EHR, create a new feature vector, and measure the drop in probability compared to the original probability.  This corresponds to asking: all else being equal, how would the probability change if this patient was not diagnosed with XYZ, prescribed drug ABC, etc? This drop in probability is considered the influence the code has on the model's decision for that patient. Demographic features are handled as follows. We zero out the age and  swap the gender to the opposite sex, and measure the respective drops in probability. Finally we sort the codes in descending order by influence, and pick the top 5 in each code category.  A random example of such a positive and false positive case are shown in Table \ref{tab:explainpositive} and \ref{tab:explainfp}.

\section{Conclusion}
We demonstrate that routinely collected EHR data can be used to create a system that prioritizes patients for follow up for palliative care \fixed{Seems stronger to say "We demonstrate that...". Also unsure if "model" or "system" would be better for this audience.}. In our preliminary analysis we find that it is possible to create a model for all-\fixed{typo}cause mortality prediction and use that outcome as a proxy for the need of a palliative care consultation\fixed{as a proxy for palliative care need?}. The resulting model is currently being piloted for daily, proactive outreach to newly admitted patients. We will collect objective outcome data (such as rates of palliative care consults, and rates of goals of care documentation) resulting from the use of our model \fixed{from the use of the model}. We also demonstrate a novel method of generating explanations from complex deep learning models that helps build confidence of practitioners to act on the recommendations of the system.

\section*{Acknowledgment}
We thank the Stanford Research IT team for their  support and help in this project. Research IT, and the Stanford Clinical Data Warehouse (CDW) are supported by the National Center for Research Resources and the National Center for Advancing Translational Sciences, National Institutes of Health, through grant UL1 TR001085. The content of studies done using the CDW is solely the responsibility of the authors and does not necessarily represent the official views of the NIH.