Ali Ziat

Edouard Delasalles

Ludovic Denoyer

Patrick Gallinari

Spatio-Temporal Neural Networks for Space-Time Series Forecasting and Relations Discovery

come   L'archive ouverte pluridisciplinaire

I. INTRODUCTION

Time series exhibiting spatial dependencies are present in many domains including ecology, meteorology, biology, medicine, economics, traffic, and vision. The observations can come from multiple sources e.g. GPS, satellite imagery, video cameras, etc. Two main difficulties when modeling spatiotemporal data come from their size -sensors can cover very large space and temporal lags -and from the complexity of the data generation process. Reducing the dimensionality and uncovering the underlying data generation process naturally leads to consider latent dynamic models. This has been exploited both in statistics [START_REF] Cressie | Statistics for spatio-temporal data, ser. Wiley series in probability and statistics[END_REF] and in machine learning (ML) [START_REF] Bahadori | Fast multivariate spatio-temporal analysis via low rank tensor learning[END_REF], [START_REF] Koppula | Learning spatio-temporal structure from rgbd videos for human activity detection and anticipation[END_REF].

Deep learning has developed a whole range of models for capturing relevant information representations for different tasks and modalities. For dynamic data, recurrent neural networks (RNN) handles complex sequences for tasks like classification, sequence to sequence prediction, sequence generation and many others [START_REF] Bengio | Neural net language models[END_REF]- [START_REF] Li | Gated graph sequence neural networks[END_REF].

These models are able to capture meaningful features of sequential data generation processes, but the spatial structure, essential in many applications, has been seldom considered in Deep Learning. Very recently, convolutional RNNs [START_REF] Shi | Convolutional lstm network: A machine learning approach for precipitation nowcasting[END_REF], [START_REF] Srivastava | Unsupervised learning of video representations using lstms[END_REF] and video pixel networks [START_REF] Kalchbrenner | Unsupervised learning of video representations using lstms[END_REF] have been used to handle both spatiality and temporality, but for video applications only.

We explore a general class of deep spatio-temporal models by focusing on the problem of time series forecasting of spatial processes for different types of data. Although the target in this paper is forecasting, the model can easily be extended to cover related tasks like spatial forecasting (kriging [START_REF] Stein | Interpolation of spatial data: some theory for kriging[END_REF]) or data imputation.

The model, denoted Spatio-Temporal Neural Network (STNN), has been designed to capture the dynamics and correlations in multiple series at the spatial and temporal levels. This is a dynamical system model with two components: one for capturing the spatio-temporal dynamics of the process into latent states, and one for decoding these latent sates into actual series observations.

The model is tested and compared to state of the art alternatives, including recent RNN approaches, on spatiotemporal series forecasting problems: disease prediction, traffic forecasting, meteorology and oceanography. Besides a comparative evaluation on forecasting tasks, the ability of the model to discover relevant spatial relations between series is also analyzed.

The paper is organized as follow: in section II we introduce the related work in machine learning and spatio-temporal statistics. The model is presented in section III with its different variants. The experiments are described in section V for both forecasting V-B and relations discovery V-C.

II. RELATED WORK

The classical topic of time series modeling and forecasting has given rise to an extensive literature, both in statistics and machine learning. In statistics, classical linear models are based on auto-regressive and moving average components. Most assume linear and stationary time dependencies with a noise component [START_REF] De Gooijer | 25 years of time series forecasting[END_REF]. In machine learning, non linear extensions of these models based on neural networks were proposed as early as the nineties, opening the way to many other non linear models developed both in statistics and ML, like kernel methods [START_REF] Muller | Using support vector machines for time series prediction[END_REF] for instance.

Dynamical state space models, such as recurrent neural networks, have been used for time series forecasting in different contexts since the early nineties [START_REF] Connor | Recurrent neural networks and robust time series prediction[END_REF]. Recently, these models have witnessed important successes in different areas of sequence modeling problems, leading to breakthrough in domains like speech [START_REF] Graves | Speech recognition with deep recurrent neural networks[END_REF], language generation [START_REF] Sutskever | Generating text with recurrent neural networks[END_REF], translation [START_REF] Cho | Learning phrase representations using rnn encoder-decoder for statistical machine translation[END_REF] and many others. A model closely related to our work is the dynamic factor graph model [START_REF] Mirowski | Dynamic factor graphs for time series modeling[END_REF] designed for multiple series modeling. Like ours, it is a generative model with a latent component that captures the temporal dynamics and a decoder for predicting the series. However, spatial dependencies are not considered in this model, and the learning and inference algorithms are different.

Recently, the development of non parametric generative models has become a very popular research direction in Deep Learning, leading to different families of innovative and promising models. For example, the Stochastic Gradient Variational Bayes algorithm (SGVB) [START_REF] Kingma | Auto-encoding variational bayes[END_REF] provides a framework for learning stochastic latent variables with deep neural networks, and has recently been used by some authors to model time series [START_REF] Bayer | Learning stochastic recurrent networks[END_REF]- [START_REF] Krishnan | Deep kalman filters[END_REF]. In our context, which requires to model explicitly both spatial and temporal dependencies between multiple time series, variational inference as proposed by such models is still intractable, especially when the number of series grows, as in our experiments.

Spatio-temporal statistics have already a long history [START_REF] Cressie | Statistics for spatio-temporal data, ser. Wiley series in probability and statistics[END_REF], [START_REF] Wikle | A general science-based framework for dynamical spatio-temporal models[END_REF]. The traditional methods rely on a descriptive approach using the first and second-order moments of the process for modeling the spatio-temporal dependencies. More recently, dynamical state space models, where the current state is conditioned on the past have been explored [START_REF] Wikle | Modern perspectives on statistics for spatio-temporal data[END_REF]. For these models, time and space can be either continuous or discrete. The usual way is to consider discrete time, leading to the modeling of time series of spatial processes. When space is continuous, models are generally expressed by linear integrodifference equations, which is out of the scope of our work. With discrete time and space, models come down to general vectorial autoregressive formulations. These models face a curse of dimensionality in the case of a large number of sources. Different strategies have been adopted to solve this problem, such as embedding or parameter reduction. This leads to model families that are close to the ones used in machine learning for modeling dynamical phenomena, and incorporate a spatial components. An interesting feature of these approaches is the incorporation of prior knowledge inspired from physical models of space-time processes. This consists in taking inspiration from prior background of physical phenomena, e.g. diffusion laws in physics, and using this knowledge as guidelines for designing dependencies in statistical models. In climatology, models taking into account both temporal and geographical components have also been used such as Gaussian Markov Random Fields [START_REF] Rue | Gaussian Markov random fields: theory and applications[END_REF] In the machine learning domain, spatio-temporal modeling has been seldom considered, even though some spatio-temporal models have been proposed [START_REF] Ceci | Predictive modeling of pv energy production: How to set up the learning task for a better prediction?[END_REF]. [START_REF] Bahadori | Fast multivariate spatio-temporal analysis via low rank tensor learning[END_REF] introduce a tensor model for kriging and forecasting. [START_REF] Koppula | Learning spatio-temporal structure from rgbd videos for human activity detection and anticipation[END_REF] use conditional random fields for detecting activity in video, where time is discretized at the frame level and one of the tasks is the prediction of future activity. Brain Computer Interface (BCI) is another domain for spatio-temporal data analysis with some work focused on learning spatio-temporal filters [START_REF] Dornhege | Optimizing spatio-temporal filters for improving braincomputer interfacing[END_REF], [START_REF] Ren | Convolutional deep belief networks for feature extraction of eeg signal[END_REF], but this is a very specific and different topic.

III. THE STNN MODEL

A. Notations and Task

Let us consider a set of n temporal series, m is the dimensionality of each series and T their length 1 . m = 1 means that we consider n univariate series, while m > 1 correspond to n multivariate series each with m components. We will denote X the values of all the series between time 1 and time T . X is then a R T ×n×m 3-dimensional tensor, such that X t,i,j is the value of the j-th component of series i at time t. X t will denote a slice of X at time t such that X t ∈ R n×m denotes the values of all the series at time t.

For simplicity, we first present our model in a monorelational setting. An extension to multi-relational series where different relations between series are observed is described in section III-D. We consider that the spatial organization of the sources is captured through a matrix W ∈ R n×n . Ideally, W would indicate the mutual influence between sources, given as a prior information. In practice, it might be a proximity or similarity matrix between the sources: for geo-spatial problems, this might correspond to the inverse of a physical distance -e.g. geodesic -between sources. For other applications, this might be provided through local connections between sources using a graph structure (e.g. adjacency matrix for connected roads in a traffic prediction application or graph kernel on the web). In a first step, we make the hypothesis that W is provided as a prior on the spatial relations between the series. An extension where weights on these relations are learned is presented in section IV.

We consider in the following the problem of spatial time series forecasting i.e predicting the future of the series, knowing their past. We want to learn a model f : R T ×n×m × R n×n → R τ ×n×m able to predict the future at τ time-steps of the series based on X and on their spatial dependency.

B. Modeling Time Series with Continuous Latent Factors

Let us first introduce the model in the simpler case of multiple time series prediction, without considering spatial relations. The model has two components.

The first one captures the dynamic of the process and is expressed in a latent space. Let Z t be the latent representation, or latent factors, of the series at time t. The dynamical component writes Z t+1 = g(Z t ). The second component is a decoder which maps latent factors Z t onto a prediction of the actual series values at t: Xt = d(Z t ), Xt being the prediction computed at time t. In this model, both the representations Z t and the parameters of the dynamical and decoder components are learned. Note that this model is different from the classical RNN formulations [START_REF] Cho | Learning phrase representations using rnn encoder-decoder for statistical machine translation[END_REF], [START_REF] Hochreiter | Long short-term memory[END_REF]. The state space component of a RNN with self loops on the hidden cells writes Z t+1 = g(Z t , X t ), where X t is the ground truth X t during training, and the predicted value Xt during inference. In our approach, latent factors Z t are learned during training and are not an explicit function of past inputs as in RNNs: the dynamics of the series are then captured entirely in the latent space.

