inequaity comes from the fact that $\sum_{j\in
Q^{{}^{\prime}}}\tilde{p}(A_{j})\in[\frac{{\gamma^{{}^{\prime}}}^{\frac{1}{2}}\delta^{\frac{1}{4}}}{2},2{\gamma^{{}^{\prime}}}^{\frac{1}{2}}\delta^{\frac{1}{4}}]$.
Here again by choosing $\delta$ small enough, particularly
$\delta\leq\left(\frac{\varepsilon}{2R{\gamma^{{}^{\prime}}}^{\frac{1}{2}}}\right)^{4}$,
we can ensure the change in length is smaller than $\varepsilon$, hence we
achieved our goal.

Proof of Proposition 3.8: Let $h^{*}\in\mathcal{H}$ be an optimal solution of
the Relaxed Minimax Problem and the Relaxed Primary Problem, when the
optimization is over all the measurable functions from $\mathcal{X}$ to
$\mathbb{R}$. Let us denote the corresponding optimal solution to the outer
minimization of the Relaxed Minimax Problem by $f^{*}$. Recall that
$\mathcal{F}$ is a finite dimensional affine class of functions. That is to
say there exists a vector $\boldsymbol{\beta^{*}}$ such that
$f^{*}(x)=\langle\boldsymbol{\beta^{*}},\Phi(x)\rangle\equiv
f_{\boldsymbol{\beta^{*}}}(x)$. Since $h^{*}$ is a solution of the Relaxed
Primary Problem so it should be conditionally valid with respect to class
$\mathcal{F}$. Therefore, as a result of Lemma 3.5, we should have,

$\displaystyle\nabla_{\boldsymbol{\beta}}\;g_{\alpha}(f_{\boldsymbol{\beta}},h^{*})\bigg{|}_{\boldsymbol{\beta^{*}}}=\vec{0}.$
(27)

Now define
$\text{obj}(\boldsymbol{\beta})=\underset{\theta\in\Theta}{\max}\;g_{\alpha}(f_{\boldsymbol{\beta}},h_{\theta})$.
Since $\Theta$ is a compact set and
$g_{\alpha}(f_{\boldsymbol{\beta}},h_{\theta})$ is convex (in fact linear) in
$\boldsymbol{\beta}$ so the Danskin’s theorem applies to
$\text{obj}(\boldsymbol{\beta})$. Therefore, as a result of (27),

$\displaystyle\vec{0}\in\partial_{\boldsymbol{\beta}}\;\text{obj}(\boldsymbol{\beta})\bigg{|}_{\boldsymbol{\beta^{*}}}.$
(28)

Now again as a result of Danskin’s theorem, $\text{obj}(\boldsymbol{\beta})$
is convex in $\boldsymbol{\beta}$. Therefore, as a consequence of (28),
$\boldsymbol{\beta^{*}}$ is a minimizer of $\text{obj}(\boldsymbol{\beta})$.
That is to say, $(f^{*},h^{*})$ is an optimal solution to the Relaxed Minimax
Problem. On the other hand, $h^{*}$ is also a solution to the Relaxed Primary
Problem (as it is a solution to the Relaxed Primary Problem when solved over
all the measurable functions). Putting everything together, $h^{*}$ is a joint
optimal solution to both the Relaxed Minimax Problem and the Relaxed Primary
Problem, hence we have strong duality.

Proof of Proposition 3.4: Let us start by recalling the definitions,

$\displaystyle
g_{\alpha}(f,C)=\operatorname{\mathbb{E}}\left[f(X)\bigg{\\{}{\boldsymbol{1}}[Y\in
C(X)]-(1-\alpha)\bigg{\\}}\right]-\operatorname{\mathbb{E}}\int_{\mathcal{Y}}{\boldsymbol{1}}[y\in
C(X)]dy$

where

$\displaystyle C_{f}(x)=\\{y\in\mathcal{Y}\;|\;f(x)p(y|x)\geq 1\\}.$

Now, fixing $f\in\mathcal{F}$, all we have to show is that
$g_{\alpha}(f,C_{f}(x))-g_{\alpha}(f,C(x))\geq 0$ for every
$C(x):\mathcal{X}\rightarrow 2^{\mathcal{Y}}$.

###### Claim 2

$g_{\alpha}(f,C_{f}(x))-g_{\alpha}(f,C(x))\geq 0$ for every
$C(x):\mathcal{X}\rightarrow 2^{\mathcal{Y}}$

proof: We have,

$\displaystyle g_{\alpha}(f,C_{f}(x))-$ $\displaystyle
g_{\alpha}(f,C(x))\stackrel{{\scriptstyle(a)}}{{=}}\operatorname{\mathbb{E}}_{X}\int_{\mathcal{Y}}\left(f(X)p(y|X)-1\right)\bigg{\\{}\left({\boldsymbol{1}}[y\in
C_{f}(X)]-{\boldsymbol{1}}[y\in C(X)]\right)\bigg{\\}}dy$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{=}}\operatorname{\mathbb{E}}_{X}\int_{\mathcal{Y}}$
$\displaystyle\left(f(X)p(y|X)-1\right)\bigg{\\{}{\boldsymbol{1}}[y\in
C_{f}(X)\backslash C(X)]-{\boldsymbol{1}}[y\in C(X)\backslash
C_{f}(X)]\bigg{\\}}dy$
$\displaystyle\stackrel{{\scriptstyle(c)}}{{=}}\operatorname{\mathbb{E}}_{X}\int_{\mathcal{Y}}$
$\displaystyle\bigg{|}\left(f(X)p(y|X)-1\right)\bigg{|}\bigg{\\{}{\boldsymbol{1}}[y\in
C_{f}(X)\Delta C(X)]\bigg{\\}}dy$ $\displaystyle\geq 0,\hskip 19.91692pt$

where, (a) follows from the definitions, (b) follows from the definition of
set difference operation ($\backslash$) where, $A\backslash B=\\{x\mid x\in A$
and $x\notin B\\}$, and (c) comes from the definition of $C_{f}(x)$. The proof
is complete now. One can define,

$\mathbb{E}_{X}[{\rm
d}(C_{f},C)]\equiv\operatorname{\mathbb{E}}_{X}\int_{\mathcal{Y}}\bigg{|}\left(f(X)p(y|X)-1\right)\bigg{|}\bigg{\\{}{\boldsymbol{1}}[y\in
C_{f}(X)\Delta C(X)]\bigg{\\}}dy,$

and rewrite the above calculation as,

$g_{\alpha}(f,C(x))=g_{\alpha}(f,C_{f}(x))-\mathbb{E}_{X}[{\rm d}(C_{f},C)].$

This reformulation is very intuitive. at a high level, it says maximizing
$g_{\alpha}(f,C(x))$ over $C$, is equivalent to minimizing a distance between
$C$ and $C_{f}$.

Proof of Proposition 3.1: Let us start by restating the proposition.

###### Proposition D.2

The Primary Problem and the Minimax Problem are equivalent. Let
$(f^{*},C^{*}(x))$ be the optimal solution of Minimax Problem. Then, $C^{*}$
is also the optimal solution of the PP. Furthermore, $C^{*}$ has the following
form:

$C^{*}(x)=\\{y\in\mathcal{Y}\;|\;f^{*}(x)p(y|x)\geq 1\\}$ (29)

Let us also recall the definition of Primary Problem,  Primary Problem (PP):
$\displaystyle\underset{C(x)}{\text{Minimize}}$
$\displaystyle\operatorname{\mathbb{E}}\left[\text{len}(C(X))\right]$ subject
to $\displaystyle\mathbb{E}\left[f(X)\bigg{\\{}{\boldsymbol{1}}[Y\in
C(X)]-(1-\alpha)\bigg{\\}}\right]=0,\quad\forall f\in\mathcal{F}$ Recall that
in this paper we assume that
$\mathcal{F}=\\{\langle\beta,\Phi(x)\rangle|\beta\in\mathbb{R}^{d}\\}$ is a
finite dimensional affine class of functions (see section 2.1).

Assuming $\Phi=[\phi_{1},\phi_{2},\cdots,\phi_{d}]$. The Primary Problem can
be equivalently written in terms of $d$ linear constraints on the prediction
sets $C(x)$.  Primary Problem- finite constraints:
$\displaystyle\underset{C(x)}{\text{Minimize}}$
$\displaystyle\operatorname{\mathbb{E}}\left[\text{len}(C(X))\right]$ subject
to $\displaystyle\mathbb{E}\left[\phi_{i}(X)\bigg{\\{}{\boldsymbol{1}}[Y\in
C(X)]-(1-\alpha)\bigg{\\}}\right]=0,\quad\forall i\in[1,d].$ We now take a
closer look at the main optimization variable, i.e. the prediction sets
$C(x)$, and put it in the proper format. The prediction sets $C(x)$ can be
equivalently represented by a function
$C:\mathcal{X}\times\mathcal{Y}\rightarrow\\{0,1\\}$ such that
$C(x,y)={\boldsymbol{1}}[y\in C(x)]$. Furthermore, we can expand the
optimization domain from
$C(x,y):\mathcal{X}\times\mathcal{Y}\rightarrow\\{0,1\\}$ to
$C(x,y):\mathcal{X}\times\mathcal{Y}\rightarrow[0,1]$. One should pay
attention that this change does not affect the optimal solutions as the
solutions will be integer values. Putting everything together we can write the
following problem,

Linear PP
$\displaystyle\underset{C(x,y):\mathcal{X}\times\mathcal{Y}\rightarrow[0,1]}{\text{Minimize}}$
$\displaystyle\int_{\mathcal{X}\times\mathcal{Y}}C(x,y)p(x)\nu(y)dxdy$ subject
to
$\displaystyle\mathbb{E}\left[\phi_{i}(X)\bigg{\\{}C(X,Y)-(1-\alpha)\bigg{\\}}\right]=0,\quad\forall
i\in[1,d].$

Here $\nu$ is the lebesgue measure; recall from Section 2.1 that we defined
length through the lebesgue measure on $\mathcal{Y}=\mathbb{R}$, which is
equivalent to cardinality of a set in the case where $\mathcal{Y}$ is a
discrete and finite set. This problem is a Linear program in terms of $C$ with
finitely many constraints. Also, the Linear PP, includes the Primary Problem
in the sense that any solution to the Primary Problem is a feasible solution
for Linear PP. That is to say, to prove Proposition 3.1, we just have to prove
it for Linear PP.

For $d=1$, Linear PP can be seen through the lens of the Neyman-Pearson Lemma
in Hypothesis Testing [68]. As a result, our analysis of the Linear PP (for
general $d$) can be considered as an extension of this lemma which will be
done using the KKT conditions.

Now, let us rewrite the Linear PP:

Minimize
$\displaystyle\int_{\mathcal{X}\times\mathcal{Y}}C(x,y)p(x)\nu(y)dxdy$ subject
to:
$\displaystyle\int_{\mathcal{X}\times\mathcal{Y}}\phi_{i}(x)C(x,y)p(x,y)dxdy-(1-\alpha)=0,\quad\forall
i\in[1,d]$ $\displaystyle C(x,y)\in[0,1]\quad\forall
x\in\mathcal{X},y\in\mathcal{Y}$

The above is a standard linear program on $C(x,y)$. In what follows, we will
find a closed-form solution using the dual of this program. Here, the
“optimization variable” $C(x,y)$ belongs to an infinite-dimensional space.
Hence, in order to be fully rigorous, we will need to use the duality theory
developed for general linear spaces that are not necessarily finite-
dimensional. For a reader who is less familiar with infinite-dimensional
spaces, what appears below is a direct extension of the duality theory (i.e.
writing the Lagrangian) for the usual linear programs in finite-dimensional
spaces.

Let $\mathcal{F}$ be the set of all measurable function defined on
$\mathcal{X}\times\mathcal{Y}$. Note that $\mathcal{F}$ is a linear space. Let
$\Omega$ be the set of all the measurable functions on
$\mathcal{X}\times\mathcal{Y}$ which are bounded between $0$ and $1$; I.e.

$\Omega=\left\\{C\in\mathcal{F}\text{ s.t.
}C:\mathcal{X}\times\mathcal{Y}\to[0,1]\right\\}$ (30)

Note that $\Omega$ is a convex set. We can then rewrite our linear program as
follows:

Minimize
$\displaystyle\int_{\mathcal{X}\times\mathcal{Y}}C(x,y)p(x)\nu(y)dxdy$ subject
to:
$\displaystyle\int_{\mathcal{X}\times\mathcal{Y}}\phi_{i}(x)C(x,y)p(x,y)dxdy-(1-\alpha)=0,\quad\forall
i\in[1,d]$ $\displaystyle C\in\Omega$

