I think it depends on how it is to be used.
Just for reference for other readers, if $P$ and $Q$ are probability measures, then the Jensen-Shannon Divergence is
$$
J(P,Q) = \frac{1}{2} \big( D(P \mid\mid R) + D(Q\mid\mid R) \big)
$$
where $R = \frac{1}{2} (P + Q)$ is the mid-point measure and $D(\cdot\mid\mid\cdot)$ is the Kullback-Leibler divergence.
Now, I would be tempted to use the square root of the Jensen-Shannon Divergence since it is a metric, i.e. it satisfies all the "intuitive" properties of a distance measure.
For more details on this, see
Endres and Schindelin, A new metric for probability distributions, IEEE Trans. on Info. Thy., vol. 49, no. 3, Jul. 2003, pp. 1858-1860.
Of course, in some sense, it depends on what you need it for. If all you are using it for is to evaluate some pairwise measure, then any monotonic transformation of JSD would work. If you're looking for something that's closest to a "squared-distance", then the JSD itself is the analogous quantity.
Incidentally, you might also be interested in this previous question and the associated answers and discussions.
The midpoint measure $\newcommand{\bx}{\mathbf{x}} \newcommand{\KL}{\mathrm{KL}}M$ is a mixture distribution of the two multivariate normals, so it does not have the form that you give in the original post. Let $\varphi_p(\bx)$ be the probability density function of a $\mathcal{N}(\mu_p, \Sigma_p)$ random vector and $\varphi_q(\bx)$ be the pdf of $\mathcal{N}(\mu_q, \Sigma_q)$. Then the pdf of the midpoint measure is
$$
\varphi_m(\bx) = \frac{1}{2} \varphi_p(\bx) + \frac{1}{2} \varphi_q(\bx) \> .
$$
The Jensen-Shannon divergence is
$$
\mathrm{JSD} = \frac{1}{2} (\KL(P\,\|M)+ \KL(Q\|M)) = h(M) - \frac{1}{2} (h(P) + h(Q)) \>,
$$
where $h(P)$ denotes the (differential) entropy corresponding to the measure $P$.
Thus, your calculation reduces to calculating differential entropies. For the multivariate normal $\mathcal{N}(\mu, \Sigma)$, the answer is well-known to be
$$
\frac{1}{2} \log_2\big((2\pi e)^n |\Sigma|\big)
$$
and the proof can be found in any number of sources, e.g., Cover and Thomas (1991), pp. 230-231. It is worth pointing out that the entropy of a multivariate normal is invariant with respect to the mean, as the expression above shows. However, this almost assuredly does not carry over to the case of a mixture of normals. (Think about picking one broad normal centered at zero and another concentrated normal where the latter is pushed out far away from the origin.)
For the midpoint measure, things appear to be more complicated. That I know of, there is no closed-form expression for the differential entropy $h(M)$. Searching on Google yields a couple potential hits, but the top ones don't appear to give closed forms in the general case. You may be stuck with approximating this quantity in some way.
Note also that the paper you reference does not restrict the treatment to only discrete distributions. They treat a case general enough that your problem falls within their framework. See the middle of column two on page 1859. Here is where it is also shown that the divergence is bounded. This holds for the case of two general measures and is not restricted to the case of two discrete distributions.
The Jensen-Shannon Divergence has come up a couple of times recently in other questions on this site. See here and here.
Addendum: Note that a mixture of normals is not the same as a linear combination of normals. The simplest way to see this is to consider the one-dimensional case. Let $X_1 \sim \mathcal{N}(-\mu, 1)$ and $X_2 \sim \mathcal{N}(\mu, 1)$ and let them be independent of one another. Then a mixture of the two normals using weights $(\alpha, 1-\alpha)$ for $\alpha \in (0,1)$ has the distribution
$$
\varphi_m(x) = \alpha \cdot \frac{1}{\sqrt{2\pi}} e^{-\frac{(x+\mu)^2}{2}} + (1-\alpha) \cdot
\frac{1}{\sqrt{2\pi}} e^{-\frac{(x-\mu)^2}{2}} \> .
$$
The distribution of a linear combination of $X_1$ and $X_2$ using the same weights as before is, via the stable property of the normal distribution is
$$
\varphi_{\ell}(x) = \frac{1}{\sqrt{2\pi\sigma^2}} e^{-\frac{(x-(1-2\alpha)\mu)^2}{2\sigma^2}} \>,
$$
where $\sigma^2 = \alpha^2 + (1-\alpha)^2$.
These two distributions are very different, though they have the same mean. This is not an accident and follows from linearity of expectation.
To understand the mixture distribution, imagine that you had to go to a statistical consultant so that she could produce values from this distribution for you. She holds one realization of $X_1$ in one palm and one realization of $X_2$ in the other palm (though you don't know which of the two palms each is in). Now, her assistant flips a biased coin with probability $\alpha$ out of sight of you and then comes and whispers the result into the statistician's ear. She opens one of her palms and shows you the realization, but doesn't tell you the outcome of the coin flip. This process produces the mixture distribution.
On the other hand, the linear combination can be understood in the same context. The statistical consultant merely takes both realizations, multiplies the first by $\alpha$ and the second by $(1-\alpha)$, adds the result up and shows it to you.
Best Answer
I found one paper that uses the JS-Divergence of multiple distributions to estimate the hardness of a query (in the area of Information Retrieval). The paper can be found here and they themselves refer a paper called "Divergence measures based on the shannon entropy (1991)". They also give a little bit different mathematical expression for it.
As for the interpretation, they explain it as follows:
As for the Python code, I couldn't find any package that implements the JSD for more than two distributions. But there is already one quite straightforward code example on crossvalidated (see here)
.