Alternative CD formula for MOPRL of Mixed Type] Another Christoffel–Darboux Formula  for Multiple Orthogonal Polynomials of Mixed Type
Gerardo Araznibarreta
Departamento de Física Teórica II (Métodos Matemáticos de la Física), Universidad Complutense de Madrid, 28040-Madrid, Spain
Manuel Mañas

An alternative expression for the Christoffel–Darboux formula  for multiple orthogonal polynomials of mixed type is derived from the $LU$ factorization of the
moment matrix of a given measure and two sets of weights. We use the action of the generalized Jacobi matrix $J$, also responsible for the recurrence relations, on the linear forms and their duals to obtain the  result.

§ INTRODUCTION

In this paper we address a natural question that arises from the $LU$ factorization approach to multiple orthogonality  [8].  The Gauss–Borel factorization of a Hankel matrix, which plays the role of a moment matrix, leads in the classical case to a natural description of algebraic facts regarding orthogonal polynomials on the real line (OPRL) such as recursion relations and Christoffel–Darboux formula. In that case we have a chain of orthogonal polynomials $\{P_l(x)\}_{l=0}^\infty$ of increasing degree $l$. In [8] we extended that approach to the multiple orthogonality scenario, and the Gauss–Borel factorization of an appropriate moment matrix leaded to sequences of families of multiple orthogonal polynomials in the real line (MOPRL), $\big\{Q_{[\vec\nu_{1}(l);\vec\nu_{2}(l-1)]}^{(\rII,a_1(l))}\big\}_{l=0}^\infty$ and $\big\{\bar Q^{(\rI,a_1(l))}_{[\vec\nu_{2}(l);\vec\nu_{1}(l-1)]}\big\}_{l=0}^\infty$.
The recursion relations are relations constructed in terms of the linear forms in these sequences. However, the Christoffel–Darboux formula given  in Proposition <ref>, that  was re-deduced by linear algebraic means (Gauss–Borel factorization), was not expressed in terms of linear forms belonging to the mentioned sequences. This situation is rather different to the OPRL case, in that classical case the Christoffel-Darboux formula is expressed in terms of orthogonal polynomials in the sequence. The aim of this paper is to show that, within that scheme, we can deduce  an alternative but equivalent MOPRL Christoffel–Darboux formula constructed in terms of linear forms in the sequences $\big\{Q_{[\vec\nu_{1}(l);\vec\nu_{2}(l-1)]}^{(\rII,a_1(l))}\big\}_{l=0}^\infty$ and $\big\{\bar Q^{(\rI,a_1(l))}_{[\vec\nu_{2}(l);\vec\nu_{1}(l-1)]}\big\}_{l=0}^\infty$ as in OPRL situation.

§.§ Historical background

Simultaneous rational approximation starts back in 1873 when Hermite proved the transcendence of the Euler number e [24].  Later,  K. Mahler delivered  at the University of Groningen several lectures [29] where he settled down the  foundations of this theory, see also [15] and [25].  Simultaneous rational approximation when expressed in terms of Cauchy transforms leads to multiple orthogonality of polynomials. Given an interval $\Delta \subset \R$ of the real line, let ${\mathcal{M}}(\Delta)$ denote all the finite positive Borel measures with support containing infinitely many points in $\Delta$. Fix $\mu \in {\mathcal{M}}(\Delta)$, and let us consider  a system of weights $\vec w=(w_1,\ldots,w_p)$ on $\Delta$, with $p \in {\mathbb{N}}$; i.e. $w_1,\dots,w_p$ being  real integrable functions on $\Delta$ which does not change sign on $\Delta$. Fix a multi-index $\vec \nu=(\nu_1,\ldots, \nu_p) \in {\mathbb{Z}}_+^p,$ ${\mathbb{Z}}_+=\{0,1,2,\ldots\}$, and denote $|\vec \nu|=\nu_1+\cdots+\nu_p$.
Then, there exist polynomials, $A_1,\ldots, A_p$  not all identically equal to zero which satisfy the following orthogonality relations

\begin{align}\label{tipoI}
\int_{\Delta}  x^{j} \sum_{a=1}^{p} A_a(x)w_{a} (x)\d\mu (x)&=0,  & \deg
A_{a}&\leq\nu_{a}-1,&   j&=0,\ldots, |\vec \nu|-2.
\end{align}

Analogously, there exists a polynomial $B$ not identically equal to zero, such that

\begin{align}\label{tipoII}
\int_{\Delta} x^{j} B (x) w_{b} (x) \d\mu (x)&=0, & \deg B &\leq|\vec \nu|,&  j&=0,\ldots, \nu_b-1, \quad b=1,\ldots,p.
\end{align}

These are the so called multiple orthogonal polynomials of type I and type II, respectively, with respect to the combination $(\mu, \vec w, \vec \nu)$ of the measure $\mu$, the systems of weights $\vec w$ and the multi-index $\vec \nu.$ When $p=1$ both definitions coincide with standard orthogonal polynomials on the real line. Given a measure $\mu \in {\mathcal{M}}(\Delta)$ and a system of weights $\vec w$ on $\Delta$ a multi-index $\vec \nu$ is called type I or type II normal if $\deg A_a$ must equal to $\nu_a-1,$ $a=1,\ldots,p$, or $\deg B$ must equal to $|\vec \nu|-1$, respectively. When for a pair $(\mu, \vec w)$ all the multi-indices are type I or type II normal, then the pair is called type I perfect or type II perfect respectively. Multiple orthogonal of polynomials have been employed in several proofs of irrationality of numbers. For example, in [12] F. Beukers shows that Apery's proof  [10] of the irrationality of $\zeta(3)$ can be placed in the context of a combination of type I and type
II multiple orthogonality which is called mixed type multiple orthogonality of polynomials. More recently, mixed type approximation has appeared in random matrix and non-intersecting Brownian motion theories, [13], [16], [27].  Sorokin [38] studied a simultaneous rational approximation construction which is closely connected with multiple orthogonal  polynomials of mixed type. In [41]  a Riemann–Hilbert problem was found that characterizes multiple orthogonality of type I and II, extending in this way the result previously found in [23] for standard orthogonality. In [16] mixed type multiple orthogonality was analyzed from this perspective.

§.§ Perfect systems and MOPRL of mixed type

In order to introduce  multiple orthogonal polynomials of mixed type we consider two systems of weights  $\vec w_1=(w_{1,1},\dots,w_{1,p_1})$  and $\vec w_2=(w_{2,1},\dots,w_{2,p_2})$ where  $p_1,p_2\in\N$,  and two multi-indices $\vec \nu_1=(\nu_{1,1},\dots,\nu_{1,p_1})\in{\mathbb{Z}}_+^{p_1}$ and $\vec \nu_2=(\nu_{2,1},\dots,\nu_{2,p_2})\in {\mathbb{Z}}_+^{p_2}$ with $|\vec \nu_1|=|\vec \nu_2|+1$. There exist polynomials $A_1,\ldots, A_{p_1},$ not all identically zero, such that $\deg A_s < \nu_{1,s}$, which satisfy the following relations

\begin{equation}\label{orth}
\int_{\Delta} \sum_{a=1}^{p_1} A_a(x)w_{1,a} (x)w_{2,b}(x)x^{j} d\mu (x)=0, \quad  j=0,\ldots, \nu_{2,b}-1,\quad b=1,\ldots,p_2.
\end{equation}

They are called mixed multiple-orthogonal polynomials with respect to the combination $(\mu,\vec w_1,\vec w_2,\vec \nu_1,\vec\nu_2)$ of the measure $\mu,$ the systems of weights $\vec w_{1}$ and $\vec w_2$ and the multi-indices $\vec \nu_1$ and $\vec \nu_2.$  It is easy to show that finding the polynomials $A_1,\ldots, A_{p_1}$ is equivalent to solving  a system of $|\vec \nu_2|$ homogeneous linear equations  for the $|\vec \nu_1|$ unknown coefficients of the polynomials. Since $|\vec \nu_1|=|\vec \nu_2|+1$  the system always has a nontrivial solution.   The matrix of this system of equations is the so called moment matrix, and the study of its Gauss–Borel factorization will be the cornerstone of this paper. Observe that when $p_1=1$ we are in the type II case and if $p_2=1$ in  type I case. Hence in general we can find a solution of  (<ref>) where there is an $a \in \{1,\ldots,p_1\}$ such that $\deg A_a< \nu_{1,a}-1$. When given a combination $(\mu,\vec w_1,\vec w_2)$ of a measure $\mu \in {\mathcal{M}
}(\Delta)$ and systems of weights $\vec w_1$ and $\vec w_2$ on $\Delta$ if for each pair of multi-indices $(\vec \nu_1,\vec \nu_2)$ the conditions (<ref>) determine that $\deg A_a=\nu_{1,a}-1,$ $a=1, \ldots, p_1$, then we say that the combination $(\mu,\vec w_1,\vec w_2)$ is perfect. In this case

we can determine a unique system of mixed type orthogonal polynomials $\big(A_1, \ldots, A_{p_1}\big)$ satisfying (<ref>) requiring for  $a_1 \in \{1, \ldots p_1\}$ that $A_{a_1}$ monic. Following [16] we say that we have a type II normalization and denote the corresponding system of polynomials by   $A_a^{(\rII,a_1)},$ $j=1, \ldots, p_1$. Alternatively, we can proceed as follows, since the system of weights is perfect from (<ref>) we deduce that

\begin{align*}
\int x^{\nu_{1,r_1}} \sum_{a=1}^{p_1} A_a(x)w_{1,a} (x)w_{2,b}(x) \d\mu (x)\neq 0.
\end{align*}

Then, we can determine a unique system of mixed type of multi-orthogonal polynomials
$(A_1^{(\rI,a_2)},\ldots,A_{p_2}^{(\rI,a_2)})$ imposing that

\begin{align*}
\int x^{\nu_{1,a_2}} \sum_{a=1}^{p_1} A_a^{(\rI,a_2)}(x)w_{1,a} (x)w_{2,b}(x) \d\mu (x)=1,
\end{align*}

which is a type I normalization. We will use the notation $A_{[\vec\nu_1;\vec\nu_2],a}^{(\rII,a_1)}$ and $A_{[\vec\nu_1;\vec\nu_2],a}^{(\rI,a_2)}$ to denote these multiple orthogonal polynomials with type II and I normalizations, respectively.

A known illustration of perfect combinations $(\mu, \vec w_1, \vec w_2)$ can be constructed with an arbitrary positive finite Borel measure $\mu$ and  systems of weights formed with exponentials:

\begin{align}\label{exponentials}
&  (\Exp{\gamma_1x},\ldots,\Exp{\gamma_px}),& \gamma_i &\neq \gamma_j, &i &\neq j,& i,j& = 1,\ldots,p,
\end{align}

or by binomial functions

\begin{align}\label{binomials}
&((1-z)^{\alpha_1},\ldots,(1-z)^{\alpha_p}),& \alpha_i -\alpha_j &\not\in {\mathbb{Z}}, &i &\neq j,& i,j& =
\end{align}

or combining both classes, see [33]. Recently, in [22] the authors were able to prove perfectness for a wide class of systems of weights. These systems of functions, now called Nikishin systems, were introduced by E.M. Nikishin  [33] and initially named MT-systems (after Markov and Tchebycheff).