Moreover, let us define the functional $F:\mathcal{F}\to\mathbb{R}$ as

$F(C)=\int_{\mathcal{X}\times\mathcal{Y}}C(x,y)p(x)\nu(y)dxdy,$ (31)

and also define, for $i\in[1,d]$, the functional
$G_{i}:\mathcal{F}\to\mathbb{R}$ as

$G_{i}(C)=\int_{\mathcal{X}\times\mathcal{Y}}\phi_{i}(x)C(x,y)p(x,y)dxdy-(1-\alpha).$
(32)

Finally, we define the mapping ${\boldsymbol{G}}:\mathcal{F}\to\mathbb{R}^{d}$
as

${\boldsymbol{G}}(C)=\left[G_{1}(C),G_{2}(C),\cdots,G_{d}(C)\right].$

Note that $G$ is a linear (and hence convex) mapping from $\mathcal{F}$ to the
Euclidean space $\mathbb{R}^{d}$.

Using the above-defined notation, our linear program becomes:

Minimize $\displaystyle F(C)$ subject to:
$\displaystyle{\boldsymbol{G}}(C)=\mathbf{0}$ $\displaystyle C\in\Omega$

where $\mathbf{0}\in\mathbb{R}^{d}$ is the all-zero vector.

Note that the feasibility set of the above program is non-empty, as
$C(x,y)=1-\alpha$, for all $(x,y)\in\mathcal{X}\times\mathcal{Y}$, is a
feasible point. We can now use the duality theory of convex programs in vector
spaces (See Theorem 1, Section 8.3 of [67]; Also see Problem 7 in Chapter 8 of
the same reference). Specifically, let OPT be the optimal value achievable in
the above linear program. Then, there exists a vector
$\mathbf{{\boldsymbol{\beta}}}\in\mathbb{R}^{d}$ such that the following
holds:

$\text{OPT}=\inf_{C\in\Omega}\left\\{F(C)-\langle{\boldsymbol{\beta}},{\boldsymbol{G}}(C)\rangle\right\\},$
(33)

where $\langle{\boldsymbol{\beta}},{\boldsymbol{G}}(C)\rangle$ denotes the
Euclidean inner-product of the two vectors
${\boldsymbol{\beta}},{\boldsymbol{G}}(C)\in\mathbb{R}^{d}$. Here, note that
the vector ${\boldsymbol{\beta}}$ is the usual Lagrange multiplier.

By denoting
${\boldsymbol{\Phi}}(x)=(\phi_{1}(x),\phi_{2}(x),\cdots,\phi_{d}(x))$, and
using (32), we can write

$\langle{\boldsymbol{\beta}},{\boldsymbol{G}}(C)\rangle=\int_{\mathcal{X}\times\mathcal{Y}}\langle{\boldsymbol{\beta}},{\boldsymbol{\Phi}}(x)\rangle
C(x,y)p(x,y)dxdy-(1-\alpha)\langle{\boldsymbol{\beta}},\mathbf{1}\rangle,$

where $\mathbf{1}\in\mathbb{R}^{d}$ is the all-ones vector. As a result, by
using (31), in order to solve the optimization in (33) we need to solve the
following optimization:

$\displaystyle\inf_{C\in\Omega}\left\\{\int_{\mathcal{X}\times\mathcal{Y}}C(x,y)\left(p(x)\nu(y)-\langle{\boldsymbol{\beta}},{\boldsymbol{\Phi}}(x)\rangle
p(x,y)\right)dxdy\right\\}+(1-\alpha)\langle{\boldsymbol{\beta}},\mathbf{1}\rangle.$

And by noting the fact the $p(x,y)=p(x)p(y\mid x)$, and removing the term
$(1-\alpha)\langle{\boldsymbol{\beta}},\mathbf{1}\rangle$ which is independent
of $C$, our dual optimization becomes:

$\displaystyle\inf_{C\in\Omega}\left\\{\int_{\mathcal{X}\times\mathcal{Y}}C(x,y)\biggl{(}\nu(y)-\langle{\boldsymbol{\beta}},{\boldsymbol{\Phi}}(x)\rangle
p(y|x)\biggr{)}p(x)dxdy\right\\}$ (34)

Now, it is easy to see that as $C(x,y)\in[0,1]$ (due to the constraint
$C\in\Omega$), the above optimization problem has a closed-form solution
$C^{*}$:

$C^{*}(x,y)=\begin{cases}1&\quad\text{if
}\langle{\boldsymbol{\beta}},{\boldsymbol{\Phi}}(x)\rangle p(y|x)>\nu(y),\\\
0&\quad\text{if }\langle{\boldsymbol{\beta}},{\boldsymbol{\Phi}}(x)\rangle
p(y|x)<\nu(y),\\\ \in\\{0,1\\}&\quad\text{otherwise}.\end{cases}$ (35)

Finally, note that in this paper the measure $\nu$ is considered to be the
Lebesgue measure – see Section 2.1 – hence, up to a constant normalization
factor that can be absorbed into ${\boldsymbol{\beta}}$, we have $\nu(y)=1$
for all $y\in\mathcal{Y}$.

The structure given in (35) is also sufficient. I.e. for any pair
$({\boldsymbol{\beta}},C^{*})$ such that (i) $C^{*}$ has the form given in
(35); and (2) $C^{*}$ satisfies the coverage-constraints
${\boldsymbol{G}}(C^{*})=0$, then $C^{*}$ is an optimal solution of the
Primary Problem(LABEL:primary_proper). This sufficiency result follows again
from the duality theory of linear programs in linear spaces; E.g. see Theorem
1 in Section 8.4 of [67]. This necessary and sufficient condition is known as
Strong Duality, which then results in the equivalence of Minimax Problem and
Primary Problem, i.e we can change the order of min and max.

##  Appendix E Proofs of Section 4

Proof of Theorem 4.4: Let us recall the algorithm.  CPL:
$\displaystyle\underset{f\in\mathcal{F}}{\text{Minimize}}\;\underset{h\in\mathcal{H}}{\text{Maximize}}\;\tilde{g}_{\alpha,n}(f,h),$
$\displaystyle\text{where,}\;C_{h}^{S}(x)=\\{y\in\mathcal{Y}\;|\;S(x,y)\leq
h(x)\\}.$ For the ease of notation let us call
$\tilde{g}(f,h)=\tilde{g}_{\alpha,n}(f,h)$. let us also call the stationary
solution to CPL by $h_{\rm CPL}^{*}$ and $f_{\rm CPL}^{*}$. Now since
$\tilde{g}$ is a smooth function with respect to both of its arguments we have
the following optimality condition,

$\displaystyle\frac{d}{d\varepsilon}\tilde{g}(f_{\rm CPL}^{*}+\varepsilon
f,h_{\rm CPL}^{*})\bigg{|}_{\varepsilon=0}=0,\quad\text{for
every}\;f\in\mathcal{F}.$ (36)

Recall that $\mathcal{F}$ is a $d$-dimensional affine class of functions over
the basis $\Phi=[\phi_{1},\cdots,\phi_{d}]$. We can then rewrite 36 as what
follows.

$\displaystyle\frac{d}{d\varepsilon}\tilde{g}(f_{\rm
CPL}^{*}+\varepsilon\phi_{j},h_{\rm
CPL}^{*})\bigg{|}_{\varepsilon=0}=0,\quad\text{for every}\;j\in[1,\cdots,d].$
(37)

Now looking at the definition of $g(f,h)=$

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\left[f(x_{i})\bigg{\\{}\tilde{{\boldsymbol{1}}}[S(x_{i},y_{i}),h(x_{i})]-(1-\alpha)\bigg{\\}}\right]-\frac{1}{n}\sum_{i=1}^{n}\int_{\mathcal{Y}}\bigg{(}\tilde{{\boldsymbol{1}}}[S(x_{i},y)\leq
h(x_{i})]-(1-\alpha)\bigg{)}dy,$

We can take the derivative with respect to $f$. Hence we can rewrite 37 as
what follows:

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\left[\phi_{j}(x_{i})\bigg{\\{}\tilde{{\boldsymbol{1}}}[S(x_{i},y_{i}),h_{\rm
CPL}^{*}(x_{i})]-(1-\alpha)\bigg{\\}}\right]=0,\quad\text{for
every}\;j\in[1,\cdots,d].$ (38)

One can think of the mathematical term above as the smoothed version of
coverage under covariate shift $\phi_{j}$. Now we can apply Lemma C.3.
Therefore, fixing $j\in[1,\cdots,d]$, with probability $1-\delta$ we have,

$\displaystyle\bigg{|}\frac{1}{n}\sum_{i=1}^{n}\bigg{[}\phi_{j}(x_{i})\bigg{\\{}\tilde{{\boldsymbol{1}}}[S(x_{i},y_{i}),h^{*}(x_{i})]-(1-\alpha)\bigg{\\}}\bigg{]}$
$\displaystyle-\operatorname{\mathbb{E}}\left[\phi_{j}(x)\bigg{\\{}\tilde{{\boldsymbol{1}}}[S(x,y),h_{\rm
CPL}^{*}(x)]-(1-\alpha)\bigg{\\}}\right]\bigg{|}$
$\displaystyle\leq\frac{2B\varepsilon}{\sqrt{2\pi}\sigma}+\frac{\sqrt{2}B\sqrt{\ln\left(\frac{2\mathcal{N}(\mathcal{H},d_{\infty},\varepsilon)}{\delta}\right)}}{\sqrt{n}}\quad\text{for
any}\;\varepsilon>0.$

Combining with (38) we get,

$\displaystyle\bigg{|}\operatorname{\mathbb{E}}\left[\phi_{j}(x)\bigg{\\{}\tilde{{\boldsymbol{1}}}[S(x,y),h_{\rm
CPL}^{*}(x)]-(1-\alpha)\bigg{\\}}\right]\bigg{|}\leq\frac{2B\varepsilon}{\sqrt{2\pi}\sigma}+\frac{\sqrt{2}B\sqrt{\ln\left(\frac{2\mathcal{N}(\mathcal{H},d_{\infty},\varepsilon)}{\delta}\right)}}{\sqrt{n}}\quad\text{for
any}\;\varepsilon>0.$

Union bounding over (38) we have with probability $1-\delta$ for any
$j\in[1,\cdots,d]$,

$\displaystyle\bigg{|}\operatorname{\mathbb{E}}\left[\phi_{j}(x)\bigg{\\{}\tilde{{\boldsymbol{1}}}[S(x,y),h_{\rm
CPL}^{*}(x)]-(1-\alpha)\bigg{\\}}\right]\bigg{|}\leq\frac{2B\varepsilon}{\sqrt{2\pi}\sigma}+\frac{\sqrt{2}B\sqrt{\ln\left(\frac{2d\mathcal{N}(\mathcal{H},d_{\infty},\varepsilon)}{\delta}\right)}}{\sqrt{n}}\quad\text{for
any}\;\varepsilon>0.$ (39)

In other words, we proved that the expected smoothed version of coverage is
bounded. The last step that we have to take is then to prove that the expected
smoothed coverage is actually close to the expected actual coverage. The
following claim makes this precise.

###### Claim 3

The following inequality holds for every $j\in[1,\cdots,d]$.

$\displaystyle\bigg{|}\operatorname{\mathbb{E}}\left[\phi_{j}(x)\bigg{\\{}\tilde{{\boldsymbol{1}}}[S(x,y),h_{\rm
CPL}^{*}(x)]-(1-\alpha)\bigg{\\}}\right]$
$\displaystyle-\operatorname{\mathbb{E}}\left[\phi_{j}(x)\bigg{\\{}{\boldsymbol{1}}[S(x,y)\leq
h_{\rm CPL}^{*}(x)]-(1-\alpha)\bigg{\\}}\right]\bigg{|}$ $\displaystyle\leq
BL\sigma\sqrt{\frac{2}{\pi}}$

Proof.