This formal definition makes the model more flexible than RNNs since not only the dynamic transition function g(.), but also the state representations Z t are learned from data. A similar argument is developed in [START_REF] Mirowski | Dynamic factor graphs for time series modeling[END_REF]. It is similar in spirit to Hidden Markov models or Kalman filters.

a) Learning problem: Our objective is to learn the two mapping functions d and g together with the latent factors Z t , directly from the observed series. We formalize this learning problem with a bi-objective loss function that captures the dynamics of the series in the latent space and the mapping from this latent space to the observations. Let L(g, d, Z) be this objective function:

L(d, g, Z) = 1 T t ∆(d(Z t ), X t ) (i) +λ 1 T T -1 t=1 ||Z t+1 -g(Z t )|| 2 (ii) (1)
The first term (i) measures the ability of the model to reconstruct the observed values X t from the latent factor Z t . It is based on loss function ∆ which measures the discrepancy between predictions d(Z t ) and ground truth X t . The second term (ii) aims at capturing the dynamicity of the series in the latent space. This term forces the system to learn latent factors Z t+1 that are as close as possible to g(Z t ). Note that in the ideal case, the model converges to a solution where Z t+1 = g(Z t ), which is the classical assumption made when using RNNs. The hyper-parameter λ is used here to balance this constraint and is fixed by cross-validation. The solution d * , g * , Z * to this problem is computed by minimizing L(d, g, Z):

d * , g * , Z * = arg min d,g,Z L(d, g, Z) (2)
b) Learning algorithm: In our setting, functions d and g, described in the next section, are differentiable parametric functions. Hence, the learning problem can be solved endto-end with Stochastic Gradient Descent (SGD) techniques2 directly from [START_REF] Kalchbrenner | Unsupervised learning of video representations using lstms[END_REF]. At each iteration, a pair (Z t , Z t+1 ) is sampled, and Z t , Z t+1 , g and d are updated according to the gradient of [START_REF] Cressie | Statistics for spatio-temporal data, ser. Wiley series in probability and statistics[END_REF]. Training can also be performed via mini-batch, meaning that for each iteration several pairs (Z t , Z t+1 ) are sampled, instead of a single pair. This results in a high learning speed-up when using GPUs which are the classical configuration for running such methods. c) Inference: Once the model is learned, it can be used to predict future values of the series. The inference method is the following: the latent factors of any future state of the series is computed using the g function, and the corresponding observations is predicted by using d on these factors. Formally, let us denote Zτ the predicted latent factors at time T + τ . The forecasting process computes Zτ by successively applying the g function τ times on the learned vector Z T :

Zτ = g • g • ... • g(Z T ) (3)
and then computes the predicted outputs Xτ :

Xτ = d( Zτ ) (4)

C. Modeling Spatio-Temporal Series

Let us now introduce a spatial component in the model. We consider that each series has its own latent representation at each time step. Z t is thus a n × N matrix such that Z t,i ∈ R N is the latent factor of series i at time t, N being the dimension of the latent space. This is different from approaches like [START_REF] Mirowski | Dynamic factor graphs for time series modeling[END_REF] or RNNs for multiple series prediction, where Z t would be a single vector common to all the series. The decoding and dynamic functions d and g are respectively mapping R n×N to R n×m and R n×N to R n×N .

The spatial information is integrated in the dynamic component of our model through a matrix W ∈ R n×n + . In a first step, we consider that W is provided as prior information on the series' mutual influences. In IV, we remove this restriction, and show how it is possible to learn the weights of the relations, and even the spatial relations themselves, directly from the observed data. The latent representation of any series at time t + 1 depends on its own latent representation at time t (intradependency) and on the representations of the other series at t (inter-dependency). Intra-dependency will be captured through a linear mapping denoted Θ (0) ∈ R N ×N and inter-dependency will be captured by averaging the latent vector representations of the neighboring series using matrix W , and computing a linear combination denoted Θ (1) ∈ R N ×N of this average. Formally, the dynamic model g(Z t ) is designed as follow:

Z t+1 = h(Z t Θ (0) + W Z t Θ (1) ) (5)
Here, h is a non-linear function. In the experiments we set h = tanh but h could also be a more complex parametrized function like a multi-layer perceptron (MLPs) for examplesee section V. The resulting optimization problem over d, Z, Θ (0) and Θ (1) writes:

d * , Z * , Θ (0) * , Θ (1) * = arg min d,Z,Θ (0) ,Θ (1) 1 T t ∆(d(Z t ), X t )+ λ 1 T T -1 t=1 ||Z t+1 -h(Z t Θ (0) + W Z t Θ (1) )|| 2 with Z t ∈ R n×N (6)

D. Modeling different types of relations

The model in section III-C considers that all the spatial relations are of the same type (e.g. source proximity). For many problems, we will have to consider different types of relations. For instance, when sensors correspond to physical locations and the target is some meteorological variable, the relative orientation or position of two sources may imply a different type Each possible relation type is denoted r and is associated to a matrix W (r) ∈ R n×n + . For now, and as before, we consider that the W (r) are provided as prior knowledge. Each type of relation r is associated to a transition matrix Θ (r) . This learned matrix captures the spatio-temporal relationship between the series for this particular type of relation. The model dynamics writes:

Z t+1 = h(Z t Θ (0) + r∈R W (r) Z t Θ (r) ) ( 7
)
where R is the set of all possible types of relations. The learning problem is similar to equation ( 6) with Z t+1 replaced by the expression in [START_REF] Shi | Convolutional lstm network: A machine learning approach for precipitation nowcasting[END_REF]. The corresponding model is illustrated in figure 1. This dynamic model aggregates the latent representations of the series for each type of relation, and then applies Θ (r) on this aggregate. Each Θ (r) is able to capture the dynamics specific to relation (r).

IV. LEARNING THE RELATION WEIGHTS AND CAPTURING SPATIO-TEMPORAL CORRELATIONS

In the previous sections, we made the hypothesis that the spatial relational structure and the strength of influence between series were provided to the model through the W (r) matrices. We introduce below an extension of the model where weights on these relations are learned. This model is denoted STNN-R(efining). We further show that, with a slight modification, this model can be extended to learn both the relations and their weights directly from the data, without any prior. This extension is denoted STNN-D(iscovering).

We will first introduce the STNN-R extension. Let Γ (r) ∈ R n×n be a matrix of weights such that Γ (r) i,j is the strength of the relation between series i and j in the relation r. Let us extend the formulation in Equation ( 7) as follows:

Z t+1 = h(Z t Θ (0) + r∈R (W (r) Γ (r) )Z t Θ (r) ) (8)
where Γ (r) is a matrix to be learned, W (r) is a prior i.e a set of observed relations, and is the element-wise multiplication between two matrices. The learning problem can be now be written as:

d * , Z * , Θ * , Γ * = arg min d,Z,Γ 1 T t ∆(d(Z t ), X t ) + γ|Γ| +λ 1 T T -1 t=1 ||Z t+1 -h( r∈(R) (W (r) Γ (r) )Z t Θ (r) )|| 2 (9)
where |Γ (r) | is a l 1 regularizing term that aims at sparsifying Γ (r) . We thus add an hyper-parameter γ to tune this regularization factor.

If no prior is available, then simply removing the W (r) s from equation ( 8) leads to the following model:

Z t+1 = h(Z t Θ (0) + r∈R Γ (r) Z t Θ (r) ) ( 10
)
where Γ (r) is no more constrained by the prior W (r) so that it will represent both the relational structure and the relation weights. Both models are learned with SGD, in the same way as described in III-B. The only difference is that a gradient step on the Γ (r) s is added.

V. EXPERIMENTS

Experiments are performed on a series of spatio-temporal forecasting problems representative of different domains. We consider predictions within a +5 horizon i.e. given a training series of size T , the evaluation of the quality of the model will be made over T + 1 to T + 5 time steps. The different model hyper-parameters are selected using a time-series crossvalidation procedure called rolling origin as in [START_REF] Taieb | Boosting multi-step autoregressive forecasts[END_REF], [START_REF] Ganeshapillai | Learning connections in financial time series[END_REF]. This protocol makes use of a sliding window of size T : on a series of length T , a window of size T shifted several times in order to create a set of train/test folds. The beginning of the T window is used for training and the remaining for test. The value of T is fixed so that it is large enough to capture the main dynamics of the different series. Each series was re-scaled between 0 and 1.

We performed experiments with the following models: (i) Mean: a simple heuristic which predicts future values of a series with the mean of its observed past values computed on the T training steps of each training fold. (ii) AR: a classical univariate Auto-Regressive model. For each series and each variable of the series, the prediction is a linear function of R past lags of the variable, R being a hyper-parameter tuned on a validation set. (iii) VAR-MLP: a vectorial auto-regressive model where the predicted values of the series at time t + 1 depend on the past values of all the series for a lag of size R. The predictive model is a multi-layer perceptron with one hidden layer. Its performance were uniformly better than a linear VAR. Here again the hidden layer size and the lag R were set by validation (iv) RNN-tanh: a vanilla recurrent neural network with one hidden layer of recurrent units and tanh non-linearities. As for the VAR-MLP, one considers all the series simultaneously, i.e. at time t the RNN receives as input X t-1 the values of all the series at t -1 and predicts X t . A RNN is a dynamical state-space model but its latent state Z t explicitly depends through a functional dependency both on the preceding values of the series X t-1 and on the preceding state Z t-1 . Note that this model has the potential to capture the spatial dependencies since all the series are considered simultaneously, but does not (v) RNN-GRU: same as the RNN-tanh, but recurrent units is replaced with gated recurrent units (GRU) units, which are considered state of the art for many sequence prediction problems today 4 . We have experimented with several architectures, but using more than one layer of GRU units did not improve the performance, so we used 1 layer in all the experiments. (vi) Dynamic Factor Graph (DFG): the model proposed in [START_REF] Mirowski | Dynamic factor graphs for time series modeling[END_REF] is the closest to ours but uses a joint vectorial latent representation for all the series as in the RNNs, and does not explicitly model the spatial relations between series.