§.§ Borel–Gauss factorization and multiple orthogonality of mixed type. A remainder

Orthogonal polynomials and the theory of integrable systems has been connected in several ways in the mathematical literature. We are particularly interested in the one based in the Gauss–Borel factorization that was developed in [1]-[5], and applied further in [7]-[9].  These papers set the basis for the method we use in this paper to get an alternative Christoffel–Darboux formula for MOPRL of mixed type.

In the following we extract from  [8] the necessary material for the construction of the mentioned alternative Christoffel–Darboux formula. We introduce the moment matrix and recall how the Borel–Gauss factorization leads to multiple orthogonality. Then, we outline how the recursion relations appears by introducing a Jacobi type semi-infinite  matrix and recall the reader the Chistoffel–Darboux formula [17, 16].

§.§.§ The moment matrix

Fixed a composition $\vec n_\alpha$, $\alpha=1,2$, any given $l\in\Z_+:=\{0,1,2,\dots\}$, see [40],  determines uniquely the following non-negative
integers $k_{\alpha}(l)\in\Z_+$, $a_{\alpha}(l)\in \{1,2,\dots,p_{\alpha}\}$ and $r_{\alpha}(l)$ such that
$0\leq r_{\alpha}(l)<n_{\alpha,a_{\alpha}(l)}$ and

\begin{align}\label{i}
k_{\alpha}(l)|\vec n_{\alpha}|+n_{\alpha,1}+\dots+n_{\alpha,a_{\alpha}(l)-1}+r_{\alpha}(l), &  a_{\alpha}(l)\neq1,\\
k_{\alpha}(l)|\vec n_{\alpha}|+r_{\alpha}(l),  & a_{\alpha}(l)=1.
\end{cases}
\end{align}

We define now the monomial strings as vectors that may be understood as sequences of monomials according to the
composition $\vec n_{\alpha}$, $\alpha=1,2$,  introduced previously.

\begin{align*}
\chi_{\alpha}&:=\begin{pmatrix}
\chi_{\alpha,[0]}\\
\chi_{\alpha,[1]}\\
\vdots\\
\chi_{\alpha,[k]}\\
\vdots
\end{pmatrix} &\mbox{where}& &
\chi_{\alpha,[k]}:=\begin{pmatrix}
\chi_{\alpha,[k],1}\\
\chi_{\alpha,[k],2}\\
\vdots\\
\chi_{\alpha,[k],a_{\alpha}}\\
\vdots\\
\chi_{\alpha,[k],p_{\alpha}}
\end{pmatrix} & &\mbox{and}& &
\chi_{\alpha,[k],a_{\alpha}}:=\begin{pmatrix}
\vdots\\
\end{pmatrix}.
\end{align*}

In a similar manner for  $\alpha=1,2$ we  define the weighted monomial strings

\begin{align*}
\xi_{\alpha}&:=\begin{pmatrix}
\xi_{\alpha,[0]}\\
\xi_{\alpha,[1]}\\
\vdots\\
\xi_{\alpha,[k]}\\
\vdots
\end{pmatrix} &\mbox{where}& &
\xi_{\alpha,[k]}:=\begin{pmatrix}
\vdots\\
\end{pmatrix}.
\end{align*}

For any given $l\in\Z_+$ and $a_{\alpha}:={1,2,\dots,p_{\alpha}}$ we define

\begin{align*}
\nu_{\alpha,a_{\alpha}}(l):=\begin{cases}
k_{\alpha}(l)|\vec n_{\alpha}|+n_{\alpha,a_{\alpha}}-1, & a_{\alpha}<a_{\alpha}(l),\\
k_{\alpha}(l)|\vec n_{\alpha}|+r_{\alpha}(l), & a_{\alpha}=a_{\alpha}(l),\\
k_{\alpha}(l)|\vec n_{\alpha}|-1,& a_{\alpha}>a_{\alpha}(l).
\end{cases}
\end{align*}

Notice  that $\nu_{\alpha,a_{\alpha}}(l)$ is the hightest degree of all the monomials
of type $a_{\alpha}$ up to the component $\chi_{\alpha}^{(l)}$ included, of the monomial string. Actually

\begin{align*}
\chi_{\alpha}^{(l)}=x^{\nu_{\alpha,a_{\alpha}(l)}}.
\end{align*}

Given $l\geq 1$ and  $a_{\alpha}=1,\cdots, p_{\alpha}$ the $+$ ($-$) associated integer is  the smallest (largest) integer
$ l_{\{+,a_{\alpha}\}}$ ($l_{\{-,a_{\alpha}\}}$) such that  $l_{\{+,a_{\alpha}\}} \geq l$ ($l_{\{-,a_{\alpha}\}} \leq l$ )
and $a( l_{\{+,a_{\alpha}\}})=a_{\alpha}$ ($a( l_{\{-,a_{\alpha}}\})=a_{\alpha}$).
It can be shown that

\begin{align}\label{ai}
\begin{aligned}
l_{\{-,a_{\alpha}\}}&:=\begin{cases} k_{\alpha}(l)|\vec n_{\alpha}|+\sum_{i=1}^{a_{\alpha}}n_{\alpha,i}-1, & a_{\alpha}<a_{\alpha}(l), \\
l, &a_{\alpha}=a_{\alpha}(l),\\
k_{\alpha}(l)|\vec n_{\alpha}|-\sum_{i=a_{\alpha}+1}^{p_{\alpha}}n_{\alpha,i}-1, & a_{\alpha}>a_{\alpha}(l-1), \end{cases}\\
(  k_{\alpha}(l)+1)|\vec n_{\alpha}|+\sum_{i=1}^{a_{\alpha}-1}n_{\alpha,i},& a_{\alpha}<a_{\alpha}(l),\\
( k_{\alpha}(l)+1)|\vec n_{\alpha}|-\sum_{i=a_{\alpha}}^{p_{\alpha}}n_{\alpha,i},& a_{\alpha}>a_{\alpha}(l).
\end{cases}
\end{aligned}
\end{align}

Finally given the weighted monomials $\xi_{\vec n_{\alpha}}$, associated to the compositions $\vec n_{\alpha}$,  $\alpha=1,2$, we introduce
the moment matrix in the following manner

The moment matrix is given by

\begin{align}
\label{compact.g}
g_{\vec n_1,\vec n_2}:=\int \xi_{\vec n_1}(x)\xi_{\vec n_2}(x)^\top\d \mu(x).
\end{align}

§.§.§ Multiple Orthogonality of mixed type

For a given a perfect combination  $(\mu, \vec w_1,\vec w_2)$ we define

* The  Gauss–Borel factorization  (also known as   $LU$  factorization) of a semi-infinite moment matrix $g$,
determined by $(\mu, \vec w_1,\vec w_2)$,  is the
problem of finding the solution of

\begin{align}\label{facto}
g&=S^{-1}\bar S, & S&=\begin{pmatrix}
1&0&0&\cdots \\
\vdots&\vdots&\vdots&\ddots
\end{pmatrix}\in G_{-}, &
\bar S^{-1}&=\begin{pmatrix}
\bar S_{0,0}'&\bar S_{0,1}'&\bar S_{0,2}'&\cdots\\
0&\bar S_{1,1}'&\bar S_{1,2}'&\cdots\\
0&0&\bar S_{2,2}'&\cdots&\\
\vdots&\vdots&\vdots&\ddots
\end{pmatrix}\in G_{+},
\end{align}

where $S_{i,j},\bar S'_{i,j}\in\R.$

* In terms of these matrices we construct the polynomials

\begin{align} \label{defmops}
\end{align}

where the sum $\sum'$ is taken for a fixed $a=1,\dots,p_1$ over those $i$
such that $a=a_1(i)$ and $i\leq l$. We also construct the dual polynomials

\begin{align} \label{defdualmops}
\bar A^{(l)}_b:={\sum}'_jx^{k_2(j)}\bar S_{j,l}',
\end{align}

where the sum $\sum'$ is taken for a given $b$ over those $j$ such that $b=a_2(j)$ and $j\leq l$.

* Strings of linear forms and dual linear forms associated with  multiple ortogonal polynomials and their duals  are
defined by

\begin{align}\label{linear forms S}
Q:=  \begin{pmatrix}
\vdots
\end{pmatrix}&=S\xi_{1},&
\bar Q:=\begin{pmatrix}
\bar Q^{(0)}\\
\bar Q^{(1)}\\
\vdots
\end{pmatrix}&=(\bar S^{-1})^\top\xi_{2},
\end{align}

* The linear  forms and their duals, introduced in Definition <ref>, are given by

\begin{align}\label{linear.forms}
Q^{(l)}(x)&:= \sum_{a=1}^{p_1}A^{(l)}_{a}(x)w_{1,a}(x),&
\bar Q^{(l)}(x)&:= \sum_{b=1}^{p_2}\bar  A^{(l)}_{b}
\end{align}

* The orthogonality relations

\begin{align}\label{linear.form.orthogonality}
\begin{aligned}
\int  Q^{(l)}(x)w_{2,b}(x)x^k\d \mu(x)&=0,&0&\leq k\leq \nu_{2,b}(l-1)-1,&b&=1,\dots,p_2,\\
\int  \bar Q^{(l)}(x)w_{1,a}(x)x^k\d \mu(x)&=0,&0&\leq k\leq \nu_{1,a}(l-1)-1,&a&=1,\dots,p_1,
\end{aligned}
\end{align}

are fulfilled.

* The following multiple bi-orthogonality relations among linear forms and their duals

\begin{align}\label{biotrhoganility}
\int  Q^{(l)}(x)\bar Q^{(k)}(x)\d \mu(x)&=\delta_{l,k},& l,k\geq 0,
\end{align}

* We have the following identifications

\begin{align*}
A^{(l)}_a&=A_{[\vec\nu_{1}(l);\vec\nu_{2}(l-1)],a}^{(\rII,a_1(l))}, &
\bar A^{(l)}_b&= A^{(\rI,a_1(l))}_{[\vec\nu_{2}(l);\vec\nu_{1}(l-1)],b},
\end{align*}

in terms of multiple orthogonal polynomials of mixed type with two  normalizations $\rI$ and $\rII$, respectively.

§.§.§ Functions of the second kind

The Cauchy transforms of the linear forms (<ref>) play a crucial role in the
Riemann–Hilbert problem associated with the multiple orthogonal polynomials of mixed type
[16].  Observe that the construction of multiple orthogonal polynomials  performed so far is synthesized in the following
strings of multiple orthogonal polynomials and their duals

\begin{align}\label{mop-s}
\begin{aligned}
\A_{a}&:=  \begin{pmatrix}
\vdots
\end{pmatrix}=S\chi_{1,a},&
\bar\A_{b}&:= \begin{pmatrix}
\bar A^{(0)}_{b}\\
\bar A^{(1)}_{b}\\
\vdots
\end{pmatrix}=(\bar S^{-1})^\top\chi_{2,b},&
\end{aligned}
\end{align}

for $a=1,\dots,p_1$ and $b=1,\dots,p_2$.
In order to complete these formulae and in terms of $  \chi^*_a:=z^{-1}\chi_a(z^{-1})$
let us introduce the  following formal semi-infinite vectors