$\displaystyle\bigg{|}\operatorname{\mathbb{E}}\left[\phi_{j}(x)\bigg{\\{}\tilde{{\boldsymbol{1}}}[S(x,y),h_{\rm
CPL}^{*}(x)]-(1-\alpha)\bigg{\\}}\right]-\operatorname{\mathbb{E}}\left[\phi_{j}(x)\bigg{\\{}{\boldsymbol{1}}[S(x,y)\leq
h_{\rm CPL}^{*}(x)]-(1-\alpha)\bigg{\\}}\right]\bigg{|}$
$\displaystyle=\bigg{|}\operatorname{\mathbb{E}}\left[\phi_{j}(x)\tilde{{\boldsymbol{1}}}[S(x,y),h_{\rm
CPL}^{*}(x)]\right]-\operatorname{\mathbb{E}}\left[\phi_{j}(x){\boldsymbol{1}}[S(x,y)\leq
h_{\rm CPL}^{*}(x)]\right]\bigg{|}$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}\operatorname{\mathbb{E}}\left[\bigg{|}\phi_{j}(x)\bigg{|}\bigg{|}\left(\tilde{{\boldsymbol{1}}}(S(x,y),h_{\rm
CPL}^{*}(x))-{\boldsymbol{1}}[S(x,y)\leq h_{\rm
CPL}^{*}(x)]\right)\bigg{|}\right]$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{\leq}}B\operatorname{\mathbb{E}}\left[\bigg{|}\left(\tilde{{\boldsymbol{1}}}(S(x,y),h_{\rm
CPL}^{*}(x))-{\boldsymbol{1}}[S(x,y)\leq h_{\rm
CPL}^{*}(x)]\right)\bigg{|}\right]$
$\displaystyle=B\operatorname{\mathbb{E}}_{X}\operatorname{\mathbb{E}}_{S|X}\left[\bigg{|}\left(\tilde{{\boldsymbol{1}}}(S(x,y),h_{\rm
CPL}^{*}(x))-{\boldsymbol{1}}[S(x,y)\leq h_{\rm
CPL}^{*}(x)]\right)\bigg{|}\right]$
$\displaystyle\stackrel{{\scriptstyle(c)}}{{\leq}}B\operatorname{\mathbb{E}}_{X}L\sigma\sqrt{\frac{2}{\pi}}$
$\displaystyle=BL\sigma\sqrt{\frac{\pi}{2}},$

where (a) comes from triangle inequality, (b) is derived by the definition of
$B=\max_{i\in[1,\cdots,d]}\sup_{x\in\mathcal{X}}\Phi_{i}(x)$, and (c) is
followed by assumption 2 and Lemma C.2. Now combining Claim 3 and (39) we have
with probability $1-\delta$,

$\displaystyle\bigg{|}\operatorname{\mathbb{E}}\left[\phi_{j}(x)\bigg{\\{}{\boldsymbol{1}}[S(x,y)\leq
h_{\rm CPL}^{*}(x)]-(1-\alpha)\bigg{\\}}\right]\bigg{|}\leq$ $\displaystyle
BL\sigma\sqrt{\frac{\pi}{2}}+\frac{2B\varepsilon}{\sqrt{2\pi}\sigma}+\frac{\sqrt{2}B\sqrt{\ln\left(\frac{2d\mathcal{N}(\mathcal{H},d_{\infty},\varepsilon)}{\delta}\right)}}{\sqrt{n}}\quad\text{for
any}\;\varepsilon>0.$

Putting $\sigma=\frac{1}{\sqrt{n}}$ and $\varepsilon=\frac{1}{n}$ we have for
every $j\in[1,\cdots,d]$,

$\displaystyle\bigg{|}\operatorname{\mathbb{E}}\left[\phi_{j}(x)\bigg{\\{}{\boldsymbol{1}}[S(x,y)\leq
h_{\rm
CPL}^{*}(x)]-(1-\alpha)\bigg{\\}}\right]\bigg{|}\leq\frac{\sqrt{2}B\sqrt{\ln\left(\frac{2d\mathcal{N}(\mathcal{H},d_{\infty},\frac{1}{n})}{\delta}\right)}+BL\sqrt{\frac{\pi}{2}}+\frac{2B}{\sqrt{2\pi}}}{\sqrt{n}}$

Let us also remind that each element $f\in\mathcal{F}$ can be represented by a
$\boldsymbol{\beta}\in\mathbb{R}^{d}$, where we use the notation
$f(x)=\langle{\boldsymbol{\beta}},\Phi(x)\rangle\equiv
f_{{\boldsymbol{\beta}}}(x)$ (look at section 2.1 for more details). By
linearity of the class $\mathcal{F}$ we can conclude with probability
$1-\delta$ for every $f_{\boldsymbol{\beta}}\in\mathcal{F}$,

$\displaystyle\bigg{|}\operatorname{\mathbb{E}}\left[f(x)\bigg{\\{}{\boldsymbol{1}}[S(x,y)\leq
h_{\rm
CPL}^{*}(x)]-(1-\alpha)\bigg{\\}}\right]\bigg{|}\leq\frac{||{\boldsymbol{\beta}}||_{1}\sqrt{2}B\sqrt{\ln\left(\frac{2d\mathcal{N}(\mathcal{H},d_{\infty},\frac{1}{n})}{\delta}\right)}+||{\boldsymbol{\beta}}||_{1}BL\sqrt{\frac{\pi}{2}}+||{\boldsymbol{\beta}}||_{1}\frac{2B}{\sqrt{2\pi}}}{\sqrt{n}}.$

##  Appendix F Refrences for 11 datasets for Section 5.1

Here we list all the datasets.

  * •

MEPS-19 [69]

  * •

MEPS-20 [70]

  * •

MEPS-21 [70]

  * •

blog feedback (blog-data)[71]

  * •

physicochemical properties of protein tertiary structure (bio) [72]

  * •

bike sharing (bike) [73]

  * •

community and crimes (community) [74]

  * •

Tennessee’s student teacher achievement ratio (STAR) [75]

  * •

concrete compressive strength (concrete) [76]

  * •

Facebook comment volume variants one (facebook-1) [77]

  * •

Facebook comment volume variants two (facebook-2) [78]


# Quantifying the Effects of the 2008 Recession using the Zillow Dataset

Arunav Gupta  Lucas Nguyen  Camille Dunning  Ka Ming Chan

( Halıcıoğlu Data Science Institute
December 2019 )

##  1 Background of Crisis

The 2008 Recession, also known as the “subprime mortgage crisis” was an
economic recession that started in the US and quickly had global economic
implications. From a period from late 2007 to mid-2009, 8.4 million Americans
lost their jobs, 1 in every 54 homes filed for foreclosure, and the US GDP
fell by 4.3 percent.

Most economists agree that the crisis has its roots in the overuse of subprime
mortgages, which were a type of security usually given to those who exhibit
high financial risk and have low credit scores. Starting in the early 2000s,
the housing industry was booming, and as houses steadily increased in value,
more banks felt they could assume the risk of subprime mortgages and thus
started giving them out en masse, even to those who would have traditionally
been denied a loan. However, in 2007, the housing bubble burst and home values
fell by almost 31.8%. As demand for housing plateaued, home values dropped,
and people who had taken out subprime loans found that they were no longer
able to pay off the high interest payments associated with a subprime
mortgage. As thousands defaulted on their mortgages, the lending institutions
lost money, too. Compounding on the current problems, many of these
institutions had traded mortgage-backed securities (MBSs), which were backed
by these risky subprime loans, to other institutions seeking a profit when
housing prices increased. When that didn’t happen, MBSs lost value as well,
causing many banks like the Lehman Brothers to go bankrupt. The collapse of
the real estate industry then caused banks and businesses to lose trust in
each other, driving stock prices down. As publicly-traded companies then saw
decreased valuations, businesses shut down, and unemployment skyrocketed.

This “perfect storm” of financial disasters – subprime mortgages, housing
value decline, mistrust in the stock market – hurt those at the bottom of the
economic ladder the most, as without a job or a good credit rating to fall
back on, they were unable to buy a home or provide for their families.

###  1.1 Maps

To visualise the recession effect on home values, we plotted the average ZHVI
for each state per year, as well as the differences between each year.

(a) Figure A

(b) Figure B

Figure 1: ZHVI, 2007-2015

The first map is normalized to a scale of 0 to 550,000 second map is
normalized to a scale of -50,000 to +150,000. Overall we see a significant
drop and then recovery reflected in the maps.

##  2 Description of Dataset

Our focus dataset was the ‘Zillow Economics Data’ found on Kaggle, which
includes records of transactions, contracts, and specs about public properties
throughout neighbourhoods, cities, statistical metropolitan areas, counties
and states. The data are in a time series format starting from 1996 to around
2014.

The observations in this dataset are derived from property listings and user
behaviour on Zillow, and statistics were computed for metrics including the
sales price, listing price, rent price. These statistics were acquired for
properties per number of bedrooms, property type, property tier, and overall.
Some other features include, but aren’t limited to, the number of days the
property was listed on Zillow, the raw inventory, the percentage of annual
increase or decrease in property value, sales turnover, and raw sales.

We decided to conduct our analyses on the property statistics across
metropolitan areas. In our data, metropolitan areas are specified by a CBSA
code. This code corresponds to a certain, CBSA, or core-based statistical
area, which is characterised by one or more counties that are anchored by an
urban centre. The Zillow Economics Data include both micropolitan and
metropolitan CBSAs, where the urban centre of a micropolitan CBSA is between
10,000 and 50,000 and the metropolitan urban centre has greater than 50,000
people.

For reasons described in our section on the method of PCA, we decided not to
use statistics for individual metrics as our dependent variables, but instead
another provided measure, the Zillow Home Value Index (ZHVI). The ZHVI is
based on the Zestimate home valuation model; it takes the median estimate in a
geographic area on a given day. The median Zestimate is more sensible to use
as it handles extreme values much better. Thus, we proceeded to use the ZHVI
as the basis for our AUB and ARIMA implementations. There is an equivalent for
rental spaces, the Zillow Rent Index (ZRI), but we didn’t run analyses for
these properties.

##  3 Methodologies

###  3.1 PCA

The first methodology our group has performed for trial is called the
Principal Component Analysis (PCA). The original dataset has approximately 95
columns of variables for us to build our analysis on. For our study, we would
only like to extract a few useful features in the purpose of reducing the
amount of excessive information. One way to reduce our high dimensional 95
columns of data into low dimensional representations of it is by using the
method of PCA. In summary, the PCA selects a few important features from a
sparse vector of data and compresses it by ignoring components which are not
meaningful. Therefore, the data can be recovered and summarized as few
dimensions as possible. The process can simply be described by figure 3 below:

Figure 2: Visualization of PCA

As shown in the process, the way we select the important features is by
keeping the dimensions with the highest variance, and discarding the
dimensions with the lowest variance. The highest variance dimensions maximize
the amount of “randomness” that gets preserved in the compressed data. The
method for determining the compressed data is simply by minimizing the Mean
Square Error (MSE) between the original data and the new compressed line.

We apply the PCA method onto our dataset and have reduced the 95 variables
into 2. Here is a plot of the result in figure 4:

Figure 3: Visualization of PCA on US nationwide data (n=3)

Unfortunately this graph does not provide much unique insights to us, as there
is not much useful information from it that we can extract. The graph is
messy, and not scalable because in the original dataset, most of the regions
do not have complete data in them. In conclusion, the first methodology our
group has tried, the PCA, does not reinforce us in finding any interesting
correlation between our data and the recession.

###  3.2 Area Under Baseline (AUB)

The Area Under Baseline metric seeks to answer the question: “How much did the
recession affect the ZHVI of a city?” To answer this question, we must look at
the total impact the recession had on the city in question. There are multiple
parts to the procedure: (a) transform the ZHVI data into a moving average
(hereafter noted as $ZHVI_{MA}$), (b) find the “recession window” for the
city, and (c) find the area between the $ZHVI_{MA}$ trend and the baseline
ZHVI across the recession window.

Below, we’ll illustrate the procedure on Aberdeen, WA, which has a ZHVI trend
graph that is very conducive to the process. Figure 5 shows Aberdeen’s
$ZHVI_{MA}$, computed with a window = 5. The red lines denote the recession
window: the beginning of the window is the greatest local maximum of the
$ZHVI_{MA}$ after Jan. 1, 2007. The end of the window is the point in time
where the $ZHVI_{MA}$ intersects the value of $ZHVI_{MA}$ at the start point.
The intuition here is that the recession is considered “recovered” once the
ZHVI has returned to its pre-recession value. If the ZHVI never reaches its
pre-recession value, the end date of the recession window is set at the last
available value in the dataset. Figure 5 also shows the baseline (green line),
which is defined as the $ZHVI_{MA}$ value at the start point of the recession
window. Finally, we can compute the residual between the baseline and the
$ZHVI_{MA}$ for each point in time between the start and end date of the
recession window, then sum those values up to get the AUB for Aberdeen, WA.

Figure 4: $ZHVI_{MA}$ for Aberdeen, WA

####  3.2.1 Theory and Notation