(vii) STNN: our model where g is the function described in equation ( 7), h is the tanh function, and d is a linear function.

Note that other architectures for d and g have been tested (e.g. multi-layer perceptrons) without improving the quality of the prediction. The λ value has been set by cross validation.

(viii and ix) STNN-R and STNN-D: For the forecasting experiments the γ value of the L 1 penalty (see equation ( 9)) 4 We also performed tests with LSTM and obtained similar results as GRU.

wer set to 0 since higher value decreased the performance, a phenomena often observed in other models such as L 1regularized SVMs. The influence of γ on the discovered spatial structure is further discussed and illustrated in figure 6.

The complete set of hyper-parameters values for the different models is given in appendix.

A. Datasets

The different forecasting problems and the corresponding datasets are described below. The dataset characteristics are provided in table V-A.

• Disease spread forecasting: The Google Flu dataset contains for 29 countries, about ten years of weekly estimates of influenza activity computed by aggregating Google search queries (see http://www.google.org/flutrends). We extract binary relations between the countries, depending on whether or not they share a border, as a prior W . is set to 1 if and only if source j is located 2 degree at the north of source i (the pixel just above on the satellite image).

• Car Traffic Forecasting: The goal is to predict car traffic on a network of streets or roads. We use the Beijing dataset provided in [START_REF] Yuan | Driving with knowledge from the physical world[END_REF], [START_REF] Yuan | T-drive: driving directions based on taxi trajectories[END_REF] which consists of GPS trajectories for ∼ 10500 taxis during a week, for a total of 17 millions of points corresponding to road segments in Beijing. From this dataset, we extracted the traffic-volume aggregated on a 15 min window for 5,000 road segments.

The objective is to predict the traffic at each segment. We connect two sources if they correspond to road segments with a shared crossroad.

For all the datasets but PST (i.e. Google Flu, GHO, Wind and Bejing), we defined the relational structure using a simple adjacency matrix W . Based on this matrix, we defined K different relations by introducing the powers of this matrix: W (1) = W , W (2) = W ×W , etc. In our setting K took values from 1 to 3 and the optimal value for each dataset has been selected during the validation process.

B. Forecasting Results

A quantitative evaluation of the different models and the baselines, on the different datasets is provided in table II. All the results are average prediction error for T + 1 to T + 5 predictions. The score function used is the Root Mean Squared Error (RMSE). A first observation is that STNN and STNN-R models which make use of prior spatial information significantly outperform all the other models on all the datasets. For example, on the challenging PST dataset, our models increase by 23% the performance of the GRU-RNN baseline. The increase is more important when the number of series is high (geo-spatial and traffic datasets) than when it is small (disease datasets). In these experiments, STNN-D is on par with RNN-GRU. The two models do not use prior information on spatial proximity. STNN makes use of a more compact formulation than RNN-GRU for expressing the series mutual dependency but the results are comparable. Vectorial AR logically improves on mono-variable AR (not shown here) and non linear MLP-VAR improves on linear VAR.

We also provide in table III the score for each of the 25 diseases of the GHO dataset. STNN-R obtained the best performance compared to STNN and STNN-D. It outperforms state-of-the-art methods in 20 over 25 datasets, and is very close to the RNN-GRU model on the 5 remaining diseases where RNN-GRU performs best. It thus shows that our model is able to benefit from the neighbour information in the proximity graph.

Figures 2 and3 illustrate the prediction of STNN-R and RNN-GRU on the meteorology and on the oceanography datasets along with the ground truth. Clearly on these datasets, STNN qualitatively performs much better than RNNs by using explicit spatial information. STNN is able to predict fine details corresponding to local interactions when RNNs produce a much more noisy prediction. These illustrations are representative of the general behavior of the two models.

We also provide the performance of the models at different prediction horizons T + 1, T + 2, ....T + 13 on figure 5 for the Google Flu dataset. Results show that STNN performs better than the other approaches for all prediction horizons and is thus able to better capture longer-term dependencies.

Figure 4 illustrates the RMSE of the STNN-R model when predicting at T +1 on the Google Flu dataset for different values of λ. One can see that the best performance is obtained for an average value of λ: low values corresponding to weak temporal

C. Discovering the Spatial Correlations

In this subsection, we illustrate the ability of STNN to discover relevant spatial correlations on different datasets. Figure 6 illustrates the values of Γ obtained by STNN-D where no structure (e.g. adjacency matrix W ) is provided to the model on the PST dataset. Each pixel corresponds to a particular time series and the figure shows the correlation Γ i,j discovered between each series j with a series i roughly located at the center of the picture. The darker a pixel is, the higher the absolute value of Γ i,j is (note that black pixels correspond to countries and not sea). Different levels of sparsity are illustrated from low (up) to high (down). Even if the model does not have any knowledge about the spatial organization of the series (no W matrix provided), it is able to re-discover this spatial organization by detecting strong correlations between close series, and low ones for distant series.

Figure 7 illustrates the correlations discovered on the PST dataset. We used as priors 8 types of relations corresponding to the 8 cardinal directions (South, South-West, etc...). In this case, STNN-R learns weights (i.e Γ (r) ) for each relation based on the prior structure. For each series, we plot the direction with the highest learned weight. The strongest direction for each series is illustrated by a specific color in the figure. For instance, a dark blue pixel indicates that the stronger spatial correlation learned for the corresponding series is the North-West direction. The model extracts automatically relations corresponding to temperature propagation directions in the pacific, providing relevant information about the spatio-temporal dynamics of the system.

The model can be adapted to different situations. Figure 8 represents the temporal evolution of the spatial relations on the PST dataset. For this experiment, we have slightly changed the STNN-R model by making the Γ (r) time dependent according to:

Γ (r) t,j,i = f r (Z i t ) (11)
This means that with this modified model, the spatial relation weights depend on the current latent state of the corresponding series and may evolve with time. In the experiment, f r is a logistic function. On figure 8, the different plots correspond to successive time steps. The color represent the actual sea surface temperatures, and the arrows represent the direction of the stronger relation weights Γ (r) t among the eight possible directions (N, NE, etc). One can see that the model captures coherent dynamic spatial correlations such as global currents directions or rotating motions that gradually evolve with time.

VI. CONCLUSION

We proposed a new latent model for addressing multivariate spatio-temporal time series forecasting problems. For this model the dynamics are captured in a latent space and the prediction makes use of a decoder mechanism. Extensive experiments on datasets representative of different domains show that this model is able to capture spatial and temporal dynamics, and performs better than state of the art competing models. This model is amenable to different variants concerning the formulation of spatial and temporal dependencies between the sources.

For the applications, we have concentrated on forecasting (time based prediction). The same model could be used for interpolating (space based prediction or kriging) or for data imputation when dealing with time-series with missing values.

Fig. 1 :

 1 Fig. 1: Architecture of the STNN model as described in Section III-D .

Fig. 2 :Fig. 3 :

 23 Fig. 2: Prediction of wind speed over around 500 stations on the US territory. prediction is shown at time-step T + 1 for RNN-GRU (centre) and STNN-R (right).

•Fig. 4 :Fig. 5 :

 45 Fig. 4: RMSE on Google Flu w.r.t λ

Fig. 6 :

 6 Fig. 6: Illustrations of correlations Γ discovered by the STNN-D model, with γ in {0.01, 0.1, 1} ( from top to bottom).

Fig. 7 :

 7 Fig. 7: Spatial correlations extracted by the STNN-R model on the PST dataset. The color of each pixel correspond to the principal relation extracted by the model.

Fig. 8 :

 8 Fig. 8: Dynamic spatio-temporal relations extract from the PST dataset on the training set. The color represents the actual sea surface temperature. The arrows represent the extracted spatial relations that evolve through time.

TABLE I :

 I Datasets statistics. n is the number of series, m is the dimension of each series, timestep corresponds to the duration of one time-step and #folds corresponds to the number of temporal folds used for validation. For each fold, evaluation has been made on the next 5 values at T + 1, T + 2, ..., T + 5. The relation columns specifies the number of different relation types used in the experiments i.e the number of W (r) matrices used in each dataset. 1 to 3 means that the best among 1 to 3 relations was selected using cross validation of dependency between the sources. In the experimental section, we consider problems with relations based on the relative position of sources, north, south, west, east, .... The multirelational framework generalizes the previous formulation of the model, and allows us to incorporate more abstract relations, like different measures of proximity or similarity between sources. For instance, when sources are spatiality organized in a graph, it is possible to define different graph kernels, each one of them modeling a specific similarity. The following multi-relational formulation is based on adjacency matrices, and can directly incorporate different graph kernels.

TABLE II :

 II Average RMSE for the different datasets computed for T+1, T+2,...,T+5. Standard deviation was computed by re-training the models on different seeds.

	Disease / Model	AR	VAR-MLP	RNN-GRU Mean	DFG	STNN-R
	All causes	0.237	0.228	0.199	0.35	0.291	0.197
	Tuberculosis	0.407	0.418	0.37	0.395	0.421	0.377
	Congenital syphilis	0.432	0.443	0.417	0.459	0.422	0.409
	Diphtheria	0.406	0.396	0.387	0.404	0.419	0.385
	Malignant neoplasm of esophagus	0.355	0.341	0.341	0.363	0.372	0.345
	Malignant neoplasm of stomach	0.44	0.434	0.431	0.455	0.452	0.43
		0.267	0.254	0.282	0.303	0.301	0.253
	Malignant neoplasm of intestine	0.281	0.29	0.278	0.314	0.305	0.275
	Malignant neoplasm of rectum	0.501	0.499	0.481	0.504	0.509	0.498
	Malignant neoplasm of larynx	0.321	0.313	0.32	0.314	0.329	0.310
	Malignant neoplasm of breast	0.375	0.375	0.382	0.394	0.38	0.36
	Malignant neoplasm of prostate	0.111	0.113	0.109	0.184	0.138	0.109
	Malignant neoplasm of skin	0.253	0.243	0.227	0.264	0.256	0.221
	Malignant neoplasm of bones	0.103	0.099	0.097	0.204	0.173	0.08
	Malignant neoplasm of all other and	0.145	0.157	0.147	0.164	0.169	0.156
	unspecified sites
	Lymphosarcoma	0.15	0.132	0.13	0.231	0.135	0.122
	Benign neoplasms	0.366	0.362	0.332	0.398	0.331	0.331
	Avitaminonsis	0.492	0.474	0.449	0.571	0.58	0.414
	Allergic disorders	0.208	0.217	0.221	0.342	0.24	0.202
	Multiple sclerosis	0.061	0.057	0.061	0.242	0.152	0.056
	Rheumatic fever	0.325	0.31	0.287	0.345	0.313	0.256
	Diseases of arteries	0.302	0.301	0.269	0.345	0.328	0.238
	Influenza	0.141	0.141	0.155	0.23	0.217	0.125
	Pneumonia	0.119	0.128	0.1	0.187	0.187	0.1
	Pleurisy	0.246	0.246	0.247	0.29	0.272	0.245
	Gastro-enteritis	0.386	0.369	0.291	0.394	0.398	0.295
	Disease of teeth	0.344	0.312	0.305	0.413	0.361	0.302

TABLE III :

 III RMSE of STNN-R over the 25 datasets in GHO for T+1, T+2,...,T+5

We assume that all the series have the same dimensionality and length. This is often the case for spatio-temporal problems otherwise this restriction can be easily removed.

In the experiments, we used the Nesterov's Accelerated Gradient (NAG) method[START_REF] Sutskever | On the importance of initialization and momentum in deep learning[END_REF].

VII. ACKNOWLEDGMENTS

Locust project ANR-15-CE23-0027-01, funded by Agence Nationale de la Recherche.

APPENDIX

Hyper-parameters selection

We list in this section the set of hyper-parameters that have been chosen by cross-validation.

• For RNN-models:

-Size of the hidden state = (20, 50, 80, 150, 300, 500)

• For AR:

-Number of lags R ∈ {1, 2, 5, 10, 15, 25}


PhD Pritam Das
email: pritamdas005@gmail.com

PhD Andries D Van Der Meer
email: andries.vandermeer@utwente.nl

Aisen Vivas Msc

Yusuf B Arik Msc

PhD Jean-Christophe Remigy
email: remigy@chimie.ups-tlse.fr

PhD Jean-François Lahitte

PhD Rob G H Lammertink
email: r.g.h.lammertink@utwente.nl

PhD Patrice Bacchin
email: bacchin@chimie.ups-tlse.fr

Tunable microstructured membranes in organ-on-chip to monitor trans-endothelial hydraulic resistance

Tissue engineering is an interdisciplinary field, where scientists from different backgrounds collaborate to address the challenge of replacing damaged tissues and organs via the in vitro fabrication of functional and transplantable biological structures. Because the development and optimization of tissue engineering strategies relies on the complex interaction of cells, materials and the physical-chemical tissue microenvironment, there is a need for experimental models that allow controlled studies of these aspects. Organs-onchips have recently emerged as in vitro models that capture the complexity of human tissues in a controlled manner, while including functional read-outs related to human organ physiology. Organs-on-chips consist of multiple microfluidic cell culture compartments, which are interfaced by porous membranes or hydrogels in which human cells can be cultured, thereby providing a controlled culture environment that resembles the microenvironment of a certain organ including mechanical, biochemical and geometrical aspects. Because organs-on-chips provide both a well-controlled microenvironment and functional read-outs, they provide a unique opportunity to incorporate, evaluate and optimize materials for tissue engineering. Here we introduce a polymeric blend membrane with a three-dimensional double-porous morphology prepared from a poly(ε-caprolactone)/chitosan blends (PCL/CHT) by a modified liquid induced phase inversion technique. The membranes have different physico-chemical, microstructural and morphological properties depending on different PCL/CHT ratios. Big surface pores (macrovoids) provide a suitable microenvironment for the incorporation of cells or growth factors, while an interconnected small porous (macroporous) network allows transfer of essential nutrients, diffusion of oxygen and removal of waste. Human umbilical vein endothelial cells (HUVECs) were seeded on the blend membranes embedded inside an organ-on-chip device. The cellular hydraulic resistance was evaluated by perfusing culture medium at a realistic trans-endothelial pressure of 20 cmH 2 O or 2 kPa at 37 °C after one day and three days post-seeding. By introducing and increasing CHT weight-percentage, the resistance of the cellular barrier after three days was significantly improved. The high tuneability over the membrane physico-chemical and architectural

3

characteristics might potentially allow studies of cell-matrix interaction, cell-transportation and barrier function for optimization of vascular scaffolds using organs-on-chips.

Impact statement (less than 100 words)

Organ-on-chips (OOCs) offer interesting potential for progress in the treatment of diseases and injury in the growing field of tissue engineering and regenerative medicine. The paper presents a new way to develop polymer membrane with a tunable microstructured morphology and to implement this biomaterial inside an OOCs device. The reader should find measurements of the trans-endothelial hydraulic resistance in real time during endothelial cells culture: a simple and controlled way of mimicking human physiological condition for vascular tissue regeneration. This combination of novel biomaterial inside an OOC will explore innovative ideas in tissue engineering field.

Keywords: Transendothelial hydraulic resistance, Transendothelial pressure (TEP), Microstructured 3D porous membrane, Organs-on-chips, Microfluidics 1. Introduction:

To evaluate specific tissue engineered constructs, current in vitro methods are relatively limited, because they either rely on simple monolayer cell culture systems (e.g. for the evaluation of cytotoxicity) or on low-throughput bioreactor systems. Because of this, there is a clear need of more advanced in vitro models that can be used as tools to evaluate and optimize tissue engineered constructs. Organs-on-chips (OOCs) are considered to better predict cellular responses and potentially mimic a particular organ by introducing organlike features such as fluid flow and mechanical stress with the corresponding cellular preferences (1)(2)(3).

Generally, OOCs consist of polydimethylsiloxane (PDMS)-based structures fabricated using soft lithography [START_REF] Xia | Soft Lithography[END_REF][START_REF] Huh | Microengineered physiological biomimicry: organs-on-chips[END_REF][START_REF] Capulli | Approaching the In Vitro Clinical Trial: Engineering Organs on Chips[END_REF] defining microfluidic channels that are often interfaced through a porous membrane material in which human cells can be cultured. Depending on the application, the membrane normally requires micron-sized pores to allow cellular interaction, while thickness and biocompatibility need to closely mimic biological conditions of the human body [START_REF] Huh | From 3D cell culture to organs-on-chips[END_REF][START_REF] Huh | Microfabrication of human organs-on-chips[END_REF]. Moreover, a straightforward way to resemble the microenvironment of a certain organ by tuning mechanical, biochemical and geometrical aspects [START_REF] Meer Ad Van Der | Organs-on-chips: breaking the in vitro impasse[END_REF] leads to better predictions [START_REF] Bhatia | Microfluidic organs-on-chips[END_REF]. Like all in vitro models, OOCs hold the potential of parallelization and high-throughput screening with very limited resources.

In spite of the substantial progress in OOC technologies in general (better sensors, pumps and microfluidics) with a rapidly expanding library of OOCs of different organs, the progress in the development of suitable porous, polymeric membranes is very slow. To date a handful of polymeric membrane such as polycarbonate (PC) and Polyethylene terephthalate (PET), are applied for interfacial support in OOCs. These membranes are generally prepared by track-etching and are mainly obtained from commercial Transwell inserts [START_REF] Van Der Helm | Direct quantification of transendothelial electrical resistance in organs-on-chips[END_REF][START_REF] Abhyankar | A Reversibly Sealed, Easy Access, Modular (SEAM) Microfluidic Architecture to Establish In Vitro Tissue Interfaces. Eddington DT[END_REF][START_REF] Lee | A pumpless multi-organ-on-a-chip (MOC) combined with a pharmacokinetic-pharmacodynamic (PK-PD) model[END_REF][START_REF] Loskill | WATon-a-chip: a physiologically relevant microfluidic system incorporating white adipose tissue[END_REF][START_REF] Shim | Microfluidic gut-on-a-chip with three-dimensional villi structure[END_REF][START_REF] Zhu | A vertical-flow bioreactor array compacts hepatocytes for enhanced polarity and functions[END_REF][START_REF] Pocock | Low-temperature bonding process for the fabrication of hybrid glass-membrane organ-on-a-chip devices[END_REF][START_REF] Maoz | Organs-on-Chips with combined multi-electrode array and transepithelial electrical resistance measurement capabilities[END_REF] or filter membranes [START_REF] Van Der Helm | Direct quantification of transendothelial electrical resistance in organs-on-chips[END_REF][START_REF] Wang | Microfluidic blood-brain barrier model provides in vivo-like barrier properties for drug permeability screening[END_REF](20)[START_REF] Ma | Microdroplet chain array for cell migration assays[END_REF][START_REF] Lee | Construction of 3D multicellular microfluidic chip for an in vitro skin model[END_REF][START_REF] Jie | Integrated microfluidic system for cell co-culture and simulation of drug metabolism[END_REF][START_REF] Brown | Metabolic consequences of inflammatory disruption of the blood-brain barrier in an organ-onchip model of the human neurovascular unit[END_REF]. Depending on the application, they often provide only a basic, non-physiological cell culture microenvironment due to their limited biocompatible, structural and mechanical properties. All the material properties of the membrane like surface microstructure [START_REF] Chen | Geometric control of cell life and death[END_REF][START_REF] Hulsman | Analysis of high-throughput screening reveals the effect of surface topographies on cellular morphology[END_REF][START_REF] Esch | Characterization of in vitro endothelial linings grown within microfluidic channels[END_REF][START_REF] Broaders | Coupling between apical tension and basal adhesion allow epithelia to collectively sense and respond to substrate topography over long distances[END_REF], roughness (29-32), hydrophilicity [START_REF] Lee | Interaction of Different Types of Cells on Polymer Surfaces with Wettability Gradient[END_REF][START_REF] Wala | Influence of physico-mechanical properties of elastomeric material for different cell growth[END_REF], 5 mechanical strength [START_REF] Wala | Influence of physico-mechanical properties of elastomeric material for different cell growth[END_REF][START_REF] Engler | Matrix elasticity directs stem cell lineage specification[END_REF], and porosity [START_REF] Casillo | Membrane Pore Spacing Can Modulate Endothelial Cell-Substrate and Cell-Cell Interactions[END_REF][START_REF] Choi | Three-dimensional Scaffolds for Tissue Engineering: The Importance of Uniformity in Pore Size and Structure[END_REF][START_REF] Luca | Gradients in pore size enhance the osteogenic differentiation of human mesenchymal stromal cells in three-dimensional scaffolds[END_REF][START_REF] Loh | Three-Dimensional Scaffolds for Tissue Engineering Applications: Role of Porosity and Pore Size[END_REF][START_REF] O'brien | The effect of pore size on cell adhesion in collagen-GAG scaffolds[END_REF][START_REF] Zhang | Pore size effect of collagen scaffolds on cartilage regeneration[END_REF] have a significant effect on the cellular adhesion and biocompatibity.

Recently, blends of synthetic and natural polymers such as polycaprolactone (PCL) and chitosan (CHT) were found to exhibit excellent abilities to match the needs of their specific target tissues via proper tuning of their properties [START_REF] Du | Gradient nanofibrous chitosan/poly ɛ-caprolactone scaffolds as extracellular microenvironments for vascular tissue engineering[END_REF][START_REF] Correlo | Properties of melt processed chitosan and aliphatic polyester blends[END_REF][START_REF] Malheiro | New poly(ε-caprolactone)/chitosan blend fibers for tissue engineering applications[END_REF][START_REF] Van Der Schueren | Polycaprolactone/chitosan blend nanofibres electrospun from an acetic acid/formic acid solvent system[END_REF][START_REF] Cooper | Fabrication and cellular compatibility of aligned chitosan-PCL fibers for nerve tissue regeneration[END_REF]. CHT is widely recognized for its natural ability to promote bioactive functionality, tunable degradation, and structural resemblance of native tissue extra cellular matrix (ECM) [START_REF] Teotia | Bifunctional Polysulfone-Chitosan Composite Hollow Fiber Membrane for Bioartificial Liver[END_REF][START_REF] Prabaharan | Preparation and characterization of poly(L-lactic acid)-chitosan hybrid scaffolds with drug release capability[END_REF][START_REF] Chen | Electrospun chitosan-P(LLA-CL) nanofibers for biomimetic extracellular matrix[END_REF]. PCL is approved by the Food and Drug Administration (FDA) for clinical use, and is also very suitable for tissue engineering applications [START_REF] Hong | Fabrication of electrospun polycaprolactone biocomposites reinforced with chitosan for the proliferation of mesenchymal stem cells[END_REF][START_REF] Das | Artificial membranes tuning for lymphatic wall repair[END_REF][START_REF] Li | Towards the development of polycaprolactone based amphiphilic block copolymers: molecular design, self-assembly and biomedical applications[END_REF][START_REF] Morelli | Neuronal growth and differentiation on biodegradable membranes[END_REF][START_REF] Xiao | Preparation and cytocompatibility of chitosan-modified polylactide[END_REF] by providing attractive options for the control of morphology, architecture, dynamic modulus and degradation rate to generate potential alternatives via tuning the synthetic/natural ratio in combination with natural polymers [START_REF] Du | Gradient nanofibrous chitosan/poly ɛ-caprolactone scaffolds as extracellular microenvironments for vascular tissue engineering[END_REF][START_REF] Correlo | Properties of melt processed chitosan and aliphatic polyester blends[END_REF][START_REF] Malheiro | New poly(ε-caprolactone)/chitosan blend fibers for tissue engineering applications[END_REF][START_REF] Van Der Schueren | Polycaprolactone/chitosan blend nanofibres electrospun from an acetic acid/formic acid solvent system[END_REF][START_REF] Cooper | Fabrication and cellular compatibility of aligned chitosan-PCL fibers for nerve tissue regeneration[END_REF].

Here we report a straightforward way to integrate a three-dimensional porous membrane, intended for use in vascular tissue engineering, in an OOC system. We demonstrate that endothelial cells can be embedded in the scaffolds, and their functionality can be probed by measuring their hydraulic resistance with a microfluidic assay. Our study highlights the great potential of applying of OOCs as models for in vitro testing in the field of tissue engineering.

Materials and Methods:

Organ-on-chip fabrication

The organ-on-chip consists of a porous membrane embedded inside a microfluidic device.

The chip fabrication was reported in the previous publication from our group [START_REF] Van Der Helm | Direct quantification of transendothelial electrical resistance in organs-on-chips[END_REF][START_REF] Helm Mw Van Der | Fabrication and Validation of an Organ-on-chip System with Integrated Electrodes to Directly Quantify Transendothelial Electrical Resistance[END_REF]. The membrane is made from a biodegradable PCL/CHT polymer blend with a modified phase inversion technique to allow the formation of double porous morphology. The methods for the chip and the membrane fabrication are briefly described in Supplementary Information (SI 2 and SI 3). The chip design and assembled chip are shown in 3D (Blender v2.78) in Fig.

(A)

Hydraulic resistance measurements in different pressure by Fluidic set up

To evaluate the membrane hydraulic performance, volume flow rate of the culture medium (EGM-2) through the membrane at different pressures was conducted by a fluidic setup (Fig. 1. (B)). The fluidic setup (Fluigent) controls the pressure and tracks the flow rate in real time [START_REF]FRP Flow Sensor: FLOW-RATE PLATFORM | Fluigent 2018[END_REF]. The experiments were started from 2 kPa (approximately 20 cmH 2 O: range of physiological transendothelial pressures or TEP [START_REF] Sill | Shear stress increases hydraulic conductivity of cultured endothelial monolayers[END_REF][START_REF] Pang | Shear stress regulates HUVEC hydraulic conductivity by occludin phosphorylation[END_REF][START_REF] Luckett | Hydraulic conductivity of endothelial cell monolayers cultured on human amnion[END_REF]) and continued up to 20 kPa. The experimental set up was discussed in section 2.3.1. The obtained flux was plotted against pressure (Fig. 3), and the hydraulic resistance of the membranes was determined (Table . 2) from the slope by equations below [START_REF] Desmond | Linking composition of extracellular polymeric substances (EPS) to the physical structure and hydraulic resistance of membrane biofilms[END_REF][START_REF] Desmond | Stratification in the physical structure and cohesion of membrane biofilms -Implications for hydraulic resistance[END_REF][START_REF] Tripathi | Ultralow fouling membranes by surface modification with functional polydopamine[END_REF].

𝐽 = (𝑄/𝑡)/𝐴 (1)
And, 𝐽 = ∆ µ ×

(2)

Or, 𝑅 = ∆ × µ  (3)
Where, J is Flux (m/s), (Q/t) is the volume flow rate (m 3 /s), A is membrane surface area (m 2 ), µ is 6.913 × 10 -4 Pa•s, the dynamic viscosity of water at 37 °C (all the flow experiments were done with EGM-2 but its viscosity is only 5% higher [START_REF] Sill | Shear stress increases hydraulic conductivity of cultured endothelial monolayers[END_REF][START_REF] Fröhlich | Comparison of two in vitro systems to assess cellular effects of nanoparticles-containing aerosols[END_REF], while the resistance of the membrane should not be affected by this), ∆ applied pressure difference (Pa) and R Mem is resistance (m -1 ) of the membrane.

Flux measurements before and after the cell seeding by Flow sensing set up

The consecutive experiments were done on the same batch of chips with culture medium (EGM-2) at 2 kPa for one day and three days (Fig. 4). For each batch, three chips were taken.

Flux before cell seeding

Chips were transferred in a plastic culture dish and washed with 70% EtOH thoroughly.

Then the chips were washed with sterile PBS three times to ensure removal of EtOH from inside of the channel and both the channels were filled with endothelial cell culture medium (EGM-2). One top and one bottom entrance was blocked with sealed yellow 7 pipette tips. The other top channel inlet was connected with the reservoir (30 ml of culture medium) via Tygon tubing. The other bottom channel inlet was connected to the flow rate sensor. Before connecting with the chip, all the tubing and the sensor were filled with culture medium. Then the flow was started at 2 kPa pressure and the flow rate was recorded after 2-3 minutes to ensure a steady-state measurement.

Cell seeding

After the flux measurements the channels were washed three times with sterile PBS and kept in PBS suspension for overnight at 4 °C. Then the channels were washed again with sterile PBS and filled with 0.1 mg/ml collagen I in PBS. The channels were inspected for air bubbles (any air bubbles were flushed out by Collagen I solution) and incubated for 30 minutes at 37 °C. Then the chips were taken out from the incubator and flushed with culture medium and incubated for 1 hour at 37 °C. In the meantime, HUVECs were trypsinized and collected from a confluent T175 flask. Then the chips were placed in the sterile hood and 10 µl of cell suspension (4 × 10 6 cells/ml or 1.5 × 10 5 cells/cm 2 ) was pipetted slowly inside the top channel and inspected for any bubble formation (any bubble was flushed out by the cell suspension) and incubated in 37 °C for 2 hours. The chips were taken out from the incubator, inspected for complete attachment of the HUVECs by microscope, flushed with fresh culture medium to remove non-attached cells and incubated again. The medium in the channels was changed twice every day (morning and afternoon). Every time before putting inside the incubator, one drop of culture medium was left on top of all the inlets to prevent the channel from drying.

Flux after cell seeding and staining

The chips with seeded cells were taken one by one from the incubator and the flow rate was recorded again like the previous procedure. Directly after the perfusion experiment, the chip was washed with PBS and fixed with 4% paraformaldehyde for 15 minutes. Then the channels were washed with PBS and kept in permeabilization buffer (PB; 9:1 v:v ratio of 11 mg/ml BSA with 1% Triton X-100 in PBS) at room temperature for 15 minutes to block any non-specific protein-binding sites. VE-cadherin primary antibody (1:200 in PB) was added inside the channel and incubated at room temperature for 2 hours. The channels were washed with PBS and VE-cadherin secondary antibody, Phallodin-488 and DAPI were added together in the ratio 1:400, 1 drop/ml and 1:4000 in PB, respectively, and incubated at room temperature for 1 hour. To prevent photo-bleaching, the chips were protected from light exposure from this step onwards. Then the chips were washed with PBS thoroughly and imaged with confocal microscopy (Zeiss LSM 510, Nikon) in 15× and 85× magnification (Fig. 5). The images on the membrane were taken at the intersection point of the two channels. For control, cells on the bare PDMS channel surfaces were also imaged.

Calculation of the hydraulic resistance of HUVECs

From the flow rate, flux was calculated from Eq. 1 and plotted for different membranes, before cell seeding and after cell seeding (Table . 3). In addition, the resistance of the cellular barrier (R Cell ) was calculated from the difference between the total measured resistance (membrane with the cells) and the blank membrane resistance (without cells).

The equation has shown below [START_REF] Desmond | Linking composition of extracellular polymeric substances (EPS) to the physical structure and hydraulic resistance of membrane biofilms[END_REF][START_REF] Desmond | Stratification in the physical structure and cohesion of membrane biofilms -Implications for hydraulic resistance[END_REF].

𝑅 = 𝑅 -𝑅 = ∆ µ × 𝐽 -( ∆ µ × 𝐽 ) (4)

Statistical analysis

All the experiments were done in triplicate (n=3) and the data are shown as mean ± standard deviation. Slope of the blank membrane (Table .2) was determined by using LINST equation. To find the significance of flux after cell seeding (Fig. 4), one-way ANOVA analysis with Bonferroni correction was done where; ****denotes p<0.0001, ***denotes 0.0001<p<0.001, **denotes 0.001<p<0.01 and *denotes 0.01<p<0.05.

Results and discussion:

Membrane morphological characteristics

After the membrane fabrication, the top surface and the cross section characteristics were analyzed by SEM (Fig. 2 and Table 1) according to the method presented in SI 4. It was found that the membrane surface that was covered by the track-etched membrane, contained many open big pores which are interconnected with the small pores network as evident from the cross section images. In case of pure PCL, surface big pores were around 21 µm in diameter and inside the bulk were around 30 µm. On the other hand, all the blend of PCL/CHT, surface big pores were around 14 -16 µm and inside the bulk 51 -56 µm. The size of the small pores was around 0.8 -2 µm. Thickness of all the membranes were around 55 -85 µm. Overall, by this modified phase inversion technique, a remarkable difference was found between the big pores and the small pores. The morphology also contained some pillars, which were almost the same dimension as the track-etched membrane pores. The amount of pillars was decreasing with increasing CHT %.

The commercial track-etched membrane upon the casted polymer solution triggered a heterogeneous solvent exchange during the solvent-nonsolvent exchange. The solvent exchange is faster where the nonsolvent has a direct contact with the solution through the pores of the track-etched membrane. By introducing and increasing the CHT wt% inside the PCL backbone, several material characteristics can be further tuned like surface pore size and morphology, mechanical property, hydrophilicity, biodegradable properties.

Nonetheless, all these characteristics have direct consequences on cellular viability and adhesion phenomena. As stated earlier, the goal of the current work was to study the suitability of these membranes in vascular tissue engineering by evaluating the adhesion of HUVECs, as well as their functionality in terms of generating hydraulic resistance.

Permeate flux and hydraulic resistance of the membrane at different pressures

The OOC was designed in such a way that the membrane scaffolds are embedded in between two channels where the culture medium were perfused through the top channel at realistic TEP (transendothelial pressure; 2 kPa), mimicking the inner lumen properties.

The other end of the top channel was blocked so the medium passes through the membrane porous network to the bottom channel. One exit of the bottom channel was also blocked such that the cell operates in dead-end filtration mode. Firstly, we tested the membrane stand-alone performance, starting from TEP up to 10 times TEP. After that we enhanced the resistance by seeding the endothelial cell and quantifying the additional hydraulic resistance induced by the cells. This kind of membrane scaffolds loaded with endothelial cells can be applied as grafts for disrupted vascular lumen and resist the 10 leakage of fluid from inside of the lumen to outside. Due to the biodegradable porous structure, the cell will grow inside the scaffolds by replacing the scaffolds material and hence repair the disrupted lumen.

Permeate flux and Hydraulic resistance of the membranes without cells

The permeate flux of the membranes were plotted (Fig. 3) against pressure starting from realistic TEP (2 kPa). All the membranes showed a linear response across the entire pressure range. In the blend, by increasing CHT wt%, the flux was increasing.

The hydraulic resistance of the membrane was obtained from the slopes in Fig. 3 and depicted in Table 2.

The pristine membrane does not exhibit any membrane compressibility in the TEP range (2 -16 kPa) [START_REF] Sill | Shear stress increases hydraulic conductivity of cultured endothelial monolayers[END_REF][START_REF] Pang | Shear stress regulates HUVEC hydraulic conductivity by occludin phosphorylation[END_REF][START_REF] Luckett | Hydraulic conductivity of endothelial cell monolayers cultured on human amnion[END_REF][START_REF] Madhavan | Mechanical and biocompatible characterizations of a readily available multilayer vascular graft[END_REF], suggesting in the real-time the bulk microstructure of the membrane will be intact while the cellular ingrowth should not be affected. The blend membrane has a higher resistance (or lower flux) than the control (Polyester membrane), most probably because the blend membranes have a higher hydraulic length (i.e. the membrane thickness) [START_REF] Bazyar | Liquid-liquid displacement in slippery liquid-infused membranes (SLIMs)[END_REF] and different pore morphology. These blend membranes have higher resistance even without the cellular attachment than the membrane like Polyester or Polycarbonate filter membrane which are typically used inside the trans wells inserts systems [START_REF] Sill | Shear stress increases hydraulic conductivity of cultured endothelial monolayers[END_REF][START_REF] Pang | Shear stress regulates HUVEC hydraulic conductivity by occludin phosphorylation[END_REF] or in OOC system [START_REF] Van Der Helm | Direct quantification of transendothelial electrical resistance in organs-on-chips[END_REF][START_REF] Helm Mw Van Der | Fabrication and Validation of an Organ-on-chip System with Integrated Electrodes to Directly Quantify Transendothelial Electrical Resistance[END_REF] to asses cellular barrier. The water flux of the control (polyester membrane) in the manufacturer's protocol (66) is given as -1000 ml/min/cm 2 at 10 psi pressure, corresponds to a membrane resistance of 6.25 × 10 8 m -1 , which is very close to what we find in our OOC system. This good agreement confirms that the system allows determining the resistance reliably and that the membrane sealing was intact even up to 10 times TEP.

The flux is very important to applications of membranes in vascular biology as it dictates the amount of fluid flow that can be processed for a given time and dimension at a defined applied pressure without damaging the mechanical integrity of the membrane [START_REF] Barhate | Nanofibrous filtering media: Filtration problems and solutions from tiny materials[END_REF].

Particularly in our blend membrane by increasing CHT wt% the flux was increased likely due to higher hydrophilicity and improved wetting. However, regardless of different chemical property and pore morphology, the resistance of all the blend membranes were quite close to each other. Unlike the polyester membrane, the advantage of the blend microporous structure is that they consists of big pores which are not all the way through the thickness, reducing the flux while facilitate cellular adhesion and incorporation in the same time. Starting from TEP, the applied pressure range (2 -20 kPa) was well within the range of transmembrane pressure for microfiltration membrane [START_REF] Bjorge | Performance assessment of electrospun nanofibers for filter applications[END_REF]. Copper et. al [START_REF] Cooper | Chitosan-based nanofibrous membranes for antibacterial filter applications[END_REF] evaluated the water flux of a PCL-CHT based nanofibrous membranes by electrospining.

They found their membranes exhibited flux values in the same order of magnitude with our reported data in the similar pressure range. Madhavan et. al [START_REF] Madhavan | Mechanical and biocompatible characterizations of a readily available multilayer vascular graft[END_REF] reported vascular grafts made of thin, dense nanofibrous core composed of poly-ε-caprolactone, and a thick, porous, hydrogel sleeve composed of genipin-crosslinked collagen-chitosan. They achieved a desired flux 528 ml/cm 2 /min (corresponds to 0.08 m/s) at average blood pressure of 16 kPa and mentioned flux value in that range considered as ideal for immediate grafting and fast healing. They also mentioned flux value below that range does not require any preclotting before implantation, but healing could be slower due to limited molecule transport. Whereas unlike the multilayer morphology they reported, the single layer morphology with macroporous structure like us could be an added advantage with optimum flux property and the fast healing.

Permeate flux and Hydraulic resistance of membranes with cell cultures

Flux was measured after one and three days of cell culturing and compared with the blank membranes at physiological transendothelial pressure (TEP) 2 kPa (Fig. 4). By introducing and increasing the CHT wt% in PCL, the %reduction of flux was depicted in the figure ,   where it can be found there was a significant decrease of flux took place after three days compare to one day, most likely due to added cellular resistance.

The primary point of interest of this work was to evaluate the additional hydraulic resistance induced by the HUVECs during the flux measurements. These values were derived from the mean flux after one day and three days of cell culture from Eq. 4 and listed in Table . 3.

The cellular resistance (R Cell ) increased the total cell-membrane resistance in the OOC system. After one day, by increasing CHT wt%, the cellular resistance of the membrane 12 scaffolds were decreasing perhaps due to enhanced flux property induced by the CHT moieties as discussed above. However, a notable increase in cellular resistance (R Cell ) after three days of seeding was observed presumably due to higher numbers of adherent cells and the increase is likely to be more significant by increasing CHT wt%. After three days %R Cell increase compare to one day was as follows; PCL/CHT 100/0 (+20 %), PCL/CHT 90/10 (+50 %), PCL/CHT 80/20 (+130 %) and PCL/CHT 70/30 (+200 %). This could be due to combination of better morphological property or pore spacing as evident from the SEM images (Fig. 2) where it was observed that, by increasing CHT wt% the membrane surface has more continuous and smooth space between the pores. Casillo et al. [START_REF] Casillo | Membrane Pore Spacing Can Modulate Endothelial Cell-Substrate and Cell-Cell Interactions[END_REF] grew endothelial cells on non-porous and porous SiO 2 membranes. The porous membranes had different pore sizes and corresponding pore spacing, and they found membranes with higher pore spacing give better focal adhesion formation. In contrast, they measured cellcell interactions in the form of tight junctions by ZO-1, which showed an opposite trend.

Porosity especially with smaller pore spacing, could thus limit cell adhesion, but enhance barrier function of cells.

By increasing CHT wt%, the hydrophilicity of the membrane was also increased which could have a positive influence on protein adsorption and hence cell adhesion (29,33).

Wala et al. (34) prepared poly(dimethyl siloxane) (PDMS) films with different hydrophilicity

and reported that cell attachment, spreading and growth of 3T3 fibroblasts and HaCaT keratinocytes increased with increasing hydrophilicity of the films. Others found that moderate hydrophilicity often yields the best cell response [START_REF] Lee | Interaction of Different Types of Cells on Polymer Surfaces with Wettability Gradient[END_REF][START_REF] Premnath | Tuning cell adhesion by direct nanostructuring silicon into cell repulsive/adhesive patterns[END_REF]. Adsorption of serum proteins also showed an optimum on films with moderate hydrophilicity [START_REF] Lee | Interaction of Different Types of Cells on Polymer Surfaces with Wettability Gradient[END_REF] which could be the main cause of the increased cell response. Above all, the big pores inside the scaffold are likely to play a significant role to increase the cellular resistance. During the cell seeding the HUVECs could block the big pores as they have a globular shape of around 15-20 µm in diameter which is almost similar in dimension to the big pores on the surface of the scaffold. Chen et al. [START_REF] Chen | Geometric control of cell life and death[END_REF] showed that pore sizes similar or larger than a single cell have significant influence on clusters formation by involving multiple cells together.

Microscopic observation of cells on the membrane

The cellular adhesion was observed after the flux experiments in the confocal images (Fig. 5). After one day (in 15× Mag.) all four batches of membranes showed good cellular adhesion and compatibility. The cellular adhesion was increasing by increasing CHT wt%, which could be due to an increase of surface hydrophilicity and functional property like enhanced collagen coating. Whereas in 85× magnification images, the cell-cell interactions looks better on PCL/CHT 100/0 and PCL/CHT 90/10 scaffolds by expressing VE-Cadherin stronger than the other two scaffolds. The cellular morphology was the best on PDMS by strongly expressing VE-Cadherin with confluent cell-cell interactions and cortical F-actin, presumably because the cells prefer a smoother surface during the initial attachments [START_REF] Casillo | Membrane Pore Spacing Can Modulate Endothelial Cell-Substrate and Cell-Cell Interactions[END_REF].

After three days, the cellular adhesion was clearly higher for PCL/CHT 80/20 and PCL/CHT 70/30 scaffolds than PCL/CHT 100/0 and PCL/CHT 90/10. On PDMS (control) the cellular attachments decreased compared to our membrane scaffolds. This could be due to detachment of the collagen layer from hydrophobic PDMS which eventually detaches the cellular layer. The detachment of the collagen from the PDMS likely occurred due to shear stress during medium changing followed by flux experiments. Normally to establish a strong bonding between PDMS and the collagen layer, chemical modification by aminosilane surface treatment is used inside the chips before collagen coating (2). In our case, we skipped this chemical surface modification step purposely in order to compare cellular attachment exclusively on the membrane scaffolds and the control (PDMS surface) during long-term cell culture. In 85× magnification, HUVECs on all the scaffolds strongly express junctional VE-cadherin indicating better cell-cell connectivity compared to day one. In addition, the cells were spreading on the surface by increasing CHT wt%, presumably due to its higher hydrophilicity and better pore spacing of the scaffolds which is evident in the confocal images with spreading of actin filament bundles stained by Phalloidin. For comparison, Xiao et. al. [START_REF] Xiao | The promotion of angiogenesis induced by three-dimensional porous beta-tricalcium phosphate scaffold with different interconnection sizes via activation of PI3K/Akt pathways[END_REF] seeded HUVECs in three-dimensional β-tricalcium phosphate scaffolds with big pores connected with different interconnection sizes and studied cytoskeletal proteins such as actin to illustrate the spreading and attachment of cells on scaffolds. They found the scaffolds with higher pore spacing exhibited uniform and 14 homogenous actin distribution of cytoplasm. This phenomenon is understandable, as the big pores mainly provide the space for cell growth, whereas the interconnections function as the doorway for cell-cell interaction and eventually cell growth into the scaffold [START_REF] Xiao | The promotion of angiogenesis induced by three-dimensional porous beta-tricalcium phosphate scaffold with different interconnection sizes via activation of PI3K/Akt pathways[END_REF][START_REF] Mastrogiacomo | Role of scaffold internal structure on in vivo bone formation in macroporous calcium phosphate bioceramics[END_REF]. Du,et al. (42) reported the viability of HUVECs on a gradient nano-fibrous PCL-CHT membrane scaffolds for vascular tissue engineering. They found that the HUVECs have a tendency for high proliferation and adhesion on the luminal side (CHT rich) of the gradient scaffolds due to a better ECM proteoglycan-like structure and adsorption of essential growth factors (VEGF) by the CHT matrix. Our data showed that by increasing CHT wt% in PCL-CHT membranes, the microporous structure was significantly improved promoting higher endothelial cell types growth on the scaffolds. The entire layer of endothelial cell adhesion on the membrane top surface (lumen side) and the scaffold itself as the adventitial side is a potential match to the natural structure of blood vessels.

In the past, several experiments were reported to monitor the hydraulic transportation phenomena across the endothelial cell layer seeded on membrane at TEP range [START_REF] Sill | Shear stress increases hydraulic conductivity of cultured endothelial monolayers[END_REF][START_REF] Pang | Shear stress regulates HUVEC hydraulic conductivity by occludin phosphorylation[END_REF][START_REF] Luckett | Hydraulic conductivity of endothelial cell monolayers cultured on human amnion[END_REF].

Almost all those experiments were conducted in a trans-well insert system with the conventional track-etched filter membranes. The current study is a further step forward in vascular tissue engineering by demonstrating how these three-dimensional membranes are more convenient candidates with tunable physico-chemical and microstructural properties when analyzed within OOCs.

Conclusion

We showed a straightforward and simple process to quantify the additional hydraulic resistance of HUVECs seeded on three-dimensional membrane scaffolds in an OOC system.

The membrane scaffolds were developed by a modified version of liquid induced phase inversion, using different ratios of polycaprolactone-chitosan. At first, the flux and the corresponding resistance of the membranes were evaluated in the OOC system, starting from realistic trans-endothelial pressure (TEP). The flux was found to increase with increasing CHT content in the blend. The additional hydraulic resistance induced by the living HUVECs was quantified in the same set-up seeded on the membranes after one day and three days using EGM-2 culture medium. Results suggest that by increasing the CHT 15 wt%, the microstructural morphology and the chemical properties are significantly improved leading to an increased hydraulic resistance of the cells due to enhanced cellular adhesion.

Overall, our study demonstrates that OOC systems can be used as in vitro platforms to optimize and evaluate the morphology and functionality of various tissue engineering constructs. We anticipate that the use of OOC technology in tissue engineering will increase, because the unique properties of OOCs, with their small size, engineered microenvironments and tissue-level or organ-level functional read-outs. Thus could greatly facilitate the development and optimization of new tissue engineered constructs for regenerative medicine.

Acknowledgment

The author would like to acknowledge Marinke van der Helm and Paul ter Braak (BIOS/Lab on a Chip, UTwente) for kindly providing the inset photo in was formed where the non-solvent was not able to make a direct entry (Fig. 2).

After formation of the membrane, when it turned opaque and detached from the glass plate, the track-etched membrane was removed gently and the newly formed membrane was washed several times in pure water. Then, the membrane was stored in ultrapure water in 4 °C temperature. The surface where the big pores were open towards the nonsolvent is denoted as top surface and the other surface, which was facing the glass plate is denoted as the bottom surface.

Due to the high viscosity of CHT in the solvent, we could not go beyond the ratio PCL/CHT 70/30 as it was quite difficult to make a homogeneous solution. Double porous membrane formation with pure CHT was not successful as dissolving CHT more than 4 wt.% was difficult.

SI 3. Organ-on-chip fabrication

Briefly, PDMS base agent and curing agent were mixed in a 10:1 wt% ratio. After degassing, this mixture was poured onto an SU-8 patterned silicon wafer and cured for 4 h at 60 °C. Then the PDMS with channel imprints was cut into top and bottom parts, and four inlets were punched (by 0.5 mm diameter biopsy puncher) into the top parts containing the top channel. Dust was removed using Scotch tape (3M). To assemble the two PDMS parts with an embedded membrane, a PDMS/toluene glue at 5:3 wt ratio (1,2) was prepared which ensures complete leakage-free bonding of the two PDMS parts. At first, the glue was spin-coated onto a glass coverslip (1500 rpm for 60 s, ramped at 1000 rpm s -1 ; Spin 150, Polos, the Netherlands) and a thin layer was transferred to both PDMS halves using an ink roller. Then the doubled-porous membrane of 2 mm × 2 mm was cut and placed in the centre of the bottom half. The top part was aligned and gently placed on the bottom part without applying pressure, after which the chips were baked for 3 days at 30 °C (baking was done in low temperature to avoid unwanted damage of the membrane microstructure by melting).

Polyester membrane (8 µm pore size and 10 µm thickness) (3) was also embedded inside the chip system by using the PDMS/toluene glue (Baked in 60 °C for 4 hours) and the resistance of flow was determined as a control with the same procedure. This polyester membrane is typically used in Transwell insert systems for cell culture purposes.

SI 4. Membrane characterization by scanning electron microscopy

After preparation of the membrane with double porosity, the repeatability of the membrane surface and bulk morphology was checked by scanning electron microscopy (SEM; JEOL 5600LV) (Fig. 2). Diameter of the big pores and small porous size distribution were evaluated by ImageJ software on the SEM images.

SI 5. Verification of the absence of leakage with fluorescent dye

After fabrication and before the flux experiments, the chips were verified in order to check good sealing between the membranes and the PDMS or for any unwanted leakages. To do that, a fluorescent dye (4 kDa dextran Alexa 488, 1 mg/ml in PBS) was pipetted inside the channel and the chips were inspected by performing fluorescence microscopy (EVOS FL, Life Technologies) for one day (Supp. Fig. 2).

For this experiment, chips with different batch of membrane were taken and they were imaged in EVOS, where no leakage out of the channel was found just after pipetting the dye. Then the chip were kept in room temperature for one day and imaged again. This time the fluorescent dye was slightly dispersed out of channel from inside the membrane 35 area. The normal hand pressure of pipetting inside the channel must be much higher than the pressure of interests (TEP at 2-16 kPa), suggesting at real time Transendothelial pressure (TEP) the probability of leakage inside the channel is almost none.

SI 6. Permeate flux of the membrane in Amicon cell

To assess the membrane hydraulic performance in a high-pressure filtration system, flux of all the membranes mentioned above were also measured in an Amicon cell at room temperature with deionized water at pressure 20 kPa to 60 kPa (SI Fig. 3). The hydraulic flux and the resistance for water was determined by the above equation where the dynamic viscosity of water at 25 °C was considered as 8.9 × 10 -4 Pa•s. The flux was found much lower than the OOC system at pressure 20 kPa.

The overall flux of all the membranes follow the similar trend. By increasing CHT wt% flux is increasing due to higher hydrophilicity and improved wetting. Unlike the OOC system, flux was not increasing linearly with pressure due to membrane compaction as big pores are present inside the bulk resulting collapsing of the pores through the thickness. By increasing CHT wt% membrane compaction is increasing as CHT in wet condition is mechanically less stable than PCL, wreaking the PCL structural integrity. PCL/CHT 90/10 has shown the highest compressibility due to smallest small pores. Flux of a Polyester membrane (8 µm pore size and 10 µm thickness) was also measured as a control which is around 30 -100 times higher than the blend membrane.

Due to the non-linearity of the flux vs pressure graph in the blend membrane, the hydraulic resistance was determined in each corresponding pressure (and water viscosity at 25 O C) by Eq. 3 and listed in SI. Table . 1. For the polyester membrane (as control), the resistance that obtained from OOC system was quite similar that has been reported in the manufacturer protocol whereas the value in the amicon cell was 3 times higher from the manufacturer protocol. Whereas, in case of PCL-CHT membrane, the resistance in the amicon cell was 100 times higher than that have been obtained in the fluigent system.

Comparison of hydraulic resistance obtained in Fluigent and Amicon cell

Unlike the PCL-CHT membrane, hydraulic flux in the polyester membrane was simpler due to lower hydraulic length (membrane thickness 10 µm) and smooth surface with pores all the way through the membrane thickness. In spite of that, the resistance value in the polyester membrane obtained in the Amicon cell was increasing could be due to added internal instrumental resistance. Whereas PCL-CHT membranes has more complicated surface and bulk morphology than the polyester membrane. Moreover, the blend membrane surface was heterogeneous with combination of pores and pillars where the pillar was decreasing by increasing CHT wt% (the exact proportion of the pores and pillars was unknown). Hence, the resistance that obtained from the Amicon cell is combination of pores and with some pillars. Although similar thing should happen in the OOC system but here the flux experiment was started from 10 times lower pressure and 1256 times lower membrane surface area where it possible that membrane surface majorly faced Pillar (which was giving zero reading and was discarded) and pores (which was giving reading and considered). And as the reading in the OOC was considered majorly from the pores, the resistance should be much lower than the reading of Amicon cell.

6

 6

Fig. 1 .

 1 (A). and the funding organization Erasmus mundus doctorate in membrane engineering-5 th Ed. under European commission (EACEA).

Fig. 1 :

 1 Fig. 1: (A) Organ-on-chip device (inset Picture: BIOS/Lab on a Chip, UTwente); exploded view of the PDMS chip with the top channel (TC), membrane (M), bottom channel (BC) and (PT) blocked pipette tip. The working area of the membrane suspended between two channels is 500 µm × 500 µm. (B) Fluidic set up; determination of volumetric flow rate by perfusion of the cell culture medium (EGM-2) through one top channel and collection of the permeate through the bottom channel to the sensor.

Fig. 2 .

 2 Fig. 2. SEM images of membranes on the Top surface and cross section with double porosity; (A, B) PCL/CHT 100/0, (C, D) PCL/CHT 90/10, (E, F) PCL/CHT 80/20, (G, H) PCL/CHT 70/30 respectively using track-etched membrane (10 µm). All images were taken after removing the track-etched membrane.

Fig. 3 .

 3 Fig. 3. Flux vs Pressure of four different membrane compare to Polyester membrane (control) measured (n=3) in Organ-on-chip device by fluidic set up.

Fig. 4 .

 4 Fig. 4. Flux of four different membranes before cell seeding and after one day and three days of cell seeding. All values were measured at physiological transendothelial pressure (2 kPa). The percentages denote flux reduction compared to the individual blank membrane. (**** p<0.0001, *** 0.0001<p<0.001, ** 0.001<p<0.01 and * 0.01<p<0.05 in one-way ANOVA analysis with Bonferroni correction).

Fig. 5 .

 5 Fig. 5. LCSM representative images of HUVECs on different blends of membranes after staining three different parameters (nucleus by DAPI in blue, cell cytoskeleton by phalloidin 488 in green and cell junctions by VE-Cadherin in red) taken after permeability experiments inside OOC systems after 1 day and 3 day of cell seeding. Bars; 200 µm in inset view (15× mag.) and 20 µm in enlarged view (85× mag.).

  Below is the Resistance (m -1 ) of the Polyester membrane by different report -Manufacturer protocol (3) 4.6 × 10 8 m -1 (at press 70 kPa or 10 psi and Vis. Water in 25 °C) Amicon cell 15.2 × 10 8 m -1 (from slop in Supp. Fig. 3, using Vis. Water in 25 °C) 36 In OOC 5.5 × 10 8 m -1 (from slop in Fig. 3, using Vis. Culture medium in 37 °C)

38 SI. Table. 1 .

 1 Hydraulic resistance (R Mem ) of the blend membranes with distilled water at 25 °C by the Eq. 3.

Table 2 .

 2 Table. 1. Relationship between polymer blend ratio and membrane pore size. Resistance determined by the slop of the Flux vs Pressure graph.

		P C L / C H T	PCL/CHT	PCL/CHT	PCL/CHT
		100/0	90/10	80/20	70/30
	Surface Big	21.4±5.2 14.3±4.1 15.3±3.4 15.8±5.7
	Pores Dia.
	(µm)
	Bulk Big	30.6±4.9 54.4±7.6 51.9±9.2 56.1±5.5
	Pores Dia.
	(µm)
	*Surface	2.3±1.5	0.8±0.3	1.4±0.5	1.7±0.9
	Small Pores
	Dia. (µm)
	Thickness	55	85	65	80
	(µm)

Table 3 .

 3 Resistance of HUVECs cellular barrier determined by eq. 4 and computed in terms of hydraulic conductivity L p after 1 and 3 day.

		Resistance of	Hydraulic conductivity of
		HUVECs (R Cell ) m -1	HUVECs (Lp) cm.s -1 .cmH 2 O -1 ×
		x10 8	10 -2
		1 Day	3 Day	1 Day	3 Day
	PCL/CHT 100/0	4.8±2.1	5.8±3.2	2.96±1.99	2.44±2.22
	PCL/CHT 90/10	7.7±2.5	11.8±0.9	1.84±0.70	1.21±0.09
	PCL/CHT 80/20	4.5±1.8	10.5±2.4	3.14±1.17	1.34±0.33
	PCL/CHT 70/30	2.5±0.5	7.4±2.0	5.72±1.33	1.91±0.5

SI 2. Double porous membrane scaffolds fabrication

The membranes were prepared using a modified liquid induced phase inversion technique using formic acid/acetic acid (FA/AA) as solvent and NaOH as non-solvent (Supp.