\begin{align}\label{cauchy-S}
\begin{aligned}
\Cs_b&=\begin{pmatrix}
\end{pmatrix}=\bar S\chi_{2,b}^*(z),&
\bar\Cs_a&=\begin{pmatrix}
\bar C_a^{(0)}\\\bar C_a^{(1)}\\\vdots
\end{pmatrix}=(S^{-1})^\top\chi_{1,a}^*(z),&
\end{aligned}\end{align}

for $a=1,\dots,p_1$ and $b=1,\dots,p_2$, that we call strings of second kind functions.
These objects are actually Cauchy transforms of the linear forms $Q^{(l)}$, $l \in {\mathbb{Z}}_+$, whenever the
series converge and outside the support of the measures involved.

For each  $l\in {\mathbb{Z}}_+$ the second kind functions can be expressed as follows

\begin{align}
\begin{aligned}
C_b^{(l)}(z)&=\int_\R\frac{Q^{(l)}(x)w_{2,b}(x)}{z-x}\d \mu(x),& z \in D_b^{(l)} \setminus \operatorname{supp}
(w_{1,b}\d \mu(x)),\\
\bar C_a^{(l)}(z)&=\int_\R \frac{\bar Q^{(l)}(x)w_{1,a}(x)}{z-x}\d \mu(x),&  z \in \bar D_a^{(l)} \setminus
\operatorname{supp} (w_{2,a}\d \mu(x)).
\end{aligned}
\end{align}

§.§.§ Recursion relations, a Jacobi type matrix

The moment matrix has a  Hankel type symmetry that implies the recursion relations and the Christoffel–Darboux
formula. We  consider the shift operators $\Upsilon_{\alpha}$ defined by

\begin{align}
\end{align}

Wich satisfy the following relation

\begin{align*}
\Upsilon_{\alpha}\chi_{\alpha}(x)&=x\chi_{\alpha}(x) \Longrightarrow \Upsilon_{\alpha}\xi_{\alpha}(x)=x\xi_{\alpha}(x)
\end{align*}

In terms of these shift matrices we can describe the particular Hankel symmetries for the moment matrix

The moment matrix $g$ satisfies the Hankel type symmetry

\begin{align}\label{sym2}
\Upsilon_1g  =
\end{align}

From this symmetry we see that the following is consistent

We define the matrices

\begin{align*}
J&:=S\Upsilon_1 S^{-1}=\bar S \Upsilon_2^\top\bar S^{-1}=J_++J_-,& J_+&:=(  S\Upsilon_1 S^{-1})_+, & J_-&:=(\bar S \Upsilon_2^\top\bar S^{-1})_-,
\end{align*}

where the sub-indices + and $-$denote the upper triangular and strictly lower triangular projections.

The recursion  relations follow immediately from the eigenvalue property

\begin{align}\label{rel}
J Q(x)&=x Q(x) &  \bar{Q}(x)^{\top} J&=x \bar{Q}(x)^{\top}.
\end{align}

§.§.§ Christoffel–Darboux formula

The  Christoffel–Darboux kernel is

\begin{align}
\label{def.CD}
K^{[l]}(x,y)&:=\sum_{k=0}^{l-1}Q^{(k)}(y)\bar Q^{(k)}(x)%=[\bar{Q}(x)^{\top}]^{[l]} \cdot [Q(y)]^{[l]}
\end{align}

In [17, 16] it was shown using a Riemann–Hilbert problem approach that

For  $l\geq \max(|\vec n_1|,|\vec n_2|)$ the following Christoffel–Darboux formula

\begin{align}
\begin{aligned}
\sum_{b =1}^{p_2}\bar Q^{(\rII,b)}_{[\vec\nu_2(l-1)+\vec e_{2,b};\vec \nu_1(l-1)]}(x)
Q^{(\rI,b)}_{[\vec\nu_1(l-1);\vec \nu_2(l-1)-\vec e_{2,b}]}(y)\\
&-\sum_{a=1}^{p_1}\bar Q^{(\rI,a)}_{[\vec\nu_2(l-1);\vec \nu_1(l-1)-\vec
e_{1,a}]}(x)Q^{(\rII,a)}_{[\vec\nu_1(l-1)+\vec e_{1,a};\vec \nu_2(l-1)]}(y).
\end{aligned}
\label{cd3}
\end{align}

Here $\{\vec e_{i,a}\}_{a=1}^{p_i}\subset \R^{p_i}$ stands for the vectors in the respective canonical basis, $i=1,2$.
In [8] it was given an algebraic proof of this statement not relying on analytic conditions. We refer the interested reader to [37] for a complete survey of the subject.

§ ALTERNATIVE CHRISTOFFEL–DARBOUX FORMULA FOR MULTIPLE ORTHOGONAL POLYNOMIALS OF MIXED TYPE

The result of this paper is the following

For $l\geq \max\{|\vec n_1|,|\vec n_2|\}$ the following Christofel–Darboux formula holds

\begin{align*}
(y-x)K^{[l]}(x,y)&=\smashoperator{\sum_{(i,j)\in \sigma_1[l]}}\bar{Q}(x)^{(j)} J_{j,i} Q(y)^{(i)}-
\smashoperator{\sum_{(i,j)\in\sigma_2[l]}}\bar{Q}(x)^{(j)} J_{j,i} Q(y)^{(i)}
\end{align*}

\begin{align*}
\sigma_1[l]&:=\big\{l,\dots, (l)_{\{+,r_1(a_1(l)-1)\}}\big\}\times \big\{(l-1)_{\{-,r_1(a_1(l-1)+1)\}},\dots,l-1\big\},\\
\sigma_2[l]&:=\big\{(l-1)_{\{-,r_2(a_2(l-1)+1)\}},\dots,l-1\big\}\times\big\{l,\dots,(l)_{\{+,r_2(a_2(l)-1)\}}\big\}.
\end{align*}

Splitting the eigenvalue property  (<ref>) into blocks we get

\begin{align*}
J Q(y)&=y Q(y)\Longrightarrow J^{[l]} Q(y)^{[l]}+J^{[l,\geq l]} Q(y)^{[\geq l]}=yQ(y)^{[l]} \\
\bar{Q}(x)^{\top} J&=x \bar{Q}(x)^{\top}\Longrightarrow  [\bar{Q}(x)^{\top}]^{[l]} J^{[l]}+[\bar{Q}(x)^{\top}]^{[\geq l]} J^{[\geq l,l]}=x[\bar{Q}(x)^{\top}]^{[l]}
\end{align*}

Multiply the first equation from the left by $[\bar{Q}(x)^{\top}]^{[l]}$ and the second one from the right by
$Q(y)^{[l]}$ substract both results to obtain

\begin{align*}
[\bar{Q}(x)^{\top}]^{[l]}J^{[l,\geq l]}Q(y)^{[\geq l]}-[\bar{Q}(x)^{\top}]^{[\geq l]} J^{[\geq l,l]}Q(y)^{[l]}&=
(y-x)[\bar{Q}(x)^{\top}]^{[l]}\cdot Q(y)^{[l]}\\&=(y-x)K^{[l]}(x,y)
\end{align*}

After an brief  study of the shape of $J$ we realize that even though $J^{[l,\geq l]}$ has semi-infinte
length rows, most of its elements are 0. Actually it only contains a finite number of nonzero entries that concentrate
in the lower left corner of itself. The same reasoning applies to $J^{[\geq l,l]}$. This matrix has semi infinite
length columns but again it only contains a finite number of nonzero terms concentrated in the upper right corner of itself.
Of course the number of terms involved in this expression will depend on the value of $[l]$. To be more precise we proceed as follows.  From the Euclidean division we know that for any positive integer $l\in\Z_+$ there exists unique integers $q_i, r_i$, $i=1,2$, the quotient and remainder, such that

\begin{align*}
l&=q_ip_i+r_i, & 0&\leq r_i< p_i,& i&=1,2.
\end{align*}

After a  study of the shape of $J$ we can state

For $l\geq \max\{|\vec n_1|,|\vec n_2|\}$ the  only nonzero elements of $J$ along
a given row or column are

\begin{align*}
\begin{array}{ccccccccc}
&   &       &   &         J_{(l-1)_{\{-,r_1(a_1(l-1)+1)\}},l}       &   &       &        &            \\
&   &       &   &    *    &   &       &    &                                  \\
&   &       &   & \vdots  &   &       &    &                                  \\
&   &       &   &    *    &   &       &    &                                 \\
J_{l,(l-1)_{\{-,r_2(a_2(l-1)-1)\}}} & * & \cdots & * & J_{l,l} & * & \cdots & *  & J_{l,(l+1)_{\{+,r_1(a_1(l+1)-1)\}}}\\
&   &       &   &    *    &   &       &    &                                  \\
&   &       &   & \vdots  &   &       &    &                                  \\
&   &       &   &    *    &   &       &    &                                  \\
&   &       &   &         J_{(l+1)_{\{+,r_2(a_2(l+1)-1)\}},l}    &   &       &          &                            \\
\end{array}
\end{align*}

Using this Lemma we get the desired result and the proof is complete.

In order to be more clear let us suppose that $p_1=3$ and $p_2=2$ with  $\vec n_1=(4,3,2)$ and $\vec n_2=(3,2)$. The corresponding Jacobi type matrix has the following shape

\begin{align}\label{J}
\small
J^{[12]}& J^{[12,\geq 12]}\\
J^{[\geq 12,12]} & J^{[\geq 12]}
\end{BMAT} \right)
\begin{BMAT}{cccccccccccc|ccccccccccccccc}{cccccccccccc|ccccccccccccccc}
\textcolor{blue}{*} & \boldsymbol{\textcolor{blue}{1}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \dots \\
\boldsymbol{\textcolor{blue}{*}} & \boldsymbol{\textcolor{blue}{*}} &\boldsymbol{\textcolor{blue}{1}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \dots \\
0 & \boldsymbol{\textcolor{blue}{*}} & \boldsymbol{\textcolor{blue}{\textcolor{blue}{*}}} &\boldsymbol{\textcolor{blue}{1}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \dots \\
0 & 0 & \boldsymbol{\textcolor{blue}{*}} & \boldsymbol{\textcolor{blue}{*}} & \boldsymbol{\textcolor{blue}{*}} & \boldsymbol{\textcolor{blue}{*}} &\boldsymbol{ \textcolor{blue}{*}} & \boldsymbol{\textcolor{blue}{*}} & \boldsymbol{\textcolor{blue}{*}} &\boldsymbol{\textcolor{blue}{1}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \dots \\
0 & 0 & \boldsymbol{\textcolor{blue}{*}} & \textcolor{blue}{*} & \textcolor{blue}{*} & \textcolor{blue}{*} &\textcolor{blue}{*} & \textcolor{blue}{*} & \textcolor{blue}{*} & \boldsymbol{\textcolor{blue}{*}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \dots \\
0 & 0 & \boldsymbol{ \textcolor{blue}{*}} & \boldsymbol{\textcolor{blue}{*}} & \boldsymbol{\textcolor{blue}{*}} &\textcolor{blue}{*} & \textcolor{blue}{*} & \textcolor{blue}{*} &\textcolor{blue}{*} & \boldsymbol{\textcolor{blue}{*}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \dots \\
0 & 0 & 0 & 0 & \boldsymbol{\textcolor{blue}{*}} &\textcolor{blue}{*} & \textcolor{blue}{*} & \textcolor{blue}{*} &\textcolor{blue}{*} & \boldsymbol{\textcolor{blue}{*}} & \boldsymbol{\textcolor{blue}{*}} & \boldsymbol{\textcolor{blue}{*}} & \boldsymbol{\textcolor{blue}{*}} &\boldsymbol{ \textcolor{blue}{1}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \dots \\
0 & 0 & 0 & 0 & \boldsymbol{\textcolor{blue}{*}} & \textcolor{blue}{*} & \textcolor{blue}{*} & \textcolor{blue}{*} & \textcolor{blue}{*} &\textcolor{blue}{*} & \textcolor{blue}{*} & \textcolor{blue}{*} & \textcolor{blue}{*} & \boldsymbol{\textcolor{blue}{*}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \dots \\
0 & 0 & 0 & 0 & \boldsymbol{\textcolor{blue}{*}} & \boldsymbol{\textcolor{blue}{*}} & \boldsymbol{\textcolor{blue}{*}} & \boldsymbol{\textcolor{blue}{*}} & \textcolor{blue}{*} & \textcolor{blue}{*} &\textcolor{blue}{*} & \textcolor{blue}{*} & \textcolor{blue}{*} & \boldsymbol{\textcolor{blue}{*}} & \boldsymbol{\textcolor{blue}{*}} & \boldsymbol{\textcolor{blue}{*}} &\boldsymbol{\textcolor{blue}{1}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \dots \\
0 & 0 & 0 & 0 &   0 & 0 & 0 & \boldsymbol{\textcolor{blue}{*}} & \textcolor{blue}{*} & \textcolor{blue}{*} & \textcolor{blue}{*} &\textcolor{blue}{*}& \textcolor{blue}{*} & \textcolor{blue}{*} &\textcolor{blue}{*} & \textcolor{blue}{*} & \boldsymbol{\textcolor{blue}{*}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \dots \\
0 & 0 & 0 & 0 &   0 & 0 & 0& \boldsymbol{\textcolor{blue}{*}} & \boldsymbol{\textcolor{blue}{*}} & \boldsymbol{\textcolor{blue}{*}} &\textcolor{blue}{*}& \textcolor{blue}{*} & \textcolor{blue}{*} &\textcolor{blue}{*}& \textcolor{blue}{*} & \textcolor{blue}{*} & \boldsymbol{\textcolor{blue}{*}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \dots \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \boldsymbol{\textcolor{blue}{*}} & \textcolor{blue}{*} & \textcolor{blue}{*} & \textcolor{blue}{*} & \textcolor{blue}{*} &\textcolor{blue}{*}& \textcolor{blue}{*} & \boldsymbol{\textcolor{blue}{*}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \dots \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \boldsymbol{\textcolor{blue}{*}} &\textcolor{blue}{*} & \textcolor{blue}{*} &\textcolor{blue}{*}& \textcolor{blue}{*} & \textcolor{blue}{*} &\textcolor{blue}{*}& \boldsymbol{\textcolor{blue}{*}}  &\boldsymbol{ \textcolor{blue}{*}} &\boldsymbol{\textcolor{blue}{1}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \dots \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \boldsymbol{\textcolor{blue}{*}} &\boldsymbol{\textcolor{blue}{*}} & \boldsymbol{\textcolor{blue}{*}} & \boldsymbol{\textcolor{blue}{*}} &\textcolor{blue}{*}& \textcolor{blue}{*} & \textcolor{blue}{*} &\textcolor{blue}{*} & \textcolor{blue}{*} & \boldsymbol{\textcolor{blue}{*}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \dots \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 &  0 & 0 &0 & 0 & \boldsymbol{\textcolor{blue}{*}} &\textcolor{blue}{*} &\textcolor{blue}{*} & \textcolor{blue}{*} &\textcolor{blue}{*} & \textcolor{blue}{*} &\boldsymbol{ \textcolor{blue}{*}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \dots \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 &\boldsymbol{\textcolor{blue}{*}} & \boldsymbol{ \textcolor{blue}{*}} & \boldsymbol{\textcolor{blue}{*}} & \textcolor{blue}{*} & \textcolor{blue}{*}  & \textcolor{blue}{*} & \boldsymbol{\textcolor{blue}{*}}& \boldsymbol{\textcolor{blue}{*}} & \boldsymbol{\textcolor{blue}{*}} &\boldsymbol{\textcolor{blue}{*}} &\boldsymbol{\textcolor{blue}{1}} & 0 & 0 & 0 & \dots \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 &0 & 0 & \boldsymbol{\textcolor{blue}{*}} & \textcolor{blue}{*} & \textcolor{blue}{*}  & \textcolor{blue}{*} & \textcolor{blue}{*}& \textcolor{blue}{*} & \textcolor{blue}{*} &\textcolor{blue}{*} & \boldsymbol{\textcolor{blue}{*}} & 0 & 0 & 0 & \dots \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 &0 & 0 & \boldsymbol{\textcolor{blue}{*}} & \textcolor{blue}{*} & \textcolor{blue}{*}  & \textcolor{blue}{*} & \textcolor{blue}{*}& \textcolor{blue}{*} & \textcolor{blue}{*} &\textcolor{blue}{*} & \boldsymbol{\textcolor{blue}{*}} & \boldsymbol{\textcolor{blue}{*}} & \boldsymbol{\textcolor{blue}{*}} &\boldsymbol{\textcolor{blue}{1}} & \dots \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 &\boldsymbol{\textcolor{blue}{*}} & \boldsymbol{\textcolor{blue}{*}} & \boldsymbol{\textcolor{blue}{*}}  & \boldsymbol{\textcolor{blue}{*}} & \textcolor{blue}{*}& \textcolor{blue}{*} & \textcolor{blue}{*} &\textcolor{blue}{*} & \textcolor{blue}{*} & \textcolor{blue}{*} & \textcolor{blue}{*} &\boldsymbol{ \textcolor{blue}{*}} & \dots \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \boldsymbol{\textcolor{blue}{*}} & \textcolor{blue}{*}& \textcolor{blue}{*} & \textcolor{blue}{*} &\textcolor{blue}{*} & \textcolor{blue}{*} & \textcolor{blue}{*} & \textcolor{blue}{*} & \boldsymbol{\textcolor{blue}{*}} &  \dots \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \boldsymbol{\textcolor{blue}{*}} & \boldsymbol{\textcolor{blue}{*}}& \boldsymbol{\textcolor{blue}{*}} & \textcolor{blue}{*} &\textcolor{blue}{*} & \textcolor{blue}{*} & \textcolor{blue}{*} & \textcolor{blue}{*} & \boldsymbol{\textcolor{blue}{*}}  & \dots \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 &  \boldsymbol{\textcolor{blue}{*}} & \textcolor{blue}{*} &\textcolor{blue}{*} & \textcolor{blue}{*} & \textcolor{blue}{*} & \textcolor{blue}{*} & \boldsymbol{\textcolor{blue}{*}}   & \dots \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \boldsymbol{\textcolor{blue}{*}} & \textcolor{blue}{*} &\textcolor{blue}{*} & \textcolor{blue}{*} & \textcolor{blue}{*} & \textcolor{blue}{*} &\textcolor{blue}{*} &\dots\\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \boldsymbol{\textcolor{blue}{*}} & \boldsymbol{\textcolor{blue}{*}} &\boldsymbol{\textcolor{blue}{*}} & \boldsymbol{\textcolor{blue}{*}} & \textcolor{blue}{*} & \textcolor{blue}{*} & \textcolor{blue}{*}  & \dots \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \boldsymbol{\textcolor{blue}{*}} & \textcolor{blue}{*} & \textcolor{blue}{*} & \textcolor{blue}{*} & \dots\\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 &0 &0 &0 & \boldsymbol{\textcolor{blue}{*}} & \boldsymbol{\textcolor{blue}{*}} & \boldsymbol{\textcolor{blue}{*}} & \textcolor{blue}{*}  & \dots \\
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \ddots \\
\end{BMAT}
\right),
\end{align}

where $*$ denotes a non-necessarily null real number.

In our example ($p_1=3$, $p_2=2$, $\vec n_1=(4,3,2)$ and $\vec n_2=(3,2)$)  for $[l]=[12]$ we have

\begin{align*}
(y-x)K^{[12]}(x,y)&=\left[ \sum_{i=8}^{11} \sum_{j=12}^{15} \bar{Q}(x)^{(i)}J_{i,j} Q(y)^{(j)}\right]-
\left[\sum_{i=12}^{13} \sum_{j=9}^{11} \bar{Q}(x)^{(i)}J_{i,j} Q(y)^{(j)}\right]
\end{align*}

§.§ Expressing the Jacobi type matrix in terms of factorization factors

As we have seen we can write $J$ in terms of $S$ or of $\bar S$, this means that for each term of $J$ has two different expressions, giving relations between
$S$ with $\bar S$. We are not too concerned about these relations since what we want here is the most simple expression we can get
for the elements of $J$. It is easy to realize that this is achieved if we use the expression involving $S$ in order
to calculate the upper part of $J$ and the expression involving $\bar S$ to calculate the lower part of it.  Hence, for every $J_{l,k}$ we will have expressions in terms of the
factorization matrices coefficients and the elements of their inverses –thus, in terms of the MOPRL and associated second kind functions.  The only terms from the factorization matrices (or their inverses) that will be
involved when calculating any $J_{l,k}$ are just those between the main diagonal  and the $l-|n_{1}|$ diagonal (both included)
of $S$ and those between the main diagonal and the $l+|n_2|$ diagonal (both included) of $\bar S$. And not even all of them.
As we are about to see there are three different kinds of elements in $J$. The ones
along the main diagonal, the ones along the immediate closest diagonals to the main one, and finally all the remaining diagonals.
The recursion relation coefficients $J_{k,l}$  are
ultimately related to the MOPRL and its associated  second kind functions in the following way

The elements of the recursion matrix  $J$ can be written in terms of products of the entries of the LU factorization matrices and its inverses as follows

\begin{align*}
& \begin{aligned}
J_{l,l}&=S_{l,(l-1)_{\{-,a_1(l)\}}}+S^{-1}_{(l+1)_{\{+,a_1(l)\}},l}+\sum_{\substack{a=1,\dots,p_1\\a\neq a_{1}(l)}} S_{l,(l-1)_{\{-,a\}}} S^{-1}_{(l+1)_{\{+,a\}},l},\\
&= \bar S_{l,(l+1)_{\{+,a_2(l)\}}}\bar S^{-1}_{l,l}+\bar S_{l,l}\bar S^{-1}_{(l-1)_{\{-,a_2(l)\}},l}+ \sum_{\substack{a=1,\dots,p_2\\a\neq a_{2}(l)}} \bar S_{l,(l+1)_{\{+,a\}}} \bar S^{-1}_{(l-1)_{\{-,a\}},l},
\end{aligned}
\\
J_{l,l+1}&=S^{-1}_{(l+1)_{\{+,a_1(l)\}},l+1}+ \sum_{\substack{a=1,\dots,p_1\\a\neq a_{1}(l)}} S_{l,(l-1)_{\{-,a\}}} S^{-1}_{(l+1)_{\{+,a\}},l+1}, \\
J_{l+1,l}&=   \bar S_{l+1,(l+1)_{\{+,a_2(l)\}}} \bar S^{-1}_{l,l} +\sum_{\substack{a=1,\dots,p_2\\a\neq a_{2}(l)}} \bar S_{l+1,(l+1)_{\{+,a\}}} \bar S^{-1}_{(l-1)_{\{-,a\}},l},
\end{aligned}\\
S_{l,(l-1)_{\{-,a_1\}}} S^{-1}_{(l+1)_{\{+,a_1\}},l+k} &
2&\leq k\leq (l+1)_{\{+,r_1(a_1(l+1)-1)\}}-l, \\
\bar S_{l+k,(l+1)_{\{+,a_\}}} \bar S^{-1}_{(l-1)_{\{-,a\}},l}, &  2&\leq k\leq (l+1)_{\{+,r_2(a_2(l+1)-1)\}}-l.
\end{aligned}
\end{align*}

Where, for $r,r'<p$, we have used

\begin{align*}
\sideset{}{^p}\sum_{a=r}^{r'}X_a=\begin{dcases}
\sum_{a=r}^{r'} X_a, & r\leq r',\\
\sum_{a=1}^{r'}X_a+\sum_{a=r}^pX_a, & r>r'.
\end{dcases}
\end{align*}

§ ACKNOWLEDGEMENTS

GA thanks economical support from the Universidad Complutense de Madrid  Program “Ayudas para Becas y Contratos Complutenses Predoctorales en España 2011". MM thanks economical support from the Spanish “Ministerio de Economía y Competitividad" research project MTM2012-36732-C03-01,  Ortogonalidad y aproximacion; Teoria y Aplicaciones.

[1]M. Adler and P. van Moerbeke, Group factorization, moment matrices and Toda lattices, International Mathematics Research Notices 12 (1997) 556-572.

[2] M. Adler and P. van Moerbeke, Generalized orthogonal polynomials, discrete KP and
Riemann–Hilbert  problems, Communications in  Mathematical Physics 207 (1999) 589-620.

[3] M. Adler and P. van Moerbeke, Vertex operator solutions to the discrete KP
hierarchy, Communications in Mathematical Physics 203 (1999) 185-210.

[4] M. Adler and P. van Moerbeke, The spectrum of coupled random matrices, Annals of
Mathematics 149 (1999) 921-976.

[5] M. Adler and P. van Moerbeke, Darboux transforms on band matrices, weights and
associated polynomials, International Mathematics Research Notices 18 (2001) 935-984.

[7] C. Álvarez-Fernández, U. Fidalgo and M. Mañas,
The multicomponent 2D Toda hierarchy: generalized matrix orthogonal polynomials, multiple
orthogonal polynomials and Riemann–Hilbert problems.
Inverse Problems 26 (2010) 055009   (17 pp).

[8] C. Álvarez-Fernández, U. Fidalgo and M. Mañas,
Multiple orthogonal polynomials of mixed type: Gauss-Borel factorization and the multi-component 2D Toda hierarchy.
Advances in Mathematics 227 (2011) 1451-1525.

[9] C. Álvarez-Fernández and M. Mañas,
Orthogonal Laurent polynomials on the unit circle, extended CMV ordering and 2D Toda type integrable hierarchies,
Advances in Mathematics  240 (2013) 132-193.

[10] R. Apery. Irrationalite de $\zeta(2)$ et $\zeta(3)$, Astèrisque 61 (1979) 11-13.

[12] F. Beukers, Padé approximation in number theory, Lecture Notes in Mathematics 888, Springer
Verlag, Berlin, 1981, 90-99.

[13] P.M. Bleher and A.B.J. Kuijlaars,  Random matrices with external source and multiple orthogonal
polynomials, International Mathematics Research Notices 2004 (2004), 109-129.

[15] J. Coates, On the algebraic approximation of functions, I, II, III. Indagationes Mathematicae 28
(1966) 421-461.

[16] E. Daems and A. B. J. Kuijlaars,  Multiple orthogonal polynomials of mixed type and
non-intersecting Brownian motions, Journal of Approximation Theory 146 (2007) 91-114.

[17] E. Daems and A. B. J. Kuijlaars,  A Christoffel–Darboux formula for multiple
orthogonal polynomials, Journal of Approximation Theory 130 (2004) 188-200.

[22] U. Fidalgo Prieto and G. López Lagomasino, Nikishin systems are perfect, Constructive Approximation 34  (2011) 297-356.

[23] A. S. Fokas, A. R. Its and A. V. Kitaev, The isomonodromy approach to matrix models in 2D quatum
gravity, Communications in Mathematical Physics (1992) 395-430.

[24] Ch. Hermite, Sur la fonction exponentielle, C. R. Acad. Sci. Paris 77 (1873), 18-24, 74-79,
226-233, 285-293; reprinted in his Oeuvres, Tome III, Gauthier-Villars, Paris, 1912, 150-181.

[25] H. Jager, A simultaneous generalization of the Padé table,  I-VI, Indagationes Mathematicae 26 (1964),

[27] A. B. J. Kuijlaars,    Multiple orthogonal polynomial ensembles, Contemporary Mathematics 507 (2010), 155–176

[29] K. Mahler, Perfect systems, Compositio Mathematica 19 (1968), 95-166.

[33]
E. M. Nikishin,  On simultaneous Padé approximants Matematicheskii Sbornik  113 (1980), 499–519 (Russian); English
translation in Mathematics of the USSR-Sbornik 41 (1982), 409-425.

[36]  J. A.  Shohat and J.D. Tamarkin, The problem of moments,   American Mathematical Society (1943).

[37] B. Simon, The Christoffel-Darboux Kernel,  Proceedings of Symposia in Pure Mathematics
79:“Perspectives in Partial Differential Equations, Harmonic Analysis and Applications: A Volume in
Honor of Vladimir G. Maz'ya's 70th Birthday”,  (2008) 295-336.

[38]  V.N. Sorokin, On simultaneous approximation of several linear forms, Vestnik Moskovskogo Universiteta. Seriya I. Matematika 1 (1983) 44-47.

[40] R. P. Stanley, Enumerative combinatorics, Cambridge University Press, Cambridge (1998).

[41] W. Van Assche, J. S. Geromino and A. B. J. Kuijlaars, Riemann–Hilbert problems for multiple
orthogonal polynomials in: Bustoz et al (eds.), Special Functions 2000: Current Perspectives and Future
Directions, Kluwer Academic Publishers, Dordrecht, 2001, pp 23-59.


11institutetext: Saratov State University, Astrakhanskaya str., 83, Saratov,
410012, Russia 22institutetext: Département d’Optique P. M. Duffieux, Institut
FEMTO-ST, Université Bourgogne-Franche-Comté CNRS UMR 6174, Besançon, France
33institutetext: State Medical University, B. Kazachaya str., 112, Saratov,
410012, Russia 44institutetext: Physics Department, Humboldt University,
Newtonstrasse 15, 12489 Berlin, Germany 55institutetext: Potsdam Institute for
Climate Impact Research, Telegrafenberg A31, 14473 Potsdam, Germany

# Non-invasive analysis of blood-brain barrier permeability based on wavelet
and machine learning approaches

Nadezhda Semenova 1122<EMAIL_ADDRESS>Konstantin Segreev 11 Andrei
Slepnev 11 Anastasia Runnova 1133 Maxim Zhuravlev 1133 Inna Blokhina 11
Alexander Dubrovsky 11 Oxana Semyachkina-Glushkovskaya 1144 Jürgen Kurths
114455

###### Abstract

The blood-brain barrier plays a decisive role in protecting the brain from
toxins and pathogens. The ability to analyze the BBB opening (OBBB) is crucial
for the treatment of many brain diseases, but it is very difficult to
noninvasively monitor OBBB. In this paper we analyse the EEG series of healthy
rats in free behaviour and after music-induced OBBB. The research is performed
using two completely different methods based on wavelet analysis and machine
learning approach. Both methods enable us to recognize OBBB and are in a good
agreement with each other. The comparative analysis was carried out using
F-measures and ROC-curves.

##  1 Introduction

The blood-brain barrier (BBB) is a highly selective barrier, which is formed
by microvascular endothelial cells surrounded by pericytes and perivascular
astroglia. It controls the penetration of blood-borne agents into the brain or
the release of metabolites and ions from the brain tissue to blood Abbott
(2005); Abbott et al. (2006, 2010). Therefore, the BBB plays a vital role in
protecting the brain against pathogens and toxins. The BBB disruption is
associated with aging, dementia, multiple sclerosis, Alzheimer’s disease,
stroke, brain trauma, infection and tumors Rosenberg (2012). The permeability
of BBB can be varied by neuroendocrine regulation and may play a protective
role in injury and stroke Pan et al. (1997, 1999, 2001, 2006).

The ability to analyze the BBB opening (OBBB) is crucial for the treatment of
brain diseases, and it is very difficult to noninvasively monitor OBBB
Chassidim et al. (2013); Heye et al. (2014). Therefore, there is a strong need
to create many save methods of the assessment of the BBB permeability in
clinical practice. There has been substantial progress over the past few years
Chassidim et al. (2013); Heye et al. (2014). Magnetic resonance tomography
(MRI) is an often used technique for the monitoring of OBBB Heye et al.
(2014). However, MRI is bulky and cannot be used bedside. The MRI requests the
use of contrast agents, which can be even toxic Rogosnitzky and Branch (2016);
Kaller and An (2020). The latter limits its continuous application and usage,
especially in the case of children and patients with kidney pathology
Elbeshlawi and AbdelBaki (2018); Perazella (2008). Therefore, the development
of novel promising real-time, bedside, non-invasive, label-free, economically
beneficial and readily applicable methods is of highest actual importance, and
solving this problem would open a novel era in effective diagnosis and therapy
of brain diseases that cause acute and chronical OBBB.

In this paper we consider two different approaches for analysing the OBBB by
using electroencephalogram (EEG) time series of healthy rats. The first method
is based on the wavelet analysis of EEG signals. This approach is widely used
in many biomedical papers Unser and Aldroubi (1996); Subasi (2007), partially
in epilepsy research Adeli et al. (2003); Tzallas et al. (2009), and sleep
staging Bajaj and Pachori (2013); Tripathy et al. (2020), which confirms its
applicability in these tasks. Nowadays, the wavelet analysis is one of the
main tools for processing data in real time and creating brain-computer
interface (BCI) devices Djamal and Putra (2020); Nguyen et al. (2015);
Hettiarachchi et al. (2014).

The second method uses a machine learning approach to analyze such EEG
signals. Recently, machine learning has been increasingly used in brain
activity analysis Davatzikos et al. (2005); Palaniappan and Mandic (2007); Wei
et al. (2018); yao Tian et al. (2020). The application of machine learning to
recognize features of brain activity leads to the use of artificial neural
networks (ANNs) of different types and topologies. Basically, ANNs do not
analyze the EEG signal in its pure form, but use some characteristics based on
the EEGs. For example, power spectral density Ronzhina et al. (2012),
statistical features Hassan et al. (2015), wavelet transform Hassan and
Bhuiyan (2017), combinations of energies of different frequency bands Hsu et
al. (2013) etc. Typically, deep neural networks (DNNs) are used for
classification and recognition tasks. As one of examples, the automatic
scoring of sleeping stage was recently provided by applying a convolutional
neural network Fernandez-Blanco et al. (2020). In this work we use DNN with
feedforward coupling.

Despite the fact that both methods are based on completely different types of
analysis, they show a good ability in recognizing the OBBB. BBB permeability
can be evaluated in real time, searching for the EEG signal features, which
can be found out by using both wavelet analysis and machine learning
approaches.

##  2 Experiment design and data recording

The experiments were conducted on the same five adult male Wistar-Kyoto rats
(250-280 g) in awake state and during opened BBB (OBBB). All procedures and
experiments were done in accordance with the “Guide for the Care and Use of
Laboratory Animals” Albus (2012). The experimental protocols were approved by
the Local Bioethics Commission of the Humboldt University and the Saratov
State University. The animals were kept in a light/dark environment with the
lights on from 8:00 to 20:00 and fed ad libitum with standard rodent food and
water. The ambient temperature and humidity were maintained at $24.5\pm
0.5^{\circ}C$ and 40%–60%, respectively.

Ten days before the experiments, an optical window, a polyethylene catheter,
and two epidural EEG electrodes for chronic recordings of electrocorticogram
(ECoG) in free behavior were implanted; after surgery, animals were housed
individually. The first EEG1-recording was recorded over an hour from 1 to 2
p.m. 10 days after the electrode implantation operation.

The day after the baseline EEG1-recording, the rats were underwent to
intermittent music. The detailed description of the protocol of music-induced
OBBB is described in Semyachkina-Glushkovskaya et al. (2020). To produce the
music (70-90-100 dB and 11-10000 Hz, Scorpions “Still loving you”) we used
loudspeaker (ranging of sound intensity 0-130 dB, frequencies 63-15000 Hz; 100
V, Yerasov Music Corporation, Saint Petersburg, Russia). The repetitive music
exposure was performed using a sequence of: 60 sec – music on and then 60 sec
– music off during 2h. The sound level was measured directly in a cage of
animals using the sound level meter (Megeon 92130, Russia). During the 60 min
music-off, we performed in vivo real time fluorescent microscopy of OBBB for
Evans Blue (EB) dye (i.v.), which we detected as bright fluorescence around
the cerebral microvessels in awake behaving rats. When we observed the optical
signs of EB leakage, EEG recordings in rats were performed during 60 min
(EEG2). The choice of time of EEG recordings was related to our previous
results discovering the window of music-induced OBBB that is 60-180 min after
sound exposure Semyachkina-Glushkovskaya et al. (2017); Zinchenko et al.
(2019); Semyachkina-Glushkovskaya et al. (2020). After EEG recordings, the
rats decapitated and their brains were collected for confirmation of OBBB for
ED using confocal microscopy. The EEG recordings were collected and analyzed
using wavelet and machine learning approaches.

A two-channels cortical EEG/one channel electromyogram (EMG) (Pinnacle
Technology, Taiwan) were recorded as follows. The rats were implanted two
silver electrodes (tip diameter 2–3 $\mu$m) located at a depth of 150 $\mu$m
in the coordinates (L: 2.5 mm and D: 2 mm) from Bregma on either side of the
midline under inhalation anesthesia with 2% isoflurane at 1L/min
$N_{2}O/O_{2}$ – 70:30. The head plate was mounted and small burr holes were
drilled. Afterward, EEG wire leads were inserted into the burr holes on one
side of the midline between the skull and the underlying dura. EEG leads were
secured with dental acrylic. An EMG lead was inserted in the neck muscle.
Ibuprofen (15 mg/kg) for the relief of postoperative pain was provided in
their water supply for two to three days prior to surgery and for three or
more days post-surgery. The rats were allowed 10 days to recover from surgery
prior to beginning the experiment.

All EEG recordings were supplemented with synchronous videography, on the
basis of which a neurophysiologist marked out the states of behavioral sleep
and wakefulness of animals. Since the standard sleep staging rules for rats
are not available currently, we referred the visual scoring criteria from
these studies Xie et al. (2013); Hablitz et al. (2019). Wakefulness, NREM, and
REM were visually scored in 10 s epochs. EEG activity was measured and
compared in rats in awake state, during sleep, and after OBBB on the same
rats. Wakefulness was defined as a desynchronized EEG with low amplitude and
high-frequency dynamics ($>$10%, 8-12 Hz) and relatively high-amplitude EMG.
NREM sleep was recognized as synchronized activity with high amplitude, which
is dominated by low-frequency delta waves (0–4 Hz) comprising $>$30% of EEG
waveforms/epoch and a lower amplitude EMG. REM was identified by the presence
of theta waves (5–10 Hz) comprising $>$20% of EEG waveforms/epoch with a low
EMG amplitude.

Thus, for each of the animals the EEG signals were recorded before and after
music-induced OBBB – EEG1 and EEG2, respectively.

##  3 Recognition of OBBB using time-frequency analysis

Today, one of the generally accepted nonlinear methods for processing and
detecting of oscillatory activity in biomedical signals is continuous wavelet
transformation (CWT) Lachaux and et al. (2002); Feng et al. (2017); Khare et
al. (2020). CWT is sufficiently resistant to abrupt changes in the frequency
composition of the analyzed experimental signals, which makes it possible to
adequately analyze rather short time intervals of highly nonstationary
signals. The CWT for an arbitrary signal $x(t)$ in the general form is defined
as follows:

$W(s,t_{0})=\int\limits_{-\infty}^{+\infty}x(t)\psi_{s,t_{0}}^{*}(t)dt,$ (1)

where $\psi_{s,t_{0}}(t)$ is the basic complex function, $s$ is the time scale
defining the width of the wavelet, the symbol “*” denotes complex conjugation.
Note that the time scales $s$ of the CWT allow a transition to the classical
frequencies $f$ of the Fourier spectrum. Therefore, for convenience and
simplicity of the interpretation of the results, we consider the results in
the traditional plane $(f,t_{0})$.

The basic Morlet function $\psi_{s,t_{0}}(t)$ is often used to analyze the
activity of the brain Aouinet and Adnane (2014); Ovchinnikov et al. (2011). By
analogy with the Fourier power spectrum $F(f)$ for $f$-frequencies we can
estimate the instantaneous frequency distribution of the CWT-energy

$E(f,t_{0})=|W(f,t_{0})|^{2}.$ (2)

CWT allows to analyse an experimental EEG signal $x(t)$ simultaneously in the
frequency and time domains. Moreover, the CWT excess property is well known,
which is used to search for the fine time-frequency structure of experimental
time series, which can be only poorly detected by conventional spectral
methods Montefusco (2014); Maksimenko et al. (2018a, b); Hramov et al. (2015).

In this paper to analyze the oscillatory structure of the EEG, we use the CWT
“skeleton” method Cavalier and O’Hagan (2017); Runnova et al. (2017); Grubov
et al. (2017a, b); Hramov et al. (2017), which is based on a simplification of
the entire surface $(f;t)$ of the experimental signal $x(t)$ by estimation
only the main part of oscillatory activity. For each moment of time $t_{0}$,
it is necessary to search for the maxima in the instantaneous CWT spectrum
$E(f,t)$ (2):

$\forall f\in\delta f_{max},E(f_{max},t_{0})\geq E(f,t_{0}),$ (3)

where $f_{max}=sc^{1}$ is the CWT skeleton, $\delta f_{max}$ is a $\delta$ \-
region of the $f_{max}$ frequency component. The first skeleton $f_{1}$ is the
global maximum in the surface $W(f,t_{0})$, i. e. the value of the frequency
in the signal spectrum, which correlates with the maximum of the oscillational
energy. Further, for each moment of time, we detect a local maxima
$sc^{1}>sc^{2}>\dots>sc^{j}>\dots>sc^{n_{p}}$, where $n_{p}$ is the amount of
skeletons.

Further, we limit the consideration of the oscillatory structure in EEG
signals to the frequency band $\Delta f[1;2,5]$ Hz. The choice of this “slow”
frequency interval is due to the evidence in the literature on the presence of
biomarkers in brain slow activity reflecting the BBB permeability changes
Pavlov et al. (2020); yao Tian et al. (2020). Thus, at each time moment
$t_{0}$ we estimate which skeleton’s number $N_{0}$ falls into a certain
frequency range $\Delta f$. And next we calculate the average number of $N(t)$
skeletons falling into the band $\Delta f$ over the time interval $\Delta
t=50$ sec.

In Fig. 1(a) the results of estimating the number of $N(t)$ skeletons for one
experimental animal are presented. We observe that the number of patterns in
the low frequency range $\Delta f_{1}$ increases after a powerful auditory
impact (EEG2). After evaluating the distribution of the number of patterns for
records EEG1 and EEG2, their mean values differ $\langle
N\rangle_{EEG1}\neq\langle N\rangle_{EEG2}$ (see Fig. 1(b)). A similar
situation is observed in the slow oscillatory activity of the EEG for all
experimental animals. Next, we use the detectable value $\langle
N\rangle_{EEG2}$ as the threshold value that characterizes the moment of BBB
permeability changes.

For each experimental animal, individual threshold values were obtained
$\langle N\rangle_{EEG2}=\\{0,31863;$ $0,45614;$ $0,462575;$ $0,485995;$
$0,52011\\}$. Thus, we estimated the EEG-dynamics of the brain activity based
on the $N(t)$-dynamic. Generally, the value of dependence $N(t)$ in
EEG2-record exceeded the threshold $\langle N\rangle_{EEG2}$ throughout first
third of the recording. Next, this value $N(t)$ irregularly reduced towards
the end of the EEG2-record, when, apparently, the influence of the sound
impact was compensated by the neurophysiological system, and the BBB-
permeability returned to its usual value. Based on this, we assumed that the
excess of the $N(t)$ patterns number of the $\langle N\rangle_{EEG2}$
threshold value correlated with an increase in the BBB permeability. So, when
the corresponding threshold value $N(t)\geq\langle N\rangle_{EEG2}$ is
exceeded, we assumed that the permeability of the BBB was increased.

Figure 1: Results of evaluating the wavelet-based characteristics for animal #
3. Panel (a) shows the time-dependence of the $N$ average number of CWT-
skeletons in the frequency range $\Delta f_{1}[1;2,5]$ Hz. The red line is the
calculation result for EEG1, the green line corresponds to the characteristic
for EEG2. The panel (b) represents the probability distributions of $N$
skeletons. Red color shows the distribution for EEG1 (no influence on the BBB-
permeability), and green – for EEG2 (after auditory impact). The arrows
indicate the mean values for each distribution.

##  4 ANN based method

In order to estimate the permeability of the BBB, a deep neural network is
used. The neural network was trained using an open-source library in Python
Chollet et al. (2015). The use of EEG signals in their pure form is a non-
trivial task for the ANN, since it requires a large number of neurons and
layers. This dimension problem can be solved by specially prepared statistical
characteristics calculated from these EEGs. There is a number of papers on
sleep recognition, in which the EEG signals were not used by ANN, but their
averages, variations and deviations from the mean Ronzhina et al. (2012); Hsu
et al. (2013); Zhu et al. (2014); Hassan et al. (2015); Hassan and Bhuiyan
(2017). In the same manner here we use a signal-to-noise ratio (SNR) as such a
statistical characteristic. It is the ratio between the mean of the signal
$\mu(\cdot)$ and its standard deviation $\sigma(\cdot)$:

${SNR}(x)=\frac{\mu(x)}{\sigma(x)}.$ (4)

Each point of this characteristic is prepared for some averaging window. The
window size is fixed to 30 seconds. In order to make the SNR realization
smoother, the calculations were performed with a slowly shifting time window
with a 1-second step. This leads to a new temporal determination of SNR with a
time step of 1 sec, where each point corresponds to the SNR value obtained
from a thirty-second window of EEG data. Figure 2 shows the EEG of one of the
animals (panels a, c) and the calculated SNR characteristics (panels b, d).

To get the training SNR examples, we used EEG signals before (EEG1) and after
(EEG2) artificially increased permeability of the BBB Semyachkina-
Glushkovskaya et al. (2018). During the training, the fragments of SNR for the
rat after music stimulation were marked as the expected response “1” (i.e. BBB
is open), while fragments of rats before this stimulation corresponded to the
response “0” (i.e. BBB is closed). Thus, the ANN response to an unknown signal
can be interpreted as a degree of the averaged “permeability” of BBB over the
considered period of time.

Figure 2: Examples of EEG signals (a,c) for animal # 1 and corresponding SNR
characteristics (eq. 4) shown in panels (b,d).

The input layer of our ANN consists of 90 neurons. Each time the network
receives a 90-second SNR realization as an input. A smaller number of input
neurons leads to slower training and a sharp jump into overtraining. A larger
number leads to inaccuracies in the temporal marking, since it gives the
response delay of more than 1.5 minutes. Also a few DNN configurations with
constant numbers of neurons in the hidden layers 200, 500, 1000 have been
considered. A network with a variable number of neurons in hidden layers
(500$\times$200$\times$500$\times$200) has shown the optimal results during
training and testing. The network scheme and the training process are shown in
Fig. 3. All the layers have the sigmoid activation function
$f(x)=1/(1-e^{-x})$.

Figure 3: Schematic diagram of considered DNN (panel a) and the dynamics of
accuracy during training depending on the epoch number (panel b)

For each of the rats the ANN was trained on the data of the remaining rats.
This increased the number of training and testing implementations. Each EEG
recording was done using two channels. The SNR and ANN markup were prepared
for both signals separately. To increase accuracy, the network responses for
both channels were multiplied. Thus, if the network recognises the inputs from
both channels are recognized simultaneously as “1” (BBB is open), then the
final answer is “1”. If one channel leads to the answer “0” and the second one
leads to the answer “1”, then the final response is “0”. Since the network
produces a real number from 0 to 1, the multiplication of responses from two
channels can lead to an increase in intermediate values. To eliminate this,
the threshold 0.5 has been introduced. Since then, the answers $y\geq 0.5$ are
regarded as “1” (BBB is open), while responses $y<0.5$ complies with “0” (BBB
is closed). Figure 4 shows a temporal evolution of ANN responses for an animal
with an artificially open barrier, whose signals were shown in Fig. 2(a,b).
The gray line represents the result of multiplying the signals from both
channels, and red pints show the final answer with the threshold 0.5.

The ANN found 32% of the data similar to the open barrier for implementations
after music-induced OBBB. In the case of free behaviour, this percentage was
24%. These ratios are averaged over all considered animals. In order to
exclude the peculiarities of the training, the networks have been trained 5
times for each animal with an accuracy of at least 96% (approximately 250
epochs). Table 1 shows the percentage of fragments recognized as OBBB for one
animal after 5 trainings. Afterward, the general percentage of fragments,
recognized as OBBB, has been averaged over all animals and over all five
trainings. Figure 5 shows examples of SNR obtained using fragments of EEG
signals recognized as signals with open BBB (top panels) and closed BBB
(bottom panels).

Figure 4: Temporal evolution of ANN responses for animal # 1 after music-
induced OBBB. The multiplication of network responses for both EEG channels is
shown in gray. Red points correspond to the binary classification after
introducing the threshold 0.5. Green areas show time intervals in which both
markings coincide: (1) based on ANN answers and (2) wavelet method from Sect.
3.

##  5 Comparison of two independent methods of data processing

Figure 4 shows the intersection of both methods from Sect. 3 and 4 by green
background color. We find that both methods are in a good agreement with each
other. To test this numerically, the so-called F-measure Van Rijsbergen (1979)
is applied, which is often used for statistical analysis of binary
classification. It is the harmonic mean between precision and recall. The
precision is the number of correctly identified positive results divided by
the number of all positive results, including those not identified correctly,
and the recall is the number of correctly identified positive results divided
by the number of all samples that should have been identified as positive. For
such an assessment, one can choose one of the methods as true, for example,
the wavelet method, and then count the number of matches of both methods for
the answer “yes” (true-positive, TP), the number of matches for the answer
“no” (true-negative, TN), as well as mismatches (false-positive, FP, and
false-negative, FN). Based on this, the F-measure is calculated as follows:

$P=\frac{TP}{TP+FP},R=\frac{TP}{TP+FN},F^{+}=2\frac{PR}{P+R}$ (5)

Figure 5: Examples of input SNR implementations on which the neural network
exhibited a similarity to OBBB (top panels) and closed barrier (bottom
panels).

Thus, the F-measure $F^{+}$ is an evaluation of the correctness for the
answers “yes”. If we replace all letters “P” by “N” in Eq. (5), and vice
versa, then this will be an estimate of the correctness for the answer “no”,
which we denote as $F^{-}$. Table 1 shows the results of calculating the
F-measure for five ANN trainings, as well as the averaged characteristics.
Thus, according to the calculations of the F-measure, we can conclude that the
intersection of 35 % in the answers that the barrier is open, and 75 % in the
answers about the closed barrier. These values, averaged over all animals, are
28 % and 70 %.

Training: | Tr.1 | Tr.2 | Tr.3 | Tr.4 | Tr.5 | Averaged
---|---|---|---|---|---|---
% of open barrier in realization with an open barrier | 35 | 22 | 21 | 38 | 27 | 28.5
% of open barrier in realization with free behaviour | 24 | 18 | 13 | 24 | 20 | 19.8
F-measure on responses “BBB is open”, $F^{+}$ | 0.38 | 0.35 | 0.32 | 0.38 | 0.34 | 0.35
F-measure on responses “BBB is closed”, $F^{-}$ | 0.73 | 0.77 | 0.79 | 0.72 | 0.75 | 0.75
Table 1: Percentage of fragments recognized as OBBB and not, and F-measures of
animal # 1 for five ANN trainings on fragments of the other animals.

Another measure that allows one to assess the quality of the classification is
the receiver operating characteristic (ROC-curve). The ROC curve is created by
plotting the true positive rate ($TPR=\frac{TP}{TP+FN}$) against the false
positive rate ($FPR=\frac{FP}{FP+TN}$) at various threshold settings Fawcett
(2006); Powers (2011). The classifier is considered to be good if the ROC-
curve is above the diagonal (see Fig. 6). However, the degree of success of
the classifier is roughly estimated by how far this curve is above the
diagonal, and, accordingly, how large is the area below it. Figure 6 shows the
ROC curves for different threshold values in cases of artificially opened BBB
(orange color) and free behaviour (blue color), when the rat can either sleep
or be awake. The markings obtained by the method, described in Sect. 3, were
used as a verification data. Both curves are above the diagonal, which
indicates the similarity of the methods. In addition, the figure shows the
results of calculating the F-measures for both implementations with variation
of the threshold parameter. The points corresponding to the $0.5$ threshold
are highlighted by red. As follows from the graphs of F-measures, the
threshold $0.5$ is optimal for positive and negative answers. An increase in
the threshold leads to an increase in false positives, and a decrease leads to
a loss of sensitivity when the network says that the barrier is always closed.

Figure 6: Statistical comparison of both methods using ROC-curve (panel a),
and F-scores (eq. 5) for positive (b) and negative (c) responses.

##  6 Conclusions

This work is devoted to the investigation of sustained changes in brain
activity after music-induced OBBB. BBB permeability is analyzed using two
different methods – wavelet estimation of the oscillatory characteristics of
the EEG and machine learning using SNR characteristics. The wavelet analysis
demonstrated an increase in the number of arising oscillations in the low-
frequency region (1 – 2.5 Hz), which is universal for all animals. However,
the threshold for the increase in oscillatory activity in these frequencies is
highly individual for each animal.

The second method is based on the search of special features of SNR
realisation which are common for animals with OBBB. The use of a machine
learning approach allows us to demonstrate a high level of universality of the
revealed “portrait” of the electrical activity of the brain after an
artificial increase in the permeability of the BBB.

Despite the fact that both methods are based on different characteristics,
they show a good overlap both for animals with artificially open BBB and for
animals in free behaviour without any influences. Quantification of the
methods agreement was based on ROC curves and F-measures of positive and
negative responses. Both methods recognized a non-zero probability of opening
the barrier for animals with normal behaviour. It is possible that this effect
is related to the constancy of brain activity patterns during the
physiological state of sleep. A number of authors have linked spontaneous
changes in BBB permeability to certain stages of sleep, promoting a repair and
clearance of brain tissue Semyachkina-Glushkovskaya et al. (2017); Zinchenko
et al. (2019); Pavlov et al. (2020).

In the future, we plan to increase of the experimental groups and to include
instrumental methods of optical control of the BBB permeability. The
preliminary analysis of brain activity during the sleeping state of the
animals is of our interest. Probably, some sleep stages would be accompanied
with the brain activity, similar to the results of music-induced OBBB.

##  7 Acknowledgements

This work has been supported by the RF Government Grant No. 075-15-2019-1885
in part of the biological interpretation. In the part of the development of
numeric method of data analysis this work has been supported by the Council
for Grants of the President of the Russian Federation for the State Support of
Young Russian Scientists (project no. MD-645.2020.9).

## References

  * Abbott (2005) N. Abbott, Cell. Mol. Neurobiol. 25, 5 (2005).
  * Abbott et al. (2006) N. Abbott, L. Rönnbäck, and L. Hansson, Nature Rev. Neurosci. 7, 41 (2006).
  * Abbott et al. (2010) N. Abbott, A. A. Patabendige, D. E. Dolman, S. Yusof, and D. Begley, Neurobiology of Disease 37, 13 (2010).
  * Rosenberg (2012) G. Rosenberg, J. Cereb. Blood. Flow. Metab. 32, 1139 (2012).
  * Pan et al. (1997) W. Pan, W. Banks, and A. Kastin, Exp Neurol 146, 367 (1997).
  * Pan et al. (1999) W. Pan, A. Kastin, R. Bell, and R. Olson, J. Neurosci. 19, 3649 (1999).
  * Pan et al. (2001) W. Pan, A. Kastin, L. Gera, and J. Stewart, Neurosci Lett 307, 25 (2001).
  * Pan et al. (2006) W. Pan, Y. Ding, Y. Yu, H. Ohtaki, T. Nakamachi, and A. Kastin, Exp Neurol 198, 222 (2006).
  * Chassidim et al. (2013) Y. Chassidim, R. Veksler, S. Lublinsky, G. Pell, A. Friedman, and I. Shelef, Fluids Barriers CNS 10, 9 (2013).
  * Heye et al. (2014) A. K. Heye, R. D. Culling, M. C. Valdés Hernández, M. J. Thrippleton, and J. M. Wardlaw, NeuroImage Clin 6, 262 (2014).
  * Rogosnitzky and Branch (2016) M. Rogosnitzky and S. Branch, Biometals 29 (2016).
  * Kaller and An (2020) M. Kaller and J. An, _Contrast Agent Toxicity_ (Treasure Island (FL): StatPearls Publishing, 2020).
  * Elbeshlawi and AbdelBaki (2018) I. Elbeshlawi and M. S. AbdelBaki, Pediatr Neurol 86, 27 (2018).
  * Perazella (2008) M. A. Perazella, Curr. Drug. Saf. 3, 67 (2008).
  * Unser and Aldroubi (1996) M. Unser and A. Aldroubi, PROCEEDINGS OF THE IEEE 84, 626 (1996).
  * Subasi (2007) A. Subasi, Expert Systems with Applications 32, 1084 (2007).
  * Adeli et al. (2003) H. Adeli, Z. Zhou, and N. Dadmeh, Journal of Neuroscience Methods 123 (2003).
  * Tzallas et al. (2009) A. Tzallas, M. Tsipouras, and D. Fotiadis, IEEE TRANSACTIONS ON INFORMATION TECHNOLOGY IN BIOMEDICINE 13 (2009).
  * Bajaj and Pachori (2013) V. Bajaj and R. Pachori, Comput. Methods Programs Biomed 112, 320 (2013).
  * Tripathy et al. (2020) R. Tripathy, S. Ghosh, P. Gajbhiye, and U. Acharya, Entropy 22 (2020).
  * Djamal and Putra (2020) E. C. Djamal and R. D. Putra, TELKOMNIKA Telecommunication, Computing, Electronics and Control 18, 2748 (2020).
  * Nguyen et al. (2015) T. Nguyen, A. Khosravi, D. Creighton, and S. Nahavandi, Expert Systems with Applications 42, 4370 (2015).
  * Hettiarachchi et al. (2014) I. T. Hettiarachchi, T. T. Nguyen, and S. Nahavandi, in _International Conference on Neural Information Processing_ (Springer, 2014), pp. 519–526.
  * Davatzikos et al. (2005) C. Davatzikos, K. Ruparel, Y. Fan, D. Shen, M. Acharyya, J. Loughead, R. Gur, and D. Langleben, NeuroImage 28, 663 (2005), ISSN 1053-8119, URL http://www.sciencedirect.com/science/article/pii/S1053811905005914.
  * Palaniappan and Mandic (2007) R. Palaniappan and D. P. Mandic, IEEE Transactions on Pattern Analysis and Machine Intelligence 29, 738 (2007).
  * Wei et al. (2018) J. Wei, T. Chen, C. Li, G. Liu, J. Qiu, and D. Wei, Frontiers in Human Neuroscience 12, 422 (2018), ISSN 1662-5161, URL https://www.frontiersin.org/article/10.3389/fnhum.2018.00422.
  * yao Tian et al. (2020) Z. yao Tian, L. Qian, L. Fang, X. hua Peng, X. hu Zhu, M. Wu, W. zhi Wang, W. han Zhang, B. qi Zhu, M. Wan, et al., Neuroscience 436, 170 (2020), ISSN 0306-4522, URL http://www.sciencedirect.com/science/article/pii/S0306452220300798.
  * Ronzhina et al. (2012) M. Ronzhina, O. Janousek, J. Kolarova, M. Novakova, P. Honzik, and I. Provaznik, Sleep medicine reviews 16, 251 (2012).
  * Hassan et al. (2015) A. Hassan, S. Bashar, and M. Bhuiyan, in _International Conference on Advances in Computing, Communications and Informatics (ICACCI)_ (2015), pp. 2238–2243.
  * Hassan and Bhuiyan (2017) A. Hassan and M. Bhuiyan, Computer methods and programs in biomedicine 140 (2017).
  * Hsu et al. (2013) Y.-L. Hsu, Y.-T. Yang, J.-S. Wang, and C.-Y. Hsu, Neurocomputing 104, 105 (2013).
  * Fernandez-Blanco et al. (2020) E. Fernandez-Blanco, D. Rivero, and A. Pazos, Neurocomputing 410, 220 (2020).
  * Albus (2012) U. Albus, Laboratory Animals 46, 267 (2012).
  * Semyachkina-Glushkovskaya et al. (2020) O. Semyachkina-Glushkovskaya, A. Esmat, D. Bragin, O. Bragina, A. A. Shirokov, N. Navolokin, Y. Yang, A. Abdurashitov, A. Khorovodov, A. Terskov, et al., bioRxiv p. 2020.10.03.324699 (2020), URL https://app.dimensions.ai/details/publication/pub.1131448754 and https://www.biorxiv.org/content/biorxiv/early/2020/10/05/2020.10.03.324699.full.pdf.
  * Semyachkina-Glushkovskaya et al. (2017) O. Semyachkina-Glushkovskaya, A. Abdurashitov, A. Dubrovsky, D. Bragin, O. Bragina, N. Shushunova, G. Maslyakova, N. Navolokin, A. Bucharskaya, V. Tuchind, et al., Journal of Biomedical Optics 22, 121719 (2017).
  * Zinchenko et al. (2019) E. Zinchenko, N. Navolokin, A. Shirokov, B. Khlebtsov, A. Dubrovsky, E. Saranceva, A. Abdurashitov, A. Khorovodov, A. Terskov, A. Mamedova, et al., Biomedical Optics Express 10, 4003 (2019).
  * Xie et al. (2013) L. Xie, H. Kang, Q. Xu, M. Chen, Y. Liao, M. Thiyagarajan, J. O’Donnell, D. Christensen, C. Nicholson, J. Iliff, et al., Science 342 (2013).
  * Hablitz et al. (2019) L. M. Hablitz, H. S. Vinitsky, Q. Sun, F. F. Stæger, B. Sigurdsson, K. N. Mortensen, T. O. Lilius, and M. Nedergaard, Science Advances 5 (2019), https://advances.sciencemag.org/content/5/2/eaav5447.full.pdf, URL https://advances.sciencemag.org/content/5/2/eaav5447.
  * Lachaux and et al. (2002) J. P. Lachaux and et al., Neurophysiol. Clin. 32, 157 (2002).
  * Feng et al. (2017) Q. Feng, M. Zhang, Y. Zhang, N. Jiang, and J. Zhang, Journal of Medical Imaging and Health Informatics 7, 928 (2017).
  * Khare et al. (2020) S. K. Khare, V. Bajaj, S. Siuly, and G. Sinha, _Classification of schizophrenia patients through empirical wavelet transformation using electroencephalogram signals_ (2020).
  * Aouinet and Adnane (2014) A. Aouinet and C. Adnane, Signal Processing: An International Journal (SPIJ) 8, 1 (2014).
  * Ovchinnikov et al. (2011) A. Ovchinnikov, A. Hramov, A. Luttjehann, A. Koronovskii, and G. van Luijtelaar, Technical Physics 56, 1 (2011).
  * Montefusco (2014) L. Montefusco, _Wavelets_ (Elsevier Science, 2014).
  * Maksimenko et al. (2018a) V. A. Maksimenko, A. E. Runnova, M. O. Zhuravlev, P. Protasov, R. Kulanin, M. V. Khramova, A. N. Pisarchik, and A. E. Hramov, PloS one 13, e0197642 (2018a).
  * Maksimenko et al. (2018b) V. A. Maksimenko, A. E. Runnova, N. S. Frolov, V. V. Makarov, V. Nedaivozov, A. A. Koronovskii, A. Pisarchik, and A. E. Hramov, Physical Review E 97, 052405 (2018b).
  * Hramov et al. (2015) A. E. Hramov, A. A. Koronovskii, V. A. Makarov, A. N. Pavlov, and E. Sitnikova, _Wavelets in neuroscience_ (Springer, 2015).
  * Cavalier and O’Hagan (2017) P. Cavalier and D. O’Hagan, in _International Conference on Engineering Geophysics, Al Ain, United Arab Emirates, 9-12 October 2017_ (Society of Exploration Geophysicists, 2017), pp. 128–131.
  * Runnova et al. (2017) A. E. Runnova, M. O. Zhuravlev, A. N. Pysarchik, M. V. Khramova, and V. V. Grubov, in _Dynamics and Fluctuations in Biomedical Photonics XIV_ (International Society for Optics and Photonics, 2017), vol. 10063, p. 1006319.
  * Grubov et al. (2017a) V. Grubov, V. Musatov, V. Maksimenko, A. Pisarchik, A. Runnova, and A. Hramov, Cybernetics and Physics 6, 103 (2017a).
  * Grubov et al. (2017b) V. Grubov, A. Runnova, M. Zhuravlev, V. Maksimenko, S. Pchelintseva, and A. Pisarchik, Cybernetics And Physics (6) pp. 108–113 (2017b).
  * Hramov et al. (2017) A. E. Hramov, V. A. Maksimenko, S. V. Pchelintseva, A. E. Runnova, V. V. Grubov, V. Y. Musatov, M. O. Zhuravlev, A. A. Koronovskii, and A. N. Pisarchik, Frontiers in neuroscience 11, 674 (2017).
  * Chollet et al. (2015) F. Chollet et al., GitHub (2015), URL https://github.com/fchollet/keras.
  * Zhu et al. (2014) G. Zhu, Y. Li, and P. Wen, IEEE Journal of Biomedical and Health Informatics 18, 1813 (2014).
  * Semyachkina-Glushkovskaya et al. (2018) O. Semyachkina-Glushkovskaya, V. Chekhonin, D. Bragin, O. Bragina, E. Vodovozova, A. Alekseeva, V. Salmin, A. Morgun, N. Malinovskaya, E. Osipova, et al., bioRxiv (2018), https://www.biorxiv.org/content/early/2018/12/31/509042.full.pdf, URL https://www.biorxiv.org/content/early/2018/12/31/509042.
  * Van Rijsbergen (1979) C. J. Van Rijsbergen, _Information Retrieval_ (Information Retrieval Group, University of Glasgow, 1979).
  * Fawcett (2006) T. Fawcett, Pattern Recognition Letters 27, 861 (2006).
  * Powers (2011) D. M. W. Powers, Journal of Machine Learning Technologies 2, 37 (2011).
  * Pavlov et al. (2020) A. Pavlov, A. Dubrovsky, A. Koronovskii Jr, O. Pavlova, O. Semyachkina-Glushkovskaya, and J. Kurths, Chaos, Solitons & Fractals 139, 109989 (2020).