First we find the 5 month moving averages of the ZHVI values in the recession
period The 5 month Moving Average for month $iM_{i}$ is:
$M_{i}=\frac{1}{5}\sum_{k=0}^{4}y_{i-k}$

$Y_{i}=$ ZHVI value at month i

We look for local maxes by finding dates where the moving average was positive
before but negative after. If such local max exists, we find the local max
with the largest ZHVI value:

$C=\left\\{y_{i}:M_{i+1}<0,M_{i-1}>0\right\\}$
$c=\max(C)$

The time that the local max occurs is declared the time window start.

$a=i$

If it exists, we declare the time window end b to be the date of the first
ZHVI value greater than the baseline. If not, we choose the date at the end of
the recession period

$b={\operatorname*{arg\,min}}(c,y_{z})$
$y_{z}=$ ZHVI at end of recession period

Finally we find the area under the baseline by taking the sum of the
differences between each ZHVI value in the time window and the baseline.

$\sum_{i=a}^{b}y_{i}-c$
$a=$ window start
$b=$ window end
$Y_{i}=$ ith ZHVI value in window
$c=$ baseline

####  3.2.2 Results of Analysis

Overall, the top 3 and bottom 3 AUB scores and their respective cities are
listed in the tables below. Since a higher AUB score should be interpreted as
an indicator of a higher recession impact, we will call the top 3 and bottom 3
cities “losers” and “gainers,” respectively:

Top ”Loser” Cities | AUB | Top ”Gainer” Cities | AUB
---|---|---|---
Key West, FL | $1.638\times 10^{11}$ | Mt. Vernon, IL | 17,080
Salinas, CA | $1.541\times 10^{11}$ | McAlester, OK | 132,480
Carson City, NV | $1.331\times 10^{11}$ | Norfolk, NE | 287,040
Table 1: Results of AUB analysis

###  3.3 ARIMA Model

####  3.3.1 What is Time Series Forecasting?

There are many problems in predictive modeling that involve a time component.
When we are making predictions about the outcome in the future, we are still
treating all prior observations equally. In time series analysis, we have two
different goals of either trying to understand and describe our time series
data, or making predictions, or forecasting. Descriptive time series analysis
can help with prediction, as it comprehends models to aid in identifying
underlying causes, but it is not required and can be an investment.
Forecasting calls on models to fit historical data and using that information
to predict future observations. ”The purpose of time series analysis is
generally twofold: to understand or model the stochastic mechanisms that gives
rise to an observed series and to predict or forecast the future values of a
series based on the history of that series.”

To better understand time series analysis, we could decompose a time series
into the following parts:

  1. 1.

Level \- The average value in a series.

  2. 2.

Trend \- The often linear increasing or decreasing behaviour of the series
over time. - Optional, contingent on non-stationary or stationary time series

  3. 3.

Seasonality \- Repeating patterns of cycles in behaviour over time. The ZHVI
measure provided by the Kaggle data is already smoothed and seasonally
adjusted.

  4. 4.

Noise \- Variability in observations, unexplainable by model

We can combine these components to provide an observed time series, and add
them together to form our model:

$y(t)=level_{t}+trend_{t}+seasonality_{t}+noise_{t}$

Time series data can also require ample scaling and cleaning to adjust for
uneven frequency, time spacing, outliers, missing values, etc. The ZHVI data
have been cleaned to supply the metro_data.csv table.

The autoregressive integrated moving average (ARIMA) is a time-series fitted
model designed to aid in descriptive analysis and forecasting of time-series
data. ARIMA is often applied to data that show non-stationarity. As implied by
the name, ARIMA has these key attributes:

####  3.3.2 Autoregression

ARIMA employs a simple autoregression (AR) model, in which observations from
previous time steps are used as input to a regression equation to predict the
next value. Formally, we can indicate an autoregressive model of order $p$ by
the following:

$AR(p)=X_{t}=c+\sum_{i=1}^{p}\varphi_{i}X_{t}-i+\varepsilon_{t}$,

where $c$ is a constant, $\varphi_{1},...,\varphi_{p}$ are the parameters of
the model, and $\varepsilon_{t}$ is noise. $p$, the order of the
autoregressive model, represents the number of lags, or previous observed
series values, to be included in the model.

Since we are using regression with a neighborhood of terms, we can express
this model equivalently with a backshift operator, $B$:

$AR(p)=X_{t}=c+\sum_{i=1}^{p}\varphi_{i}B^{i}X_{t}+\varepsilon_{t}$

A backshift, or lag, operator operates on an element of a time series to
produce the previous element. Let us define an arbitrary time series
$X={X_{1},X_{2},...}$. Then $BX_{t}=X_{t-1},\forall t>1$. The backshift
operator can be raised to arbitrary integer powers so that
$B^{k}X_{t}=X_{t-k}$.

####  3.3.3 Integration

The ARIMA model is integrated, meaning it uses a process known as
differencing, in which observations at consecutive time steps are subtracted.
This makes our non-stationary time series stationary, stabilising the mean by
reducing trend. In the cases where seasonality will be reduced, the time
series variance will also be stabilised.

####  3.3.4 Moving Average

We express a moving average model of order $q$:

$MA(q)=X_{t}=\mu+\varepsilon_{t}+\theta_{1}\varepsilon_{t-1}+...+\theta_{q}\varepsilon_{t-q}=\mu+\varepsilon_{t}+\sum_{i=1}^{q}\theta_{i}\varepsilon_{t-i}$,

where $\mu$ is the mean of the series, $\theta_{1},...\theta_{q}$ are the
parameters of the model, and $\varepsilon_{t-1},...,\varepsilon_{t-q}$ are
noise error terms. Write in terms of the backshift operator:

$MA(q)=X_{t}=\mu+(1+\theta_{1}B+...+\theta_{q}B^{q})\varepsilon_{t}$

Simply put, the average is represented here is represented as the central
value of our set of numbers, but it’s calculated for values of the dependent
variable at different time intervals. The order $q$ denotes the size of the
moving average window.

The ARIMA() model in Python accepts $p$, $d$, and $q$, where $d$ is the order
of differencing.

After importing all necessary libraries and table, we proceeded to the ARIMA
analysis on a sample metropolitan area, San Diego. First, we produced visualis
of the ZHVI trend and a correlogram. The correlogram, or autocorrelation plot,
plots the sample correlations of the regression for each lag value. We want to
choose a nonzero value of $p$ for our model such that the autocorrelation is
high, so we can avoid overestimation or underestimation of true values for
training of our forecasting model.

As soon as we are done fitting the model, we will have a summary of the fit.
We have also plotted the distributions for the residual errors, from which we
can maybe capture some trend information. The density plot of the residual
values show that they are Gaussian but not centred at zero. This is indicative
of a bias in the prediction.

##### We now test our model and produce a 95% prediction interval for our
forecasted results for the ZHVI in San Diego from 2017 to 2020. Note that
there is a slight overlap with the in-sample and out-sample predictions.

##### Finally, we calculate the area of the confidence interval.

##### The smaller the area of the 95% confidence interval, the less volatile
the recovery over a longer period of time, and the more certain we are that
there will be a continuing increasing trend in the ZHVI. We proceeded to
calculate the area of the 95% prediction interval for each metropolitan area.
The above shows the distribution of the normalised areas for across our
observations.

##  4 Results

###  4.1 Geographic Clustering

By plotting the cities we received from each of the two working methodologies
in the previous section (AUB and ARIMA), we can observe a couple patterns:

Figure 5: Visualization of results from AUB and ARIMA

First, we can notice that the gainer cities from both methods are primarily
clustered in middle America/Great Plains, while many loser cities are located
in the West. This defies conventional wisdom, which dictates that the West was
not hit as hard as the Great Plains region due to its burgeoning technology
industry. However, we can understand why Key West was marked as a loser, as it
has a big tourism industry that was severely hit during the Recession, when
less people could afford vacations and tourism.

###  4.2 Metrics vs. Population

We plot our AUB and ARIMA outputs against population for correlation. Figure 7
shows the two resulting graphs.

(a) Figure A

(b) Figure B

Figure 6: Regressing both methods against population statistics

Each blue dot represents a metro area. Since the r-value for both graphs are
low, we conclude that there is little to no linear correlation between
population and our algorithms. In the future we may consider transforming our
data before testing for correlation.

###  4.3 Metrics vs. Unemployment

We have also plotted our AUB and ARIMA outputs against unemployment rates in
hopes of finding any possible correlations. Figure 8 shows the two resulting
graphs.

(a) Figure A

(b) Figure B

Figure 7: my caption

Each blue dot in the graphs represents a metro area. Judging from the
regression lines of these two plots, there seems to be a correlation which
reveals that the hardest-hit cities also had the biggest job losses. But as
the R-squared values of them appear to be quite low, the correlations are not
much secured. Therefore in conclusion, our attempts of plotting the population
and the unemployment rates against our AUB and ARIMA outputs unfortunately do
not show significant results.

##  5 Conclusion

Although we were able to come up with two unique metrics from determining the
impact of the recession, neither of them seemed to agree with the universally-
accepted metrics for determining recession impact – unemployment rate and city
size. Our reasons for reaching this conclusion are because the AUB and ARIMA
scores do not correlate very well with population and unemployment rate
intuitively. This could be for a variety of reasons. The most prevalent is the
notion that perhaps housing data such as the ZHVI is not the best indicator of
a recession, after all. Another issue could be noise in the data. In many
cases, the average ZHVI of a given metropolis was much larger than that of any
other one, a feature of the dataset we could have corrected for by normalizing
each metro’s data before applying a moving average.

###  5.1 Future Steps

Moving forward, some steps we may take are as follows:

  1. 1.

An investigation of the spillover effect in our data.

    1. (a)

According to a CityLab article, the spillover effect occurs when economic
event that occurs in one city or metropolitan area will again occur in an
adjacent city or metropolitan area. The study in the article pointed out that
this was evident in a number of metros Chicago, New York, and Hartford, while
Washington, D.C., Austin, and Providence were able to ride out the recession
on their own. Essentially, we will want to check if our data exhibit possible
relationships between ZHVI trends in one city and those of its neighbouring
cities.

  2. 2.

A new algorithm that combines the AUB method and ARIMA forecasting. This may
involve tuning hyperparameters of our ARIMA model to get more detailed
predictive curves/trends on which we can then apply the Area Under Baseline.

    1. (a)

We can visualise the median algorithm outputs per each state on a map, similar
to those shown toward the beginning of the paper.

    2. (b)

With enough tweaking, we could use this algorithm to predict the effect of the
2020 recession on ZHVI values. This will involve more research on how to work
with non-stationary time data, and may possibly yield the use of another more
flexible time-series module.

  3. 3.

Comparing methodologies of other researchers in acquiring results for
recession effects on housing metrics, and seeing if the results we attained
stack up. This will provide us more insight as to what machine-learning
approach we should employ with the sort of data we are given.

## References

  * [1] The great recession, 2019.
  * [2] Kimberly Amadeo.  What caused the 2008 financial crisis and could it happen again?, 2019\.
  * [3] Richie Bernardo.  2017’s most and least recession-recovered cities, Jan 2017.
  * [4] CityLab, University of Toronto’s School of Cities, and Rotman School of Management.  Which u.s. cities suffer the most during a recession?, Jun 2016.
  * [5] Investopedia.  How the 2008 housing crash affected the american dream, Nov 2019.
  * [6] Robert Rich.  The great recession, 2019.
  * [7] Lauryn Ringwood, Philip Watson, and Paul Lewin.  A quantitative method for measuring regional economic resilience to the great recession.  Growth and Change, 50(1):381–402, 2019.


††thanks<EMAIL_ADDRESS>
# Biohybrid membrane formation by directed insertion of Aquaporin in to a
solid-state nanopore

François Sicard  A. Ozgur Yazaydin  Department of Chemical Engineering,
University College London, WC1E 7JE London, UK

###### Abstract

Technical challenges in molecule sensing and chemical detection have created
an increasing demand for transformative materials with high sensitivity and
specificity. Biohybrid nanopores have attracted growing interest as they can
ideally combine the durability of solid-state nanopores with the precise
structure of biological nanopores. Particular care must be taken to control
how biological nanopores adapt to their surroundings once in contact with the
solid-state nanopore. Two major challenges are to precisely control this
adaptability under dynamic conditions and provide predesigned functionalities
that can be manipulated for engineering applications. Here, we report on the
computational design of a distinctive class of biohybrid active membrane
layer, built from the directed-insertion of aquaporin-incorporated lipid shell
in to a silica nanopore. First, we describe in detail the mechanisms at play
in the insertion of the biological membrane in to the solid-state nanopore.
Then we analyze the structural stability of the system and demonstrate that
its water permeability is comparable to the one measured in the biological
environment. Finally, we discuss how the technology implemented could be
applicable to environmental and biomedical applications, such as water
desalination and drug discovery, where targeting and controlled permeation of
small molecules must be efficiently addressed.

biohybrid nanopore, aquaporin, nanodisc, directed insertion, permeability,
molecular dynamics simulation

Considerable attention has been devoted to the design, characterization and
development of biohybrid nanopore based platforms due to their potential
application in the areas of environment and healthcare including molecule
sensing Cao et al. (2019); Nicolaï et al. (2020), disease diagnosis Minervini
et al. (2020), drug design Wang et al. (2020), chemical detection Shoji et al.
(2019), pollutant removal Esmaeilian et al. (2015), and water desalination
Ramanathan et al. (2018). In particular, biohybrid membranes have received
growing interest over the last decade as they can ideally combine the
durability of solid-state nanopores with the precise structure and
functionality of biological nanopores Haque et al. (2013). In this
perspective, tremendous progress has been made regarding the incorporation of
biological compounds inside synthetic pores by modifying their surface
chemistry. This includes, among others, the coating of nanopores with
supported fluid lipid bilayer Yusko et al. (2011); Sun et al. (2013); Wang et
al. (2013); Ding et al. (2015); Fuwad et al. (2019) and the insertion of
biological structures in silicon based substrates such as the
$\alpha$-haemolysin protein Hall et al. (2010), DNA origami structures Bell et
al. (2012), and biological nanodiscs Farajollahi et al. (2018).

We report here on the computational design of a distinctive class of biohybrid
active membrane layer built from the directed insertion of aquaporin-
incorporated lipid shell in a model silica pore, which could pave the way to
the development of transformative devices for water filtration and drug
discovery. Aquaporins (Aqp) are transmembrane proteins that are present in the
biological membranes of organisms such as plants, mammals and bacteria Verkman
(2013). They form pores that act as water channels allowing water molecules to
move across the plasma membrane while rejecting protons, charged particles and
other solutes. Of particular interest is the bacterial Aquaporin Z (AqpZ),
expressed in Escherichia coli, which is the smallest, simplest and most robust
member in the Aqp family Calamita (2000). Along with its high selectivity and
water permeability, AqpZ has been purified in high concentration allowing it
to be used for integration into biomimetic membranes for water filtration Sun
et al. (2013); Hang et al. (2016).

In this work, AqpZ is originally embedded in a lipid nanodisc stabilized by
membrane scaffold proteins (MSPs), as shown in Fig. 1. MSPs are charged
helical-amphiphatic proteins that wrap around the lipid bilayers with
precisely defined inner and outer diameters Ritchie et al. (2009). MSP-based
lipid nanodiscs are known to provide a native membrane environment that aids
in preparing integral membrane proteins in biologically active folded forms
for structural studies Denisov and Sligar (2016). They have recently been
shown to serve as biological template that can be electrochemically inserted
into solid-state nanopore Farajollahi et al. (2018). This integration could be
a promising step towards reproducible mass production of biohybrid membranes.
However, it still faces structural limitations associated with the positioning
of the MSP-based nanodisc and its interaction with the solid-state nanopore,
which affects the sealing performance of the membrane and limits its
application in biohybrid nanopore technology.

To circumvent this limitation, we show how the AqpZ-incorporated lipid
nanodisc can be dissociated from its surrounding MSPs and inserted into
functionalized silica nanopore to form a functional biohybrid membrane layer.
We first employ non-equilibrium molecular dynamics (NEMD) simulation to study
the dissociation of the AqpZ-incorporated lipid nanodisc from the immobilised
MSPs under the application of an external pressure gradient and its subsequent
insertion into the pore. In this process, we show that the thickness of the
solid-state nanopore and the nature of the functionalizing groups must be
properly chosen with respect to the characteristics of the MSP-based lipid
nanodisc to ensure the final biohybrid nanopore is stable and it is leak proof
at the interface of the functionalized nanopore and the lipid shell. In
particular, the preferential interaction between the hydrocarbon tail of the
lipid molecules and the interior surface of the pore yields the expulsion of
the water molecules and therefore is responsible for tuning the sealing
performance of the membrane.

Figure 1: Schematic representation of (a) the AqpZ-incorporated lipid nanodisc
stabilized by the MSP1E2D1 membrane scaffold protein and (b) the cylindrical
nanopore carved out of a silica-cristobalite slab. The external and internal
surface of the slab is functionalized with hydroxyl and propyl groups,
respectively. Half of the nanopore is made transparent and half of the lipid
bilayer is omitted for clarity. The backbone of the AqpZ monomer and the MSP
are shown in cyan and dark grey, respectively. Yellow, red, blue, white, gold,
and cyan spheres represent silicon, oxygen, nitrogen, hydrogen, phosphate and
carbon atoms, respectively.

We then investigate the structural stability of the biohybrid membrane layer
in pure and saline water conditions with equilibrium MD simulations. Finally,
we study the non-equilibrium osmotic permeability of the system when an
external pressure gradient is applied and compare the results with those
measured in biological environments.

## Results and Discussion

System Characteristics. We first give the characteristics of the MSP-based
lipid nanodisc and solid-state nanopore considered in this work. A schematic
representation of the system is shown in Fig. 1. One AqpZ monomer is
originally embedded in a $1$-palmitoyl-$2$-oleoyl-sn-
glycero-$3$-phosphocholine (POPC) lipid-nanodisc stabilized by the MSP1E2D1
membrane scaffold protein. The diameter of the corresponding nanodisc is
$11.1$ nm Errasti-Murugarren et al. (2021). A cylindrical nanopore is carved
out of a silica-cristobalite slab which surface was functionalized with
hydroxyl groups (-OH). The interior of the pore is functionalized with propyl
groups ($\textrm{-C}_{3}\textrm{H}_{7}$) to render its surface hydrophybic.
The length of the functionalizing group is specifically chosen such that it
confers adequate sealing properties of the nanopore due to preferential
interaction between the hydrocarbon tails of the lipid shell and the
functionalizing groups. The comparison with shorter groups is shown in the
Supporting Information (SI). The height and diameter of the cylindrical
nanopore were specifically chosen to allow the insertion of the AqpZ-
incorporated lipid shell, while preventing the insertion of the MSPs. The
CHARMM force field Huang et al. (2016) was employed to model the biological
entities (MSPs, AqpZ, POPC) and the propyl group on the silica surface. The
INTERFACE force field Heinz et al. (2013) and TIP3P potential Jorgensen et al.
(1983) were employed to model the silica nanopore and water molecules,
respectively, as they are compatible with the CHARMM force field. The
simulations were conducted using the GROMACS software package, version 2018.8
Abraham et al. (2015). The energy of the system was first optimized using the
steepest descent minimization algorithm. The system was then equilibrated at
constant temperature, $T=300$ K, finishing with equilibration at constant
pressure, $P=1$ bar, to equilibrate the fluid density. Finally, a production
simulation was carried out in the NVT ensemble. NEMD simulations were
performed using the extension of the pull module of the GROMACS software
developed by Gräter and coworkers Herrera-Rodriguez et al. (2019). The details
of the numerical parameterization and the system equilibration are given in
the SI.

Figure 2: Schematic representation of the insertion of the AqpZ-incorporated
lipid nanodisc into the solid-state nanopore. (a) Under the application of the
hydrostatic pressure difference, the lipid nanodisc first comes in contact
with the silica slab. (b) As the system moves further down the pore, the outer
edge of the lipid nanodisc acts as a support layer that stabilizes the
membrane scaffold proteins (MSPs) and prevents them to penetrate inside the
cavity. (c) The hydrophobic tails of the lipid molecules interact favorably
with the propyl groups functionalizing the interior surface of the pore, which
yields the depletion of the water molecules. The backbone of the MSPs is shown
in dark grey. Yellow, red, blue, white, gold, and cyan spheres represent
silicon, oxygen, nitrogen, hydrogen, phosphate and carbon atoms in the lipid
molecules and silica surface, respectively. The oxygen atoms of the water
molecules are shown in blue. Only half of the pore and the water molecules in
the pore are shown for clarity.

Directed insertion of AqpZ-incorporated lipid nanodisc in to a solid-state
nanopore. The MSP-based lipid nanodisc is initially positioned above the
cylindrical pore, as shown in Fig. 2. Farajollahi et al. have recently
demonstrated that the positioning of MSP-based nanodiscs in solid state
nanopores could be effectively driven by applying concentration gradient-based
or electrophoresis-driven diffusion Farajollahi et al. (2018). A water
pressure gradient is applied across the membrane to dissociate the AqpZ-
incorporated lipid shell from the MSPs and to direct its insertion in the
cylindrical pore. Using the numerical method of Zhu et al. Zhu et al. (2002,
2004), a hydrostatic pressure difference of $\Delta P\sim 80$ MPa is created
across the pore and is applied for $10$ ns, while the position of the solid-
state membrane was restrained. The value of $\Delta P$ was chosen to achieve
the insertion of the nanodisc in a reasonably short simulation time. For
comparison, the free energy profile associated with the equilibrium
dissociation of the AqpZ monomer from the MSP-based lipid nanodisc is shown in
the SI.

Under the application of the hydrostatic pressure difference, the lipid
nanodisc first comes in contact with the silica slab. Due to preferential
interaction between the hydrophilic head of the lipid molecules and the
hydroxyl group on the surface of the slab, the outer edge of the lipid
nanodisc acts as a support layer that stabilizes the MSPs and prevents them to
penetrate inside the cavity (see Fig 2a). As the AqpZ-incorporated lipid shell
penetrates further in the solid-state membrane, the hydrophobic tails of the
lipid molecules start interacting favorably with the propyl groups. This
preferential interaction yields the depletion of the water molecules in the
pore. The depletion of water does not evenly happen around the interior
surface of the pore at first, as water molecules preferentially interact with
the hydrophilic head of the lipid molecules (see Fig 2b). As the simulation
progresses, the AQP-incorporated lipid shell gradually moves further down the
pore under the effect of the hydrostatic pressure difference. As a result, the
inserted lipid molecules evenly spreads in the solid-state membrane, leading
to the complete depletion of water around the Aqp-incorporated lipid shell
(see Fig 2c). The water pressure difference across the membrane is then
switched off and the system is relaxed in a $\sim 200$ ns simulation.

Stability of the active layer in pure water. We first study the structural
stability of the biohybrid membrane layer in pure water, i.e. in the absence
of ions. A representative snapshot of the system after relaxation is shown in
Fig. 3. We observe the increase of the thickness of the bilayer as the lipid
molecules radially position away from the center of the protein. This
evolution is quantitatively assessed in Fig. 4, where the axial position of
the lipid molecules is shown as a function of their radial position with
respect to the center of the pore. The average value of the position
distribution of the lipid molecules measured in the biological environment,
corresponding to a membrane thickness $h_{\textrm{bio}}\sim 2.5$ nm, is shown
for comparison. As the lipid molecules position closer to the surface of the
pore, the thickness of the bilayer progressively increases, departing from the
thickness measured in the biological environment. This structural deformation
is associated with the configurational change of the hydrocarbon chain of the
lipid molecules with a transition from a Gauche conformation near the center
of the nanopore to a nearly all-trans conformation at the pore surface, as
shown in Fig. 4.

Figure 3: Schematic representation of the Aqp-incorporated lipid shell
inserted in the silica nanopore after the system relaxed for $200$ ns. The
backbone of the AqpZ monomer is shown in cyan. Yellow, red, blue, white, gold,
and cyan spheres represent silicon, oxygen, nitrogen, hydrogen, phosphate and
carbon atoms, respectively. The water molecules are not shown and the atoms
composing the functionalized silica pore are shadowed for clarity. We observe
the progressive increase of the thickness of the bilayer as the lipid
molecules radially position away from the center of the protein.

This transition is due to preferential Van der Waals interaction between the
hydrophobic tail of lipids and the hydrophobic surface functionalized with the
propyl groups F.Sicard et al. (2018); Sicard and Striolo (2021).

Figure 4: Distribution of the lipid molecules within the ApqZ-incorporated
lipid shell inserted into the solid-state nanopore, as a function of their
radial and axial position. The mean value of the position distribution of the
lipid molecules in the biological environment is shown for comparison. Both
the positions of the upper and lower leaflets are shown in pure (green) and
saline (red) water. The thickness of the bilayer progressively increases as
the lipid molecules radially position away from the center of the pore. The
AqpZ monomer and the interior surface of the pore are positioned at $\rho<2$
nm and $\rho\sim 4$ nm, respectively. Uncertainties, defined as the standard
error, are represented by the shaded area.

To further assess the stability of the biohybrid membrane and its ability to
serve as a separator, we measured in Fig. 5 the density distribution of the
water molecules centered around the AqpZ monomer in the simulation box after
relaxation The density distribution of water measured in the biological
environment is shown for comparison. For $|z|\leq 1$, the profile of the water
density is similar to the one observed in the biological simulation. For
$1<|z|\leq 4$ we observe the progressive increase of the water density, which
happens at a slower rate than in the biological environment. This difference
is explained by the conformational organisation of the lipid bilayer in the
nanopore due to the preferential interaction with the hydrophobic surface, as
shown in Fig. 3. For $|z|>4$, i.e. outside the solid-state nanopore, the
system transitions abruptly to reach the bulk water density similar to the one
measured in the biological environment.

Figure 5: Density profile of the water molecules along the Z direction of the
simulation box and centered around the ApqZ monomer. Following the directed
insertion of the Aqp-incorporated nanodisc into the solid-state membrane, we
let the system relax for $200$ ns in fresh (green) and sea (red) water. The
density profile of the water molecules in the biological environment is shown
for comparison (blue). Uncertainties, defined as the standard error, are
represented by the shaded area.

Stability of the active layer in sea saline water. We extend the analysis
above to the study of the structural stability of the active membrane layer in
the presence of chloride and sodium ions at characteristic salinity of sea
water. In this condition, about $3.5\%$ of the weight of seawater comes from
the dissolved salts Cullum et al. (2016). The evolution of the thickness of
the biological membrane as a function of the radial position of the lipid
molecules is shown in Fig. 4. We observe the increase of the thickness of the
bilayer, similar to the evolution obtained in pure water, as the lipid
molecules radially position away from the center of the protein. In addition,
we show in Fig. 5 the density distribution of water molecules in the
simulation box. The density profile is similar to the one measured in pure
water. This indicates that the presence of ions does not affect the structural
stability of the Aqp-incorporated lipid shell once inserted in the solid-state
nanopore. Furthermore, we studied the density distribution of ions in the
simulation box after relaxation. The results are shown in Fig. S3a in the SI.
Visual inspection shows that the ions do not penetrate the water channel of
AqpZ or the lipid shell. This is quantitatively assessed in Fig. S3b in the
SI, where we measured the density distribution of ions in the simulation box.
In particular, the ion density for $|z|\leq 1$ is equal to zero. For
$1<|z|\leq 4$ we observe the progressive increase of the ion density due to
the increase of the thickness of the biological membrane, similar to what we
obtained for the density distribution of the water molecule. Finally, for
$|z|>4$, the system transitions abruptly to reach the bulk ion density.

Water transport across the active membrane layer. We complete the analysis
above with the study of the osmotic permeability of the biohybrid active
membrane layer. In the condition of reverse osmosis, the ability of the system
to conduct water is characterized by the ratio of the net water flow to the
hydrostatic pressure difference, $\Delta P$, across the membrane. The volume
flux, $J_{v}$, defined as the net flow of water per unit area of the membrane
is obtained as

$J_{v}=L_{P}\Delta P$ (1)

with $L_{P}$ the hydraulic permeability of the membrane. For comparaison, when
the two sides of the membrane have the same hydrostatic pressure but different
concentrations of an impermeable solute, an osmotic pressure difference is
established and water will flow from the side with lower solute concentration
to the other side. In dilute solutions, the molar water flux, $J_{m}$, is
linearly proportional to the solute concentration difference $\Delta C$

$J_{m}=P_{f}\Delta C$ (2)

with $P_{f}$ the osmotic permeability of the membrane. The water flux
generated due to a solute concentration difference is identical to that of
generated by a hydrostatic pressure difference $\Delta P=RT\Delta C$, with $R$
the gas constant and $T$ the temperature of the system. Therefore, $L_{P}$ and
$P_{f}$ are related by a constant factor

$P_{f}=(RT/V_{W})L_{P}\,,$ (3)

with $V_{W}$ the molar volume of water.

To determine the hydraulic permeability of the active membrane layer, a water
pressure gradient is applied along the direction orthogonal to the membrane.
The position of the membrane is restrained along the direction of the pressure
gradient to avoid undesirable translation of the system. Details of the
numerical implementation is given in the SI. The water flux, $J_{W}$, through
the channel in the membrane can then be measured at the atomistic level by
counting the net number of water molecules passing through the channel during
the simulation.

We consider five values for the hydrostatic pressure difference ranging from
$10$ MPa to $100$ MPa, as shown in Fig. 6. These values are sufficiently high
to asses the measure of the water flux in a reasonable simulation time and
sufficiently low to avoid the accumulation of water molecules in the vestibule
of AqpZ, which could affect the linear relation between water flux and
hydrostatic pressure difference Zhu et al. (2002). For comparaison, the
osmotic pressure of physiological solutions is usually below $10$ MPa Hendry
(1962) whereas the osmotic pressure measured for sea water is $\sim 26$ bars
and the pressure traditionally used in reverse osmosis desalination of sea
water is $\sim 60$ bars Davenport et al. (2018).

Figure 6: Dependence of the water flux on the applied pressure difference.
Error bars represent the standard deviations of the water flux. The line
corresponding to the best-fit slope for the five data points is shown.

The evolution of the volume flux, $J_{v}$, across the active membrane as a
function of the hydrostatic pressure difference, $\Delta P$, is shown in Fig.
6. Each NEMD simulation was run for $\sim 70$ ns and the hydraulic
permeability, $L_{P}$, of the membrane is determined by interpolating the
best-fit slope through the data points. Importantly, we observe that the
linear relation between the water flux and the hydrostatic pressure
difference, as given in Eq. 1, holds. We then measure the value for the
hydraulic permeability $L_{P}=7.3\pm 0.5\times
10^{-18}\textrm{cm}^{5}~{}\textrm{N}^{-1}~{}\textrm{s}^{-1}$. Applying Eq. 3
with $T=300K$ we obtain the osmotic permeability $P_{f}=1.1\pm 0.1\times
10^{-13}\textrm{cm}^{3}~{}\textrm{s}^{-1}$. The numerical result is in
excellent agreement with the experimental estimate of the osmotic water
permeability of the AqpZ monomer measured in biological conditions,
$P^{\textrm{exp}}_{f}\sim 1.0\times 10^{-13}\textrm{cm}^{3}/s$ Borgnia et al.
(1999). Very importantly, this result confirms that the biological
functionality of the AqpZ monomer is preserved after its directed insertion
into the silica nanopore.

## Conclusions

The extensive numerical simulations discussed above allowed us to explore the
design and stability of a distinctive class of biohybrid active membrane layer
built from the directed insertion of aquaporin-incorporated lipid nanodiscs
into a model silica nanopore. We described in details the mechanisms at play
during the insertion of the biological membrane into the solid-state nanopore.
In particular, we showed that the preferential interaction between the
hydrophobic chain of the POPC lipid bilayer and the alkyl group
functionalizing the interior surface of the pore played an essential role both
in the depletion of water molecules in the pore and the sealing properties of
the system. We numerically assessed the osmotic permeability of the biohybrid
membrane and compared it with the experimental permeability of the AqpZ
monomer in the biological environment. In particular, we showed that the
biological function of the AqpZ monomer, in terms of water transport, was
preserved after the directed insertion into the silica nanopore. Experimental
verification is required to test our expectations.

The biohybrid nanopore technology demonstrated in this work can be
instrumental to the development of robust, thin, and defect-free membranes,
with high selectivity and permeability. In particular, incorporating aquaporin
proteins into compatible solid-state nanopores, while ensuring industrial
production scalability, will be decisive for translating the nano-scale
science into a large range of environmental and biomedical applications,
including water desalination and aquaporin drug discovery.

In the former case, membrane technology has been dominating the water
desalination industry for decades due to its high efficiency and reliability
Pedersen et al. (2018). However, the experimental performance of biomimetic
desalination membranes remains far below that of biological membranes Porter
et al. (2020). The incorporation of aquaporin proteins into compatible solid-
state nanopores provides biohybrid technology with the potential for site-
specific genetic engineering or chemical modifications, which could
potentially be used as active membrane layer in both reverse osmosis and
forward osmosis applications to meet the increasing demand for fresh water at
lower energy consumption and operating costs.

In the latter case, aquaporins are attractive targets for the development of
novel drug therapies for life-threatening disorders, for example, kidney
disease, congestive heart failure, stroke, traumatic injury and tumor-induced
brain swelling, which currently lack adequate medical treatment Castle (2005).
Identifying new drug development frameworks for conditions associated with
aberrant water movement will meet the urgent clinical need of millions of
patients for whom no pharmacological interventions currently exist. A key
issue contributing to this situation may be the limited number of in vitro
assays suitable for high-throughput screening (HTS) of the pharmacological
regulation of aquaporin function Salman et al. (2022). In this perspective,
the possibility of stabilizing aquaporin-incorporated lipid nanodiscs in
compatible solid-state nanopore might be a promising way towards a
reproducible mass production of robust membranes with a multitude of identical
nanopores. This could pave the way to the design of transformative HTS drug
discovery platforms to advance aquaporin-targeted therapeutics.

## Acknowledements

This work was supported by the UK Engineering and Physical Sciences Research
Council (EPSRC) under grant number EP/V04804X/1. We are grateful to the UK
Materials and Molecular Modelling Hub for computational resources, which is
partially funded by the EPSRC under grant numbers EP/P020194/1 and
EP/T022213/1. FS thanks Kai Bin Yu, Turan Selman Erkal, and Fan Li for useful
discussions.

## References

  * Cao et al. (2019) Cao, C.; Cirauqui, N.; Marcaida, M.; Buglakova, E.; Duperrex, A.; Radenovic, A.; Peraro, M. D. Single-molecule sensing of peptides and nucleic acids by engineered aerolysin nanopores. _Nat. Commun._ 2019, _10_ , 4918.
  * Nicolaï et al. (2020) Nicolaï, A.; Rath, A.; Delarue, P.; Senet, P. Nanopore sensing of single-biomolecules: a new procedure to identify protein sequence motifs from molecular dynamics. _Nanoscale_ 2020, _12_ , 22743–22753.
  * Minervini et al. (2020) Minervini, C.; Cumbo, C.; Orsini, P.; Anelli, L.; Specchia, A. Z. G.; Albano, F. Nanopore Sequencing in Blood Diseases: A Wide Range of Opportunities. _Front. Genet._ 2020, _11_ , 76.
  * Wang et al. (2020) Wang, X.; Wilkinson, M.; Lin, X.; Ren, R.; Willison, K.; Ivanov, A.; Baum, J.; Edel, J. Single-molecule nanopore sensing of actin dynamics and drug binding. _Chem. Sci._ 2020, _11_ , 970–979.
  * Shoji et al. (2019) Shoji, K.; Kawano, R.; White, R. Spatially Resolved Chemical Detection with a Nanoneedle-Probe-Supported Biological Nanopore. _ACS Nano_ 2019, _13_ , 2606–2614.
  * Esmaeilian et al. (2015) Esmaeilian, A.; Mazdeh, A. M.; H. Ghaforian, A. L. A novel nanopore biopolymer multi adsorbent for simultaneous removal of anionic and cationic mixtures. _Desalination Water Treat._ 2015, _53_ , 2235–2248.
  * Ramanathan et al. (2018) Ramanathan, A.; Aqra, M.; Al-Rawajfeh, A. Recent advances in 2D nanopores for desalination. _Environ. Chem. Lett._ 2018, _16_ , 1217–1231.
  * Haque et al. (2013) Haque, F.; Li, J.; Wu, H.-C.; Liang, X.-J.; Guo, P. Solid-state and biological nanopore for real-time sensing of single chemical and sequencing of DNA. _Nano Today_ 2013, _8_ , 56–74.
  * Yusko et al. (2011) Yusko, E.; Johnson, J.; Majd, S.; Prangkio, P.; Rollings, R.; Li, J.; Yang, J.; Mayer, M. Controlling protein translocation through nanopores with bio-inspired fluid walls. _Nature Nano._ 2011, _6_ , pages 253–260.
  * Sun et al. (2013) Sun, G.; Chung, T.-S.; Jeyaseelan, K.; Armugam, A. Stabilization and immobilization of aquaporin reconstituted lipid vesicles for water purification. _Colloids and Surfaces B: Biointerfaces_ 2013, _102_ , 466–471.
  * Wang et al. (2013) Wang, H.; Chung, T.-S.; Tong, Y.; Jeyaseelan, K.; Armugam, A.; Duong, H. P.; Fu, F.; Seah, H.; Yang, J.; Hong, M. Mechanically robust and highly permeable AquaporinZ biomimetic membranes. _J. Membrane Sci._ 2013, _1_ , 130–136.
  * Ding et al. (2015) Ding, W.; Cai, J.; Yu, Z.; Wang, Q.; Xu, Z.; Wang, Z.; Gao, C. Fabrication of an aquaporin-based forward osmosis membrane through covalent bonding of a lipid bilayer to a microporous support. _J. Mater. Chem. A_ 2015, _3_ , 20118–20126.
  * Fuwad et al. (2019) Fuwad, A.; Ryu, H.; Lee, J.-H.; Kim, D.; Yoo, Y.-E.; Kim, Y.-R.; Kim, S.; Jeon, T.-J. An electrokinetic approach to fabricating aquaporin biomimetic membranes for water purification. _Desalination_ 2019, _452_ , 9–16.
  * Hall et al. (2010) Hall, A.; Scott, A.; Rotem, D.; Mehta, K.; Bayley, H.; Dekker, C. Hybrid pore formation by directed insertion of $\alpha$-haemolysin into solid-state nanopores. _Nature Nano._ 2010, _5_ , 874–877.
  * Bell et al. (2012) Bell, N.; Engst, C.; Ablay, M.; Divitini, G.; Ducati, C.; Liedl, T.; Keyser, U. DNA Origami Nanopores. _Nano Lett._ 2012, _12_ , 512–517.
  * Farajollahi et al. (2018) Farajollahi, F.; Seidenstücker, A.; Altintoprak, K.; Walther, P.; Ziemann, P.; Plettl, A.; Marti, O.; Wege, C.; Gliemann, H. Electrochemically-Driven Insertion of Biological Nanodiscs into Solid State Membrane Pores as a Basis for “Pore-In-Pore” Membranes. _Nanomaterials_ 2018, _8_ , 237.
  * Verkman (2013) Verkman, A. Aquaporins. _Curr. Biol._ 2013, _23_ , 52–55.
  * Calamita (2000) Calamita, G. The Escherichia coli aquaporin-Z water channel. _Mol Microbiol._ 2000, _37_ , 254–262.
  * Hang et al. (2016) Hang, B.; Pan, J.; Ni, D.; Zheng, Q.; Zhang, X.; Cai, J.; Huang, L.; Wei, P.; Xu, Z. High-level production of aquaporin Z in Escherichia coli using maltose-binding protein/polyhistidine dual-affinity tag fusion system. _Process Biochemistry_ 2016, _51_ , 599–606.
  * Ritchie et al. (2009) Ritchie, T.; Grinkova, Y.; Bayburt, T.; Denisov, I.; Zolnerciks, J.; Atkins, W.; ; Sligar, S. Reconstitution of Membrane Proteins in Phospholipid Bilayer Nanodiscs. _Methods Enzymol._ 2009, _464_ , 211–231.
  * Denisov and Sligar (2016) Denisov, I.; Sligar, S. Nanodiscs for structural and functional studies of membrane proteins. _Nat. Struct. Mol. Biol._ 2016, _23_ , 481–486.
  * Errasti-Murugarren et al. (2021) Errasti-Murugarren, E.; Bartoccioni, P.; Palacín, M. Membrane Protein Stabilization Strategies for Structural and Functional Studies. _Membranes_ 2021, _11_ , 155.
  * Huang et al. (2016) Huang, J.; Rauscher, S.; Nawrocki, G.; Ran, T.; Feig, M.; de Groot, B.; Grubmuller, H.; Jr., A. M. CHARMM36m: An Improved Force Field for Folded and Intrinsically Disordered Proteins. _Nature Methods_ 2016, _14_ , 71–73.
  * Heinz et al. (2013) Heinz, H.; Lin, T.-J.; Mishra, R. K.; Emami, F. Thermodynamically Consistent Force Fields for the Assembly of Inorganic, Organic, and Biological Nanostructures: The INTERFACE Force Field. _Langmuir_ 2013, _29_ , 1754–1765.
  * Jorgensen et al. (1983) Jorgensen, W.; Chandrasekhar, J.; Madura, J. Comparison of simple potential functions for simulating liquid water. _J. Chem. Phys._ 1983, _79_ , 926.
  * Abraham et al. (2015) Abraham, M.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J.; Hess, B.; Lindahl, E. GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. _SoftwareX_ 2015, _1-2_ , 19–25.
  * Herrera-Rodriguez et al. (2019) Herrera-Rodriguez, A.; Miletić, V.; Aponte-Santamaria, C.; Gräter, F. Molecular Dynamics Simulations of Molecules in Uniform Flow. _Biophys. J._ 2019, _116_ , 1579–1585.
  * Zhu et al. (2002) Zhu, F.; Tajkhorshid, E.; Schulten, K. Pressure-Induced Water Transport in Membrane Channels Studied by Molecular Dynamics. _Biophys. J._ 2002, _83_ , 154–160.
  * Zhu et al. (2004) Zhu, F.; Tajkhorshid, E.; Schulten, K. Theory and Simulation of Water Permeation in Aquaporin-1. _Biophys. J._ 2004, _86_ , 50–57.
  * F.Sicard et al. (2018) F.Sicard,; Bui, T.; Monteiro, D.; Lan, Q.; Ceglio, M.; Burress, C.; A.Striolo, Emergent Properties of Antiagglomerant Films Control Methane Transport: Implications for Hydrate Management. _Langmuir_ 2018, _34_ , 9701–9710.
  * Sicard and Striolo (2021) Sicard, F.; Striolo, A. Role of structural rigidity and collective behaviour in the molecular design of gas hydrate anti-agglomerants. _Mol. Syst. Des. Eng._ 2021, _6_ , 713–721.
  * Cullum et al. (2016) Cullum, J.; Stevens, D.; Joshi, M. Importance of ocean salinity for climate and habitability. _Proc. Natl. Acad. Sci. U.S.A._ 2016, _113_ , 4278–4283.
  * Hendry (1962) Hendry, E. The Osmotic Pressure and Chemical Composition of Human Body Fluids. _Clin. Chem._ 1962, _8_ , 246–265.
  * Davenport et al. (2018) Davenport, D.; Deshmukh, A.; Werber, J.; Elimelech, M. High-Pressure Reverse Osmosis for Energy-Efficient Hypersaline Brine Desalination: Current Status, Design Considerations, and Research Needs. _Environ. Sci. Technol. Lett._ 2018, _5_ , 467–475.
  * Borgnia et al. (1999) Borgnia, M.; Kozonoand, D.; Calamita, G.; Maloney, P.; Agre, P. Functional reconstitution and characterization of AqpZ, the E. coli water channel protein. _J. Mol. Biol_ 1999, _291_ , 1169–1179.
  * Pedersen et al. (2018) Pedersen, P.; Bjorkskov, F.; Alvisse, S.; Hélix-Nielsen, C. From channel proteins to industrial biomimetic membrane technology. _Faraday Discuss._ 2018, _209_ , 287–301.
  * Porter et al. (2020) Porter, C.; Werber, J.; Zhong, M.; Wilson, C.; Elimelech, M. Pathways and Challenges for Biomimetic Desalination Membranes with Sub-Nanometer Channels. _ACS Nano_ 2020, _14_ , 10894–10916.
  * Castle (2005) Castle, N. Aquaporins as targets for drug discovery. _Drug Discov._ 2005, _10_ , 485–493.
  * Salman et al. (2022) Salman, M.; P.Kitchen, .; Yool, A.; Bill, R. Recent breakthroughs and future directions in drugging aquaporins. _Trends Pharmacol. Sci._ 2022, _43_ , 30–42.
  * Savage et al. (2003) Savage, D.; Egea, P.; Robles-Colmenares, Y.; III, J. O.; Stroud, R. Architecture and Selectivity in Aquaporins: 2.5 ÅX-Ray Structure of Aquaporin Z. _PLoS Biol._ 2003, _1_ , e72.
  * Kelley et al. (2015) Kelley, L.; Mezulis, S.; Yates, C.; Wass, M.; Sternberg, M. The Phyre2 web portal for protein modeling, prediction and analysis. _Nat. Protoc._ 2015, _10_ , 845–858.
  * Humphrey et al. (1996) Humphrey, W.; Dalke, A.; Schulten, K. VMD - visual molecular dynamics. _J. Mol. Graph._ 1996, _14_ , 33–38.
  * Ferré and Clote (2005) Ferré, F.; Clote, P. DiANNA: a web server for disulfide connectivity prediction. _Nuc. Acids Res._ 2005, _33_ , 230–232.
  * Qi et al. (2019) Qi, Y.; Lee, J.; Klauda, J.; Im, W. CHARMM-GUINanodisc Builderfor Modeling and Simulationof Various Nanodisc Systems. _J. Comput. Chem._ 2019, _40_ , 893–899.
  * Jo et al. (2008) Jo, S.; Kim, T.; Iyer, V.; Im, W. CHARMM-GUI: A web-based graphical user interface for CHARMM. _J. Comput. Chem._ 2008, _11_ , 1859–1865.
  * (46) BIOVIA, Dassault Systémes, U., San Diego Materials Studio. https://hadoop.apache.org.
  * Emami et al. (2014) Emami, F.; Puddu, V.; Berry, R.; Varshney, V.; Patwardhan, S.; Perry, C.; Heinz, H. Force Field and a Surface Model Database for Silica to Simulate Interfacial Properties in Atomic Resolution. _Chem. Mater._ 2014, _26_ , 2647–2658.
  * Tribello et al. (2014) Tribello, G.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. PLUMED 2: New feathers for an old bird. _Comput. Phys. Comm._ 2014, _185_ , 604–613.
  * Hess et al. (1997) Hess, B.; Bekker, H.; Berendsen, H.; Fraaije, J. G. Lincs: A linear constraint solver for molecular simulations. _J. Comput. Chem._ 1997, _98_ , 1463.
  * Darden et al. (1993) Darden, T.; York, D.; Pedersen, L. _J. Chem. Phys._ 1993, _98_ , 10089–10092.
  * Sicard et al. (2019) Sicard, F.; Toro-Mendoza, J.; Striolo, A. Nanoparticles Actively Fragment Armored Droplets. _ACS Nano_ 2019, _13_ , 9498–9503.
  * Paci and Karplus (1999) Paci, E.; Karplus, M. Forced Unfolding of Fibronectin Type 3 Modules: An Analysis by Biased Molecular Dynamics Simulations. _J. Mol. Biol._ 1999, _288_ , 441–459.
  * Marchi and Ballone (1999) Marchi, M.; Ballone, P. Adiabatic bias molecular dynamics: A method to navigate the conformational space of complex molecular systems. _J. Chem. Phys._ 1999, _110_ , 3697–3702.
  * Camilloni et al. (2011) Camilloni, C.; Broglia, R.; Tiana, G. Hierarchy of folding and unfolding events of protein G, CI2, and ACBP from explicit-solvent simulations. _J. Chem. Phys._ 2011, _134_ , 045105.
  * Sicard and Striolo (2016) Sicard, F.; Striolo, A. Numerical analysis of Pickering emulsion stability: insights from ABMD simulations. _Faraday Discuss._ 2016, _191_ , 287–304.
  * Kastner (2011) Kastner, J. Umbrella sampling. _WIREs Comput. Mol. Sci._ 2011, _1_ , 932–942.
  * Rosta and Hummer (2015) Rosta, E.; Hummer, G. Free Energies from Dynamic Weighted Histogram Analysis Using Unbiased Markov State Model. _J. Chem. Theory Comput._ 2015, _11_ , 276–285.

Biohybrid membrane formation by directed insertion of Aquaporin in to a solid-
state nanopore

Supporting Information

## System characteristics

Aquaporin Z monomer. The crystal structure of E. coli AqpZ was obtained from
the Protein Data Bank (PDB), entry 1RC2 Savage et al. (2003). The PDB file
contains two protomers, from which protomer A was used to build the system
simulated in this work. Missing side chains of residues Arg3, Glu31, Ser104,
Arg230, and Asp231, as well as missing hydrogen atoms, were added using the
programs Phyre2 Kelley et al. (2015) and VMD Humphrey et al. (1996). Titrable
side chains were simulated in their default titration state, i.e. Glu and Asp
residues with a negative charge, Lys and Arg residues with a positive charge,
and all other side chains with zero charge. Therefore, the protein is
electrically neutral and no counterions were needed to ensure the
electroneutrality of the AqpZ monomer. Finally, the program DiANNA Ferré and
Clote (2005) was used to check for the presence of any disulfide bond.

AqpZ-incorporated lipid nanodisc. The program Nanodisc Builder Qi et al.
(2019) in CHARMM-GUI Jo et al. (2008) was used to build the Aquaporin-
incorporated lipid nanodisc. One AqpZ monomer was embedded in a
$1$-palmitoyl-$2$-oleoyl-sn-glycero-$3$-phosphocholine (POPC) lipid-nanodisc
stabilized by the MSP1E2D1 membrane scaffold protein. The diameter and net
charge of the corresponding nanodisc are $11.1$ nm and $-14$e, respectively
Errasti-Murugarren et al. (2021). The number of POPC lipid molecules in the
discoidal lipid bilayer is $230$. The net charge of the MSP nanodisc was
neutralized with $14$ sodium ions.

Silica nanopore. The model solid-state nanopore was built from crystalline
silica. The model of unit cell of cristobalite was obtained from the shared
library of Materials Studio BIOVIA, Dassault Systémes , which was periodically
replicated in all directions. A regular silica slab was then carved out from
the bulk cristobalite silica. Unsaturated silica atoms were removed from the
surface of the slab and the non-bridging oxygen atoms were saturated with
hydroxyl (-OH) groups. The structure was processed to saturate all valences,
followed by energy minimization and molecular dynamics simulation to achieve
structural stability. The volume of the equilibrated silica slab is
$\sim(140\times 140\times 95)~{}\AA^{3}$. The resulting silica interface has a
density of silanol groups of $6.4$ OH/$\textrm{nm}^{2}$, within the range of
experimental data Emami et al. (2014). Next, a $90~{}\AA$ diameter cylindrical
void was created in the center of the silica slab. Unsaturated silica atoms
were removed from the interior surface of the pore and the non-bridging oxygen
atoms were saturated with alkyl group $-\textrm{C}_{n}\textrm{H}_{2n+1}$,
followed by energy minimization and molecular dynamics simulation to achieve
structural stability. We considered different functionalizing group with
$n=1$, $2$, and $3$. The ideal silica membrane model was adopted in this work
to study the effect of interaction between the biological membrane and the
pore surface on sealing property

## Simulation details

MD simulations were performed with the GROMACS software package, version
2018.8 Abraham et al. (2015). Biased simulations were performed using version
2.3 of the plugin for free energy (FE) calculations, PLUMED Tribello et al.
(2014). The CHARMM force field Huang et al. (2016) was employed to model the
biological entities (MSP, AqpZ, POPC) and the alkyl groups on the silica
surface. The INTERFACE force field Heinz et al. (2013) and TIP3P potential
Jorgensen et al. (1983) were employed to model the silica nanopore and water
molecules, respectively, as they are compatible with the CHARMM force field.
The time step used in all the simulations is $0.002$ ps, and the list of
neighbors is updated every $0.04$ ps with the grid method and a cutoff radius
of $1.2$ nm. The LINear Constraint Solver (LINCS) algorithm Hess et al. (1997)
handled bond constraints while the particle-mesh Ewald scheme Darden et al.
(1993) was used to treat long-range electrostatic interactions. The nonbonded
van der Waals cutoff radius was $1.2$ nm. The initial velocities were chosen
randomly. The cutoff algorithm was applied for the non-Coulomb potentials with
a radius of $1.2$ nm.

The energy of the system was first optimized using the steepest descent
minimization algorithm. The system was then equilibrated for $250$ ps at
constant temperature, $T=300$ K, using the Berendsen thermostat with lower
restraints. The equilibration phase was then conducted within the isobaric-
isothermal (NPT) ensemble to equilibrate the fluid density. The temperature
and pressure were maintained at $T=300$ K and $P=1$ bar, respectively, using
the Berendsen thermostat and barostat for $2$ ns. We then switched to the
Nose-Hoover thermostat and the Parrinello-Rahman barostat for an additional
$2$ ns, which are considered more thermodynamically consistent algorithms.
Finally the simulation was continued in NVT conditions coupling with the
v-rescale thermostat at constant temperature $T=300$ K.

Figure S1: Schematic representation of the numerical method implemented to
produce a pressure difference in MD simulations using periodic boundary
conditions. A constant force $f$ is exerted on the oxygen atoms of the water
molecules in region III. The membranes shown in the figure are the periodic
images of each other.

To achieve the insertion of the Aqp-incorporated lipid nanodisc into the
solid-state nanopore and to measure the osmotic permeability of the system, a
hydrostatic pressure difference was generated across the membrane using the
method proposed by Zhu et al. Zhu et al. (2002, 2004) and the extension of the
pull module of the GROMACS software developed by Gräter and coworkers Herrera-
Rodriguez et al. (2019). As shown in Fig. S1, three regions (I, II, III) were
defined in the water layer with region III isolated from the two sides of the
membrane by region I and II, respectively. A constant force $f$ along the
$z$-direction is exerted on all water molecules in region III, generating a
pressure gradient in this region. This results in a pressure difference
between regions I and II

$\Delta P=P_{1}-P_{2}=nf/A\,,$

where $n$ is the number of water molecules in region III, and $A$ is the area
of the membrane. To prevent the overall translation of the whole system along
the direction of the applied forces, position constraints were applied on the
membrane in the $z$-direction Zhu et al. (2004). In particular, we applied the
constant force $(f)$ on the oxygen atoms of the water molecules in region III,
defined as a $8$ nm thick layer. We applied harmonic constraints to the
$C_{\alpha}$ atoms of the proteins, the phosphorus atoms of the lipid
molecules, and the silicon and oxygen atoms of the silica pore, with spring
constants of $50~{}\textrm{kJ}/\textrm{mol}/\textrm{nm}^{2}$,
$700~{}\textrm{kJ}/\textrm{mol}/\textrm{nm}^{2}$, and
$700~{}\textrm{kJ}/\textrm{mol}/\textrm{nm}^{2}$, respectively.

## Sealing properties of the biohybrid nanopore

We studied the effect of the length of the alkyl group
$-\textrm{C}_{n}\textrm{H}_{2n+1}$, which functionalizes the interior surface
of the silica pore, on the sealing property of the biohybrid membrane. In
particular, we considered methyl ($n=1$), butyl ($n=2$), and propyl ($n=3$)
groups. In Fig. S2 is shown the schematic representation of the system after
it relaxed for $1$ ns, when the interior surface of the pore was
functionalized with the butyl group. We observed the presence of water
molecules at the interface between the lipid molecules and the interior
surface of the silica pore, which can be explained by the abscence of
preferential interaction between the hydrophobic tail of the lipids and the
butyl groups.

Figure S2: Schematic representation of the Aqp-incorporated lipid shell
inserted in the silica nanopore functionalized with the butyl group, after the
system relaxed for $1$ ns. Red, blue, white, gold, and cyan spheres represent
oxygen, nitrogen, hydrogen, phosphate and carbon atoms in the lipid molecules,
respectively. The oxygen atoms of the water molecules are shown in blue. The
lipid molecules are shadowed and the silica nanopore and the AqpZ monomer are
not shown for clarity. In addition to the water channel which connects both
sides of the membrane, we observe the presence of water molecules at the
interface between the lipid molecules and the interior surface of the silica
pore.

## Density distribution of ions in sea water environment

In Fig. S3 is shown the density distribution of ions (sodium and chloride
ions) in the simulation box after the system relaxed for $200$ ns at constant
temperature $T=300$ K. The ion density in the Aqp-incorporated lipid shell,
i.e. $|z|\leq 1$, is exactly zero, indicating the efficiency of the membrane
to play the role of separator. For $1\leq|z|\leq 4$ we observe the progressive
increase of the ion density, similar to the one observed for the water
molecules (cf. main text). This profile is explained by the conformational
organisation of the lipid bilayer in the nanopore due to preferential
interaction with the hydrophobic surface. For $|z|\geq 4$, i.e. outside the
solid-state nanopore, the system transitions abruptly to reach the bulk
desired ion density.

Figure S3: (a) Schematic representation of the system in sea saline water. The
solid-state nanopore is not shown for clarity. The chloride and sodium ions
are shown in green and yellow, respectively. Red, blue, white, gold, and cyan
spheres represent oxygen, nitrogen, hydrogen, phosphate and carbon atoms in
the lipid molecules, respectively. The oxygen atoms of the water molecules are
shown in blue. (b) Density profile of ions (sodium and chloride ions) along
the Z direction of the simulation box centered around the ApqZ monomer,
measured after the system relaxed for $200$ ns in sea water environment.
Uncertainties, defined as the standard errors, are represented by the shaded
area.

## Equilibrium dissociation of the AqpZ monomer from the MSP-based lipid
nanodisc

The phenomenon of interest (i.e. the dissociation of the AqpZ monomer from the
MSP-based lipid nanodisc) occurs on time scales that are orders of magnitude
longer than those accessible with classical MD simulations. To overcome this
limitation, we use the numerical method based on constrained MD employed by
Sicard and coworkers, F.Sicard et al. (2018); Sicard et al. (2019) which
combined the adiabatic biased molecular dynamics (ABMD) Paci and Karplus
(1999); Marchi and Ballone (1999); Camilloni et al. (2011); Sicard and Striolo
(2016) and umbrella sampling (US) Kastner (2011) frameworks. To design the US
windows, we use the projection of the Cartesian distance between the center of
mass of the AqpZ monomer and the center of mass of the MSPs along the
$Z$-direction, $\textrm{d}_{Z}$, as a reaction coordinate (RC). The starting
configurations for the US simulations are obtained by pulling adiabatically
the system along the RC, generating $180$ windows. Each US window is
subsequently run for $2$ ns to allow equilibration, followed by additional $2$
ns of the production run. Upon completion of the US simulations, we obtain the
free energy profile (FEP) associated with the dissociation of the AqpZ monomer
using the dynamic histogram analysis method (DHAM). Rosta and Hummer (2015)
This unbiasing method uses a maximum likelihood estimate of the Markov State
Model (MSM) transition probabilities given the observed transition counts
during each biased trajectory. To produce MSMs from the biased simulations, we
discretize the RC into bins and count the number of transitions between each
pair of bins $i$ and $j$ in the US simulation $k$ at the chosen lagtime $t$,
$C_{ij}^{k}(t)$, as well as the number of times each bin is occupied during
each US simulation $k$, $n_{i}^{k}=\sum_{j}C_{ij}^{k}(t)$. These values
provide the conditional probabilities $M_{ij}(t)=P(j,t|i,0)$. For biased
simulations, where a biasing energy of $u_{i}^{k}$ is applied to state $i$
during simulation $k$, we compute the unbiased MSM from the biased data, as
given by

$M_{ij}(t)=\frac{\sum_{k}C_{ij}^{k}(t)}{\sum_{k}n_{i}^{k}~{}\exp\big{(}-(u_{j}^{k}-u_{i}^{k})/2k_{B}T\big{)}}.$

Once the MSM is constructed from the simulation data, the equilibrium
probabilities are calculated as the eigenvector corresponding to eigenvalue
$1$ of the transition matrix $M_{ij}$. In Fig. S4 is shown the FEP associated
with the dissocition of the AqpZ monomer plotted along the projection of the
Cartesian distance between the center of mass of AqpZ and the center of mass
of the MSPs. Uncertainties were determined by dividing the data into three
equal sections, determining the profiles independently, and calculating the
standard error. We measured a difference in FE between the global minimum
($\textrm{d}_{Z}\approx 0.1$) and the region where the FEP flattens
($\textrm{d}_{Z}\approx 1$), $\Delta F_{0}\approx 75~{}\textrm{kJ/mol}$.

Figure S4: FEP associated with the equilibrium dissociation of the AqpZ
monomer from the MSP-based lipid nanodisc, obtained within the US/ABMD
framework and calculated with DHAM and a number of bin of 200. The $x$-axis
corresponds to the projection of the Cartesian distance between the center of
mass of AqpZ and the center of mass of the MSPs, $\textrm{d}_{Z}$. The
difference in FE between the global minima ($\textrm{d}_{Z}\approx 0.1$) and
the region where the FEP flattens ($\textrm{d}_{Z}\approx 1$) is $\Delta
F_{0}\approx 75~{}\textrm{kJ/mol}$. Uncertainties, defined as the standard
errors, are represented by the shaded area for US data.