# A Robust Seemingly Unrelated Regressions For Row-Wise And Cell-Wise
Contamination

Giovanni Saraceno  Department of Mathematics, University of Trento, Trento,
Italy  Fatemah Alqallaf  Department of Statistics and Operations Research,
Kuwait University, Kuwait  Claudio Agostinelli  Department of Mathematics,
University of Trento, Trento, Italy

###### Abstract

The Seemingly Unrelated Regressions (SUR) model is a wide used estimation
procedure in econometrics, insurance and finance, where very often, the
regression model contains more than one equation. Unknown parameters,
regression coefficients and covariances among the errors terms, are estimated
using algorithms based on Generalized Least Squares or Maximum Likelihood, and
the method, as a whole, is very sensitive to outliers. To overcome this
problem M-estimators and S-estimators are proposed in the literature together
with fast algorithms. However, these procedures are only able to cope with
row-wise outliers in the error terms, while their performance becomes very
poor in the presence of cell-wise outliers and as the number of equations
increases. A new robust approach is proposed which is able to perform well
under both contamination types as well as it is fast to compute. Illustrations
based on Monte Carlo simulations and a real data example are provided.

Keywords: Feasible Generalized Least Squares, Outliers, Robust Statistics,
Two-Step Generalized S-estimator.

##  1 Introduction

The Seemingly Unrelated Regression (SUR) model or Seemingly Unrelated
Regression Equations (SURE), proposed by Zellner [1962], is a generalization
of a linear regression model that consists of several regression equations,
each having its own dependent variable and potentially different set of
exogenous explanatory variables. Each equation is a valid linear regression on
its own and can be estimated separately, however the errors are assumed to be
correlated across the equations.

The model can be estimated equation-by-equation using standard ordinary least
squares (OLS). Such estimates are consistent, however generally not as
efficient as the SUR method, which amounts to feasible generalized least
squares (FGLS) [Zellner, 1962] with a specific form of the variance-covariance
matrix. The SUR model is equivalent to OLS method in two particular cases:
when the errors are uncorrelated between the equations, so that they are truly
unrelated, and when each equation contains exactly the same set of regressors.
Furthermore, it can be viewed as either the simplification of the general
linear model where certain coefficients are restricted to be equal to zero, or
as the generalization of the general linear model where the regressors on the
right hand side are allowed to be different in each equation. Finally, the SUR
model can be further generalized into the simultaneous equations model, where
the regressors are allowed to be the endogenous variables as well. Other
estimation procedures besides FGLS have been proposed for SUR models, see
Kmenta and Gilbert [1968].

However, all these proposed methods are not robust, since they are based on
least squares and outliers, defined as observations separated from the bulk of
data, can highly affect the estimation procedure. Hence, robust alternatives
have been proposed. Koenker and Portnoy [1990] introduced a robust procedure
for the SUR model based on M-estimation, but it is not affine equivariant.
Then, Rousseeuw and Yohai [1984] and Lopuhaä [1989] studied S-estimators for
multivariate location and scatter, whereas Bilodeau and Duchesne [2000]
firstly introduced S-estimation in regression problems. The robust SUR
estimator proposed by Bilodeau and Duchesne [2000] results to be regression
and affine invariant, and has nice robust properties, such as high breakdown
point, but it is computationally expensive. To overcome the computational
issue, Hubert et al. [2017] introduced a fast and robust SUR method based on
Fast-S algorithm [Salibian-Barrera and Yohai, 2006, Salibian-Barrera et al.,
2006].

The robust estimators of SUR model seen so far deal with row-wise
contamination, Tukey-Huber contamination model (THCM), which assumes that a
proportion $\epsilon$ of observations can be contaminated and these
independent events are the units considered as outliers. Alqallaf et al.
[2009] consider a different contamination model for multivariate data: the
independent contamination model (ICM) or cell-wise contamination, where the
entries of an observation (or cells) can be independently contaminated.
According to this paradigm, given a fraction $\epsilon$ of contaminated cells,
the probability that at least one component of an observation is an outlier is
$1-(1-\epsilon)^{p}$, where $p$ is the dimension of observations. This number
is close to one when $p$ is large even if $\epsilon$ is small. For this reason
estimators that have breakdown point 0.5 under the THCM, may have breakdown
tending to zero under the ICM. Alqallaf et al. [2009] show that this happens
with the most popular high breakdown point equivariant estimators of
multivariate location, e.g., S-estimators[Davies, 1987], Minimum Volume
Ellipsoid [Rousseeuw, 1985], Minimum Covariance Determinant [Rousseeuw, 1985]
or the Stahel-Donoho estimators [Donoho, 1982, Stahel, 1981]. In order to cope
with both row-wise and cell-wise contamination, Agostinelli et al. [2015]
developed the 2SGS-estimators for multivariate location and scatter while
Leung et al. [2017] further extended the approach to linear models.

Here, we propose a robust estimator for the SUR model, which we will refer to
as surerob, able to deal with both row-wise contamination and cell-wise
contamination by robustfying the FGLS approach of Zellner [1962]. To the best
of our knowledge, there are no competiting estimators in literature which are
able to deal with cell-wise outliers for SUR models.

The remainder of the paper is organized as follows. The SUR model is presented
in Section 2, while Section 3 describes the proposed robust estimation method.
The performance of the estimators, in case of row-wise and cell-wise
contamination, is investigated through a simulation study in Section 4 and
with a real data example in Section 5. Concluding remarks end the paper in
Section 6.

##  2 The SUR model

Consider $m$ regression equations

$y_{ik}=\mathbf{x}_{ik}^{\top}\boldsymbol{\beta}_{i}+\varepsilon_{ik}\ ,\qquad
i=1,\ldots,m\ .$

Here, $i$ represents the equation index, while $k=1,\ldots,n$ is the
observation index and $n$ is the total number of observations. Each equation
$i$ has a single response variable
$\mathbf{y}_{i}=(y_{i1},\ldots,y_{in})^{\top}$, errors vector
$\boldsymbol{\varepsilon}_{i}=(\varepsilon_{i1},\ldots,\varepsilon_{in})^{\top}$,
an explanatory $n\times p_{i}$ matrix $X_{i}$ and $p_{i}$-vector
$\boldsymbol{\beta}_{i}$ of coefficients which lead to the form

$\mathbf{y}_{i}=X_{i}\boldsymbol{\beta}_{i}+\boldsymbol{\varepsilon}_{i}\
,\qquad i=1,\ldots,m\ .$

Finally, if we stack these $m$ vector equations on top of each other, the
system will take the form [Zellner, 1962, eq. (2.2)]

$\displaystyle\begin{pmatrix}\mathbf{y}_{1}\\\ \mathbf{y}_{2}\\\ \vdots\\\
\mathbf{y}_{m}\end{pmatrix}$ $\displaystyle=\begin{pmatrix}X_{1}&0&\ldots&0\\\
0&X_{2}&\ldots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\ldots&X_{m}\end{pmatrix}\begin{pmatrix}\boldsymbol{\beta}_{1}\\\
\boldsymbol{\beta}_{2}\\\ \vdots\\\
\boldsymbol{\beta}_{m}\end{pmatrix}+\begin{pmatrix}\boldsymbol{\varepsilon}_{1}\\\
\boldsymbol{\varepsilon}_{2}\\\ \vdots\\\
\boldsymbol{\varepsilon}_{m}\end{pmatrix}$ (1) $\displaystyle\mathbf{y}$
$\displaystyle=X\boldsymbol{\beta}+\boldsymbol{\varepsilon}\ .$

where $\boldsymbol{\beta}$ is a vector of dimension $P=\sum_{i=1}^{m}p_{i}$.
The assumption of the model is that error terms $\varepsilon_{ik}$ are
independent across observations, but may have cross-equation contemporaneous
correlations, that is, $\mathbb{E}(\varepsilon_{ir}\varepsilon_{is}|X)=0$
whenever $r\neq s$, whereas
$\mathbb{E}(\varepsilon_{ik}\varepsilon_{jk}|X)=\sigma_{ij}$. Note that
$\sigma_{ii}$ is the variance of the error term in the $i$th equation, whereas
$\sigma_{ij}$ is the covariance between the errors in equations $i$ and $j$.
Let $\Sigma$ denote the $m\times m$ covariance matrix of each observation with
entries $\sigma_{ij}$, then the covariance matrix of the stacked error terms
$\boldsymbol{\varepsilon}$ will be equal to

$\Omega=\mathbb{E}(\boldsymbol{\varepsilon}\boldsymbol{\varepsilon}^{\top}|X)=\Sigma\otimes
I_{n},$

where $I_{n}$ is the $n\times n$ identity matrix and $\otimes$ denotes the
matrix Kronecker product.

An alternative formulation of the multivariate SUR model is given by

$\tilde{Y}=\tilde{X}B+E$ (2)

where $\tilde{Y}=(\mathbf{y}_{1},\ldots,\mathbf{y}_{m})$ is the $n\times m$
response matrix, $\tilde{X}=(X_{1},\ldots,X_{m})$ is the $n\times P$ (recall
that $P=\sum_{i=1}^{m}p_{i}$) design matrix,
$B=\operatorname{diag}(\boldsymbol{\beta}_{1},\ldots,\boldsymbol{\beta}_{m})$
is a $P\times m$ block diagonal matrix, and
$E=(\boldsymbol{\varepsilon}_{1},\ldots,\boldsymbol{\varepsilon}_{m})$ with
$\operatorname{\mathbb{C}ov}(E)=\Omega=\Sigma\otimes I_{n}$.

Let $\mathbf{T}=(\mathbf{t}_{1},\ldots,\mathbf{t}_{n})^{\top}$ be a data set
of size $n$ corresponding to model (2), where
$\mathbf{t}_{k}=(\mathbf{t}_{1k},\ldots,\mathbf{t}_{mk})$ and
$\mathbf{t}_{ik}=(y_{ik},x_{ik1},\ldots,x_{ikp_{i}})$, $1\leq i\leq m$, $1\leq
k\leq n$, and $x_{ikp}$ is the value in the $k$-th row and $p$-th column of
the matrix $X_{i}$. According to THCM a proportion $\epsilon$ of observations
$\mathbf{t}_{k}$, $1\leq k\leq n$, can be contaminated, while following ICM
the entries $\mathbf{t}_{ik}$, $1\leq k\leq n$, $1\leq i\leq m$ can be
independently replaced by an outlier.

##  3 Robust estimation

We discuss the two main estimation methods used for the SUR model: the
standard FGLS method [Zellner, 1962] and the robust SUR method base on Fast-S
algorithm [Hubert et al., 2017], which we will refer to as fastSUR. Finally,
we introduce the proposed new estimator for the SUR model robust against row-
wise and cell-wise outliers.

The SUR model is usually estimated using FGLS method [Zellner, 1962]. This is
a two-step method where in the first step an ordinary least squares regression
is performed for each model equation separately. The residuals
$\hat{\boldsymbol{\varepsilon}_{i}}$ from these regressions are used to
estimate the elements of the matrix $\Sigma$ by computing the sample
covariance matrix $\hat{\Sigma}_{1}$ with components

$\hat{\sigma}_{ij}=\frac{1}{n}\hat{\boldsymbol{\varepsilon}}_{i}^{\top}\hat{\boldsymbol{\varepsilon}}_{j}\qquad
i,j=1,\ldots,m\ .$ (3)

In the second step a generalized least squares regression is performed using
the variance matrix $\hat{\Omega}=\hat{\Sigma}_{1}\otimes I_{n}$ to obtain
final estimates of the regression coefficients as

$\hat{\boldsymbol{\beta}}=\left(X^{\top}\hat{\Omega}X\right)^{-1}X^{\top}\hat{\Omega}^{-1}\mathbf{y}\
.$ (4)

Residuals can then be recomputed and used to obtain a final estimate
$\hat{\Sigma}_{2}$ of the covariance matrix $\Sigma$. The FGLS algorithm is
available in R [R Core Team, 2021] in package systemfit, see Henningsen and
Hamann [2007].

The fastSUR, which is the computationally efficient version of the robust SUR
method of Bilodeau and Duchesne [2000], is the couple $(\hat{B},\hat{\Sigma})$
so that

$(\hat{B},\hat{\Sigma})=\arg\min_{(B,S)}|S|$

under the condition

$\frac{1}{n}\sum_{k=1}^{n}\rho\left(\mathbf{e}_{i}(B)^{\top}S^{-1}\mathbf{e}_{i}(B)\right)=b$

where $|\cdot|$ is the determinant,
$B=\operatorname{diag}(\mathbf{\beta}_{1},\ldots,\mathbf{\beta}_{m})$,
$\mathbf{\beta}_{i}\in\mathbb{R}^{p_{i}}$ ($i=1,\ldots,m$),
$\mathbf{e}_{i}(B)=\tilde{Y}-\tilde{X}B$, and $S$ is an $m\times m$ symmetric
positive definite matrix. The function $\rho$ [see, e.g., Maronna et al.,
2018] is chosen so that

  * (C1)

$\rho$ is symmetric around zero and twice continuously differentiable;

  * (C2)

$\rho(0)=0$ and $\rho$ is strictly increasing on $[0,c_{0}]$ and constant on
$[c_{0},\infty]$ for some $c_{0}>0$.

The constant $b$ can be computed as $\mathbb{E}_{F_{0}}(\rho(|\mathbf{e}|))$,
where $\mathbf{e}\sim F_{0}$ and $F_{0}=N_{m}(\mathbf{0},I_{m})$ which ensures
consistency at the model with normal errors. Similar to classic robust
location and scatter estimators, robust SUR and its computationally efficient
version fastSUR are affine equivariant and their breakdown point tends to zero
as the number of equations $m$ gets larger under the Independent Contamination
Model.

We propose a robust estimator for the SUR model under both types of
contamination by robustfying the FGLS approach of Zellner [1962]. In the first
step, residuals $\hat{\boldsymbol{\varepsilon}_{i}}$ are estimated by means of
an MM-estimator of regression, as introduced in Yohai [1987], equation by
equation. Since these are univariate regression models, the procedure achieves
a breakdown of $0.5$ and the estimates of the regression coefficients are
affine equivariant. At this point, the ICM contamination scheme is considered,
therefore we construct the residual data matrix
$\hat{E}=(\hat{\boldsymbol{\varepsilon}}_{1},\ldots,\hat{\boldsymbol{\varepsilon}}_{m})$
where each cell of this matrix could be a cell-wise outlier. Let
$W_{i}=\operatorname{diag}(\mathbf{w}_{i})$ be an $n\times n$ diagonal matrix
where $\mathbf{w}_{i}$ is the vector of robust weights associated to each
observation in the $i$th equation, that is, the $k$th element $w_{ik}$ of the
vector $\mathbf{w}_{i}$ is given by

$w_{ik}=\left\\{\begin{array}[]{ll}\psi(\hat{\varepsilon}_{ik}/s_{i})/(\hat{\varepsilon}_{ik}/s_{i})&\text{if}\
\hat{\varepsilon}_{ik}\neq 0\\\ 1&\text{otherwise}\end{array}\right.$ (5)

where $\psi=\rho^{\prime}$ is the first derivative of the function $\rho$ and
$s_{i}$ is the estimated standard deviation of the errors for the $i$th
equation. A robust estimate $\hat{\Sigma}_{1}$ of the covariance matrix
$\Sigma$ is obtained using the 2SGS method based on the residual data matrix
$\hat{E}$. The final estimate of the coefficients is than obtained as

$\hat{\boldsymbol{\beta}}=\left(X^{\top}W(\hat{\Sigma}_{1}^{-1}\otimes
I_{n})WX\right)^{-1}X^{\top}W(\hat{\Sigma}_{1}^{-1}\otimes I_{n})W\mathbf{y}\
,$

where $W$ is a block diagonal matrix with $(W_{1},\ldots,W_{m})$ in the main
diagonal. New residuals can then be obtained and, applying 2SGS to them, we
get the final estimate $\hat{\Sigma}_{2}$ of the covariance matrix.

##  4 Simulation Study

The performance of the introduced method, which we will refer to as surerob,
is compared with classical FGLS algorithm, as implemented in the R package
systemfit [Henningsen and Hamann, 2007], indicated as sure, and the fastSUR
algorithm [Hubert et al., 2017]. The fastSUR algorithm [Hubert et al., 2017]
uses the Tukey’s bisquare function with constant $b$ such that the breakdown
point is $0.5$ [Rousseeuw and Yohai, 1984] as it is implemented in an R code
kindly made available by Prof. M. Hubert.

The R implementation of the surerob procedure uses the function lmrob in the R
package robustbase [Mächler et al., 2021] and the function TSGS in the R
package GSE [Leung et al., 2019] using default values for both. The code is
available in the R package robustsur provided as supplementary material. An
allustration of how to use the functions in the R package robustsur is
reported in Section SM–1 of the Supplemental Material.

The simulation has the following setting: sample size $n=100$; $p=p_{i}=5$,
$10$, equals for each equation and $m=5$, $10$, $20$. For each combination of
these factors we run $N=1000$ Monte Carlo replications. The regression
coefficients $\boldsymbol{\beta}_{i}$ are sampled from a standard Cauchy
random variable for each sample while the variances in
$\operatorname{diag}(\Sigma)$ are all equal to $1$, that is, $\Sigma$ is a
correlation matrix. To account for the lack of affine equivariance of the
proposed estimator, we consider different correlation structures. In
particular, for each sample in our simulation we create a different random
correlation matrix with condition number fixed at $\text{CN}=100$. Correlation
matrices with high condition number are less favorable for our proposed
estimator. For the details about the procedure used to obtain such random
correlation matrices, see Agostinelli et al. [2015].

Two types of outliers are considered: (i) generated by THCM and (ii) generated
by ICM. When the outliers are generated using THCM, we randomly replace $5\%$,
$10\%$, $20\%$ and $30\%$ of the cases in the errors data matrix by
$k\mathbf{v}$, where $k=0,5,10,\ldots,95,100$ and $\mathbf{v}$ is the
eigenvector corresponding to the smallest eigenvalue of $\Sigma$ with length
such that $\mathbf{v}^{\top}\Sigma^{-1}\mathbf{v}=1$. Monte Carlo experiments
in Agostinelli et al. [2015] show that the placement of outliers in this
direction, $\mathbf{v}$, is the least favorable for the 2SGS estimator. When
the outliers are generated using ICM, we randomly replace $5\%$, $10\%$,
$20\%$ or $30\%$ of the cells in the error data matrix by the value $k$ where
$k=1,5,\dots,95,100$.

The performance of a given regression coefficients estimator
$\hat{\boldsymbol{\beta}}$ is measures by Mean Square Error

$\operatorname{MSE}(\hat{\boldsymbol{\beta}},\boldsymbol{\beta})=\frac{1}{N}\sum_{r=1}^{N}(\hat{\boldsymbol{\beta}}_{r}-\boldsymbol{\beta})^{\top}(\hat{\boldsymbol{\beta}}_{r}-\boldsymbol{\beta})$

where $\hat{\boldsymbol{\beta}}_{r}$ is the estimate at the $r$-th
replication.

The performance of a given scatter estimator $\hat{\Sigma}$ is measured by the
Kullback-Leibler divergence between two Gaussian distribution with the same
mean and covariances $S$ and $\Sigma$:

$\delta(S,\Sigma)=\operatorname{trace}(S\Sigma^{-1})-\log(|S\Sigma^{-1}|)-p\
.$

This divergence also appears in the likelihood ratio test statistics for
testing the null hypothesis that a multivariate normal distribution has
covariance matrix $\Sigma$. Then, the performance of an estimator
$\hat{\Sigma}$ is summarized by

$\Delta(\hat{\Sigma},\Sigma)=\frac{1}{N}\sum_{r=1}^{N}\delta(\hat{\Sigma}_{r},\Sigma)$

where $\hat{\Sigma}_{r}$ is the estimate at the $r$-th replication.

Figures 1-2 report the results for the case $p=5$, $m=10$ for THCM, while
Figures 3-4 report the results for ICM. Results are similar for all the other
cases and they are not reported. Complete results are available in Section
SM–2 of the Supplemental Material. As expected, the sure method is sensitive
to the presence of contamination. The fastSUR has a better performance than
surerob only for the THCM with low level of contamination, says under $10\%$,
while in all other cases the procedure breaks down and often performs slightly
worse than the classical sure procedure. When ICM is considered surerob
outperforms the other estimators.

Figure 1: Mean Square Error
$\operatorname{MSE}=\operatorname{MSE}(\hat{\boldsymbol{\beta}},\boldsymbol{\beta})$
for surerob (solid line), sure (dashed line) and fastSUR (dotted line) for
different levels of contamination $\epsilon=5\%$, $10\%$, (top) $20\%$, $30\%$
(bottom) under the THCM considering dimension $p=5$ and number of equations
$m=10$.

Figure 2: Average Kullback-Leibler
$\Delta_{2}=\Delta(\hat{\Sigma}_{2},\Sigma)$ for surerob (solid line), sure
(dashed line) and fastSUR (dotted line) for different levels of contamination
$\epsilon=5\%$, $10\%$, (top) $20\%$, $30\%$ (bottom) under the THCM
considering dimension $p=5$ and number of equations $m=10$.

Figure 3: Mean Square Error
$\operatorname{MSE}=\operatorname{MSE}(\hat{\boldsymbol{\beta}},\boldsymbol{\beta})$
for surerob (solid line), sure (dashed line) and fastSUR (dotted line) for
different levels of contamination $\epsilon=5\%$, $10\%$, (top) $20\%$, $30\%$
(bottom) under the ICM considering dimension $p=5$ and number of equations
$m=10$.

Figure 4: Average Kullback-Leibler $\Delta_{2}=\Delta(\hat{\Sigma}_{2},\Sigma)$ for surerob (solid line), sure (dashed line) and fastSUR (dotted line) for different levels of contamination $\epsilon=5\%$, $10\%$, (top) $20\%$, $30\%$ (bottom) under the ICM considering $p=5$ and number of equations $m=10$. Table 1: Mean execution time (seconds) for the three methods under different contamination type and levels, for $p=5$ and $m=10$. |  | Contamination level
---|---|---
Method | Type | $0\%$ | $5\%$ | $10\%$ | $20\%$ | $30\%$
surerob | THCM | 5.89 | 5.83 | 5.74 | 5.62 | 5.82
| ICM | – | 6.82 | 7.68 | 9.62 | 10.87
sure | THCM | 0.36 | 0.36 | 0.35 | 0.36 | 0.35
| ICM | – | 0.36 | 0.36 | 0.36 | 0.35
fastSUR | THCM | 217.61 | 218.84 | 220.04 | 223.99 | 228.82
| ICM | – | 222.58 | 218.75 | 217.29 | 212.98

Table 1 reports the mean execution time, in seconds, for the three procedures
for $p=5$ and $m=10$. While the classic sure method is the fastest in all the
cases, our procedure is very competitive, instead the fastSUR is by far the
slowest procedure. We tried to extend the simulation study to the case
$n=500$, but the computational time needed by the fastSUR method was
prohibitive.

##  5 Real data example

In this section the introduced robust method is performed on data about
tourism and compared with the results obtained with the sure method, from the
R package systemfit, and the fastSUR algorithm.

Disegna and Osti [2016] collected data to study the relationship between
satisfaction and tourism expenditure, as well as the dependence among
different tourism expenditure categories. In particular, the aim was twofold:
to investigate the influence on tourism expenditure of tourists’ satisfaction
with the destination, considering some socio-demographic and trip-related
variables; second, to study the dependence among tourist expenditure on
different categories. Data were collected through a survey conducted by the
Bank of Italy (Banca d’Italia) which includes socio-demografic
characteristics, information on the trip, information on the expenditure,
level of satisfaction with different aspects of the trip and overall
satisfaction with the destination. In their study Disegna and Osti [2016]
focused on 1030 foreign visitors who visited the provinces of Bolzano, Trento
and Belluno in 2011 with the main purpose of trip being “tourism”, “holiday”
or “leisure”. For each observation, we have detailed information on the money
spent with respect to five categories: “Accommodation”, “Food and Beverages”,
“International transportation”, “Shopping” and “Other services” (such as
museum, excursions, shows and so on). For a complete description of the survey
and the descriptive analysis of the data set, see Disegna and Osti [2016].

We considered four regression equations, one for each expenditure category
excluding the “Other services” category. The explanatory variables considered
include the satisfaction with respect to 10 classes, destination, number of
nights, age and if the person is traveling alone, for a total of 20
explanatory variables. Each equation shares the same covariates except for one
which express the total tourism expenditure in the remaining categories (the
“Other services” category is included in these summations).

Figure 5: On the left, differences between the estimates obtained by surerob
and the estimates obtained by sure. On the right, differences between
estimates given by fastSUR and sure methods. Figure 6: Robust weights assigned
to observations for each equation by the surerob method.

Figure 5 shows the differences between the estimates obtained by the surerob
and sure, on the left, and the differences between estimates given by fastSUR
and sure, on the right. We added dashed red lines highlighting the interval
$(-0.1,0.1)$ on the $y$-axis. Parameters estimated by surerob are quite
similar to those given by the traditional sure, with some exceptions, while
the differences computed with respect to the fastSUR method show more
variability.

The surerob estimation procedure is constructed assigning cell-wise weights,
defined in equation 5, associated to observations for each equation, which are
displayed in Figure 6. Weights have been sorted to make the amount of
downweighted observations more clear. Assume that we flag a cell as outliers
if the corresponding weight is less than 0.5. In this case, surerob identifies
$2.3\%$ of cell-wise contamination that propagates to the $8.8\%$ of rows.

Finally, we computed $p$-values and $R^{2}$ for the estimates obtained by sure
and surerob, with respect to the whole system and single regression equations.
The fastSUR method has not been considered since this part has not been
implemented.

Table 2: P-values computed using the estimates given by sure and surerob
methods for the single regression equations.

| Eq.1 (Accommodation) | Eq.2 (Transportation) | Eq.3 (Food) | Eq.4 (Shopping)
---|---|---|---|---
| sure | surerob | sure | surerob | sure | surerob | sure | surerob
(Intercept) | 0.265 | 0.323 | 0.309 | 0.168 | 0.003 | 0.003 | 0.110 | 0.264
luoghi_visitati1 | 0.870 | 0.664 | 0.000 | 0.000 | 0.109 | 0.005 | 0.073 | 0.000
nr_notti | 0.007 | 0.000 | 0.001 | 0.000 | 0.062 | 0.000 | 0.171 | 0.000
viaggia_solo_971 | 0.026 | 0.024 | 0.000 | 0.000 | 0.084 | 0.458 | 0.654 | 0.103
voto_cortesia | 0.674 | 0.334 | 0.141 | 0.058 | 0.057 | 0.324 | 0.702 | 0.433
voto_arte | 0.769 | 0.807 | 0.948 | 0.834 | 0.731 | 0.910 | 0.604 | 0.814
voto_ambiente | 0.001 | 0.000 | 0.002 | 0.000 | 0.872 | 0.164 | 0.013 | 0.225
voto_alberghi | 0.760 | 0.738 | 0.027 | 0.023 | 0.150 | 0.144 | 0.556 | 0.984
voto_pasti | 0.398 | 0.098 | 0.014 | 0.008 | 0.214 | 0.699 | 0.559 | 0.582
voto_prezzi | 0.046 | 0.001 | 0.07640 | 0.019 | 0.116 | 0.699 | 0.660 | 0.246
voto_acquisti | 0.205 | 0.023 | 0.937 | 0.607 | 0.257 | 0.308 | 0.004 | 0.000
voto_informazioni | 0.419 | 0.251 | 0.375 | 0.112 | 0.754 | 0.627 | 0.551 | 0.201
voto_sicurezza | 0.292 | 0.250 | 0.051 | 0.091 | 0.714 | 0.991 | 0.089 | 0.007
voto_complessivo | 0.480 | 0.546 | 0.229 | 0.156 | 0.467 | 0.354 | 0.689 | 0.963
eta_cod11 | 0.140 | 0.461 | 0.115 | 0.034 | 0.427 | 0.881 | 0.046 | 0.026
eta_cod21 | 0.167 | 0.604 | 0.736 | 0.634 | 0.492 | 0.638 | 0.352 | 0.279
eta_cod31 | 0.028 | 0.049 | 0.742 | 0.687 | 0.683 | 0.906 | 0.355 | 0.296
stato_21 | 0.000 | 0.000 | 0.000 | 0.000 | 0.079 | 0.000 | 0.000 | 0.714
stato_31 | 0.420 | 0.289 | 0.236 | 0.783 | 0.315 | 0.782 | 0.063 | 0.003
stato_41 | 0.013 | 0.020 | 0.821 | 0.362 | 0.096 | 0.978 | 0.121 | 0.440
total_sum | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000

Table 3: R-squared computed using the estimates given by sure and surerob methods for the single regression equations and for the system. | sure | surerob
---|---|---
| $R^{2}$ | adj $R^{2}$ | $R^{2}$ | adj $R^{2}$
Eq.1 (Accommodation) | 0.465 | 0.454 | 0.527 | 0.518
Eq.2 (Transportation) | 0.340 | 0.331 | 0.372 | 0.359
Eq.3 (Food) | 0.344 | 0.331 | 0.378 | 0.366
Eq.4 (Shopping) | 0.154 | 0.137 | 0.188 | 0.172
System | 0.325 | - | 0.457 | -

Table 2 reports the $p$-values for each covariate considering the equations
separately. Remember that the last variables is the only one that differs
among equations. Considering a $95\%$ confidence level, significant $p$-values
are written in bold font. Table 3 reports the $R^{2}$ values computed for sure
and surerob with respect to the whole system and single equations. Robust
estimates lead to a set of significant variables slightly different from that
identified by standard sure, indeed some variables with a large $p$-value for
sure become significant for surerob and vice versa.

##  6 Conclusions

We proposed a new robust estimation method for the SUR model considering both
row-wise and cell-wise outliers. Under the THCM, our estimator outperforms the
robust competitor for high contamination levels and it remains competitive for
low levels of contamination. It is worth remarking that, even if it is slower
than the classic sure method, the surerob method is less computational
expensive than the other robust estimators. This is an appealing property when
$n$ increases. When cell-wise contamination is considered, the proposed
estimator shows the best performance.

## References

  * Agostinelli et al. [2015] C. Agostinelli, A. Leung, V.J. Yohai, and R.H. Zamar.  Robust estimation of multivariate location and scatter in the presence of cellwise and casewise contamination.  _TEST_ , 24(3):441–461, 2015.
  * Alqallaf et al. [2009] F. Alqallaf, S. Van Aelst, R.H. Zamar, and V.J. Yohai.  Propagation of outliers in multivariate data.  _The Annals of Statistics_ , 37(1):311–331, 2009\.
  * Bilodeau and Duchesne [2000] M. Bilodeau and P. Duchesne.  Robust estimation of the SUR model.  _Canadian Journal of Statistics_ , 28(2):277–288, 2000.  doi: 10.2307/3315978.
  * Davies [1987] P.L. Davies.  Asymptotic behaviour of S-estimates of multivariate location parameters and dispersion matrices.  _The Annals of Statistics_ , 15(3):1269–1292, 1987.
  * Disegna and Osti [2016] M. Disegna and L. Osti.  Tourists’ expenditure behaviour: The influence of satisfaction and the dependence of spending categories.  _Tourism Economics_ , 22(1):5–30, 2016.  doi: 10.5367/te.2014.0410.
  * Donoho [1982] D.L. Donoho.  Breakdown properties of multivariate location estimators.  _Qualifying paper, Harvard University, Boston_ , 1982.
  * Henningsen and Hamann [2007] A. Henningsen and J.D. Hamann.  systemfit: A package for estimating systems of simultaneous equations in R.  _Journal of Statistical Software_ , 23(4):1–40, 2007.  URL http://www.jstatsoft.org/v23/i04/.
  * Hubert et al. [2017] M. Hubert, T. Verdonck, and Ö. Yorulmaz.  Fast robust SUR with economical and actuarial applications: Fast robust SUR.  _Statistical Analysis and Data Mining: The ASA Data Science Journal_ , 10(2):77–88, 2017.  doi: 10.1002/sam.11313.
  * Kmenta and Gilbert [1968] J. Kmenta and R.F. Gilbert.  Small sample properties of alternative estimators of seemingly unrelated regressions.  _Journal of the American Statistical Association_ , 63(324):1180–1200, 1968.
  * Koenker and Portnoy [1990] R. Koenker and S. Portnoy.  M estimation of multivariate regressions.  _Journal of the American Statistical Association_ , 85(412):1060–1068, 1990.
  * Leung et al. [2017] A. Leung, V.J. Yohai, and R.H. Zamar.  Multivariate location and scatter matrix estimation under cellwise and casewise contamination.  _Computational Statistics & Data Analysis_, 111:59–76, 2017.
  * Leung et al. [2019] A. Leung, M. Danilov, V.J. Yohai, and R.H. Zamar.  _GSE: Robust Estimation in the Presence of Cellwise and Casewise Contamination and Missing Data_ , 2019.  URL https://cran.r-project.org/web/packages/GSE.  R package version 4.2.
  * Lopuhaä [1989] H. P. Lopuhaä.  On the relation between s-estimators and m-estimators of multivariate location and covariance.  _Annals of Statistics_ , 17:1662–1683, 1989.
  * Mächler et al. [2021] M. Mächler, P.J. Rousseeuw, C. Croux, V. Todorov, A. Ruckstuhl, M. Salibian-Barrera, T. Verbeke, M. Koller, E.L.T. Conceicao, and M. Anna di Palma.  _robustbase: Basic Robust Statistics_ , 2021.  URL http://robustbase.r-forge.r-project.org/.  R package version 0.93-8.
  * Maronna et al. [2018] R.A. Maronna, R.D. Martin, and Yohai V.J.  _Robust statistic: theory and methods_.  Wiley, Chichister, 2018.  ISBN 978-1-119-21468-7.
  * R Core Team [2021] R Core Team.  _R: A Language and Environment for Statistical Computing_.  R Foundation for Statistical Computing, Vienna, Austria, 2021.  URL https://www.R-project.org/.
  * Rousseeuw [1985] P.J. Rousseeuw.  Multivariate estimation with high breakdown point.  In W. Grossmann, G. Pflug, I. Vincze, and Wertz. W., editors, _Mathematical Statistics and Applications_ , pages 283–297. Reidel Publishing Company (co-published with Akadémiai e Kiadó, Budapest), Dordrecht, 1985.
  * Rousseeuw and Yohai [1984] P.J. Rousseeuw and V.J. Yohai.  Robust regression by means of s-estimators.  In J. Franke, W. Härdle, and R. D. Martin, editors, _Robust and Nonlinear Time Series Analysis_ , volume 26 of _Lecture Notes in Statistics_ , pages 256–272. Springer-Verlag, New York, 1984.
  * Salibian-Barrera and Yohai [2006] M. Salibian-Barrera and V.J. Yohai.  A fast algorithm for s-regression estimates.  _Journal of Computational and Graphical Statistics_ , 15:414–427, 2006.
  * Salibian-Barrera et al. [2006] M. Salibian-Barrera, S. Van Aelst, and G. Willems.  Pca based on multivariate mm-estimators with fast and robust bootstrap.  _Journal of the American Statistical Association_ , 101:1198–1211, 2006.
  * Stahel [1981] W.A. Stahel.  _Robuste Schätzungen: infinitesimale Optimalität und Schätzungen von Kovarianzmatrizen_.  PhD thesis, ETH Zürich, 1981.
  * Yohai [1987] V.J. Yohai.  High breakdown-point and high efficiency estimates for regression.  _The Annals of Statistics_ , 15:642–665, 1987.
  * Zellner [1962] A. Zellner.  An efficient method of estimating seemingly unrelated regressions and tests for aggregation bias.  _Journal of the American Statistical Association_ , 57(298):348–368, 1962.  doi: 10.2307/2281644.


# Using Glowscript to Teach Numerical Modeling in Undergraduate Biology
Education

Joshua G. Schreibeis, Olivia M. Merideth, and Gavin A. Buxton

Science Department, Robert Morris University, Moon Township, PA 15108, US.

Mathematical and numerical modeling is an increasingly important, yet often
neglected, topic for biology students. We have found Glowscript to facilitate
teaching and introducing computer simulations to students. In particular, the
built-in the graphics and graphing capabilities can provide students with
immediate feedback. Glowscript is a web-based form of visual Python that runs
in a standard browser, and students can easily embed their simulations in
external websites. Here we show various examples of how Glowscript is
implemented in an undergraduate computational biology course. We hope these
examples inspire others to adopt Glowscript in their science classrooms.

## Introduction

Glowscript (or web VPython) is a simple and versatile programming language for
introducing undergraduate biology students to mathematical modeling and
computer simulations. The increasing use of numerical and mathematical
modeling in biological research is not reflected in most undergraduate biology
programs [10]. Introducing computational biology, and computer programming in
general, to undergraduate biology students can help foster computational
thinking, familiarize students with computer research in their fields of
interest, and open up opportunities for students to use computer models in
their future careers [8, 7, 2]. However, teaching computer programming to
undergraduate biology students must stress the relevance of computational
biology through applications that compliment their other studies and
interests. Here we describe the use of Glowscript for introducing
undergraduate biology students, with no prior programming experience, to
computer programming for the first time and to the use of numerical and
mathematical modeling in biological research.

Glowscript is a web-based version of visual Python (web VPython) that includes
a rich selection of visual objects and a graphing capability that can be used
for displaying the results of a numerical simulation while the simulation is
running. Python has been a popular choice for introducing students to computer
programming due to the intuitive syntax [13, 9, 12, 11], the ease of defining
variables, and because it is an open and free software [9, 12]. That said,
Python has been criticized as it overloads operators in a manner that can be
be confusing given dynamic typing, especially when dealing with lists [11],
and is not as efficient or optimized as other more syntax-heavy languages (for
example, the use of structures and memory allocation in C can make it much
more efficient). However, Python is a wonderful computer language for
introducing students to computer programming and it forces the programmer to
implement a clear indented structure that facilitates debugging [13, 9]; this
is true for both the student and the instructor (who more often than not can
be found wandering around the room putting out fires). This makes the language
not only ideal for students to learn, but for teachers to teach. In addition,
the fact that Glowscript offers high-level graphics and graphing capabilities,
is entirely web-based with the students running the code within a standard
browser, and that the programs can be easily exported and embedded in a
website, makes Glowscript an invaluable improvement to the Python language
when it comes to both classroom instruction and student adoption outside of
class.

Here we introduce a couple of different numerical models, implemented using
Glowscript, that can be used in a computational biology classroom; models that
capture simple logistic growth, epidemiology, microtubule dynamic
instabilities and drug diffusion from a polymer nanoparticle. Throughout the
implementation of different models in this course, the students are introduced
to a basic model, which we will describe in detail here, and then they are
required to creatively adapt or improve the model in some unique way.

## Example: Simple Logistic Growth

As a simple example, the logistic growth model is a wonderful introduction to
using a computer program to capture a biological process. The logistic growth
rate allows for the growth of a population that is constrained by resources or
other limitations of its environment. The model describes the growth rate of
many smaller organisms (such as bacteria, various amoeboid organisms, and
diatoms) reasonably well [17]. The population dynamics are described as

$\frac{dN}{dt}=rN\left(\frac{K-N}{K}\right)$

Where $N$ is the number of organisms, $t$ is time, $r$ is the growth rate and
$K$ is the carrying capacity. The growth of the organism is initially
proportional to the number of organisms present, but will be ultimately
limited by the carrying capacity, the number of organisms that the environment
can support. The finite difference approximation of the above equation is
given by

$N^{t+1}=N^{t}+\Delta t\,rN^{t}\left(\frac{K-N^{t}}{K}\right)$

Where $\Delta t$ is the time step in the model and the time is indicated
through the superscript on $N$.

The implementation of this model in Glowscript is shown below.

⬇

GlowScript 3.2 VPython

The above command loads the glowscript libraries. The constants are created
easily in Python, without the need to first declare the variables or their
type.

⬇

K = 1000

r=1

N0 = 1

deltat = 0.1

N = N0

t = 0.

The data can be dynamically plotted, as the simulation is in progress. First,
the graph is created and stored as a variable plot in case we want to change
any of the properties (which we currently don’t). The title, axis, and size of
the graph are prescribed. Next, the function that we are plotting is defined
as a curve that is red with a width of 2 and labeled “Population”; even though
there is only going to be one curve, Glowscript will include a key.

⬇

plot = graph(title='Logistic Growth.', xtitle='time', ytitle='N', fast=False,
width=1000,

height=500)

funct1 = gcurve(color=color.red, width=2, label='Population')

The simulation can now be ran, which consists of iterations within a while
loop. During each time step the population is updated and plotted as a
function of time.

⬇

t=0

while t < 200*deltat:

rate(500)

t = t+deltat

oldN = N

N = oldN + deltat*r*oldN*(K - oldN)/K

funct1.plot( t, N )

Running the simulation, therefore, dynamically updates the graph of population
versus time, and concludes with the following plot.

Figure 1: A screenshot of the Glowscript implementation. The output from the
logistic growth model depicts the graphing capabilities of Glowscript. The
graph is updated dynamically as the simulation progresses allowing students to
visualize the evolution of the system.

The above logistic growth model, once implemented by the students, can be
adapted or modified in a number of ways. For example, by fitting the above
model to experimental data, or by extending the model to the generalized
logistic model, to include more than one competing population, or to allow for
time-varying carrying capacities (commonly used to capture acute hypoxia in
tumors). It is at this point that students are encouraged to be creative, and
use the model to either capture a system of interest or design their own
parameter sensitivity analysis.

## Example: Modeling the Zombie Apocalypse

The course is currently taught during the Fall semester, and in addition to
looking at more standard epidemiology models we also look at an epidemiology
model of the zombie apocalypse (known as the SIZR model). This model captures
the population dynamics of people susceptible to becoming zombies, the people
infected with a zombie virus, the number of zombies and the number of people
removed from the system (dead, but not undead) [14]. The equation that is used
to update the number of susceptible individuals is

$\frac{dS}{dt}=\Pi-\beta SZ-\delta S$

Where $S$ is the number of susceptible people, $t$ is time, $\Pi$ is the birth
rate of people, $\beta$ is the rate at which zombies infect susceptible
people, $Z$ is the number of zombies, and $\delta$ is the death rate of people
who are not zombies. The equation for updating the number of infected people
is

$\frac{dI}{dt}=\beta SZ-\rho I-\delta I$

Where $\rho$ is the rate at which the infected become zombies (inversely
related to the time it takes to transition). The number of zombies is updated
using the following equation

$\frac{dZ}{dt}=\rho I-\alpha SZ$

Where $\alpha$ is the rate at which susceptible individuals are fighting back
and destroying the zombies. Once a zombie is destroyed then it is removed from
the system. The number of removed individuals is of the form

$\frac{dR}{dt}=\delta S+\delta I+\alpha SZ$

Where the susceptible and infected that have died and the zombies killed by
susceptible individuals are included. The above equations can be solved using
simple finite difference approximations.

The Glowscript code for implementing this model is given below.

⬇

GlowScript 3.2 VPython

Pi = 0.

beta = 0.1

delta = 0

rho = 0.1

alpha = 0.01

dt = 0.1

t = 0.

The above code creates and stores variables for some of the constants in the
model. The initial populations are normalized as percentages and the zombies
initially comprise just 0.01% of the population.

⬇

S = 99.99

I = 0.

Z = 0.01

R = 0.0

In order to calculate the new values it is necessary to store the old values.

⬇

oldS = S

oldI = I

oldZ = Z

oldR = R

The syntax for plotting the data starts with declaring a variable that stores
the graph, and includes the title of the graph, the axis and size. Next the
different curves that will be plotted are defined, along with the color of the
lines and the labels that will be depicted in the key.

⬇

plot = graph(title='Zombie Apocalypse.', xtitle='time (days)',
ytitle='percentage',

fast=False, width=1000,height=500)

functS = gcurve(color=color.blue, width=2, label='Susceptible')

functI = gcurve(color=color.magenta, width=2, label='Infected')

functZ = gcurve(color=color.red, width=2, label='Zombies')

functR = gcurve(color=color.black, width=2, label='Removed')

We now turn our attention to iteratively solving the simulation. The current
implementation consists of 1000 time steps, which are contained in the
following while loop.

⬇

t=0

while t < 1000*dt:

rate(500)

t = t+dt

oldS = S

oldI = I

oldZ = Z

oldR = R

At each time step, the populations from the previous iteration are stored as
the old values and these old values are used in the finite difference
approximation when updating the new values.

⬇

S = oldS + dt*(Pi - beta*oldS*oldZ - delta*oldS)

I = oldI + dt*(beta*oldS*oldZ - rho*oldI - delta*oldI)

Z = oldZ + dt*(rho*oldI - alpha*oldS*oldZ)

R = oldR + dt*(delta*oldS + delta*oldI + alpha*oldS*oldZ)

The above code is simply the finite difference approximation of the equations
above. Once the values have been updated, all that remains is to update the
plot. This update occurs every iteration. This effectively creates an
animation of the graph evolving dynamically as the simulation progresses.

⬇

functS.plot( t, S)

functI.plot( t, I)

functZ.plot( t, Z)

functR.plot( t, R)

The final plot is shown below in Fig. 2. The concentrations for the
susceptible, infected, zombie, and removed populations are depicted as a
function of time. Given the values used here, the zombies very quickly overrun
the human population in true apocalyptic fashion.

Figure 2: A screenshot of the Glowscript implementation. The population of
susceptible, infected, zombified and removed individuals as a function of time
during a zombie apocalypse.

Once the basic model is implemented, students can modify the computer program
by either changing the equations or by changing the variables. However, rather
than randomly changing aspects of the numerical model, this system offers a
wonderful opportunity to highlight the connectedness between the mathematical
model (and the different terms and parameters in the model) and the system
that we are trying to capture. To this end, we talk about different zombie
movies, video games and television shows, and students are asked to write down
physical attributes that they believe describe the zombies in their favorite
media. Once we have these attributes we discuss how to capture these
attributes in the model. For example, the current implementation might not be
entirely unreasonable for _Z Nation_ or _World War Z_ , but would fail to
capture the dynamics of the zombies depicted in the _Night of the Living Dead_
, where the humans emerge victorious at the end. However, the timescale of
infection might be too long in this current implementation for capturing the
very fast dynamics of infection in _28 Days Later_. Through the variation of
parameters, or modification of the model, students gain a deeper understanding
of the different terms in the model and how they affect the dynamics of the
zombie apocalypse. This is an important illustration for students as to how
computer models essentially mimic reality, and how the model can be “fitted”
to meet our expectations.

## Example: Microtubule Dynamic Instabilities

In order to highlight some of the graphical capabilities of Glowscript, we
consider an example involving microtubule dynamic instabilities [4].
Microtubules, as part of the cytoskeleton, provide structure and facilitate
cell motility. However, at any given time a subset of microtubules will be
growing while others shrink. This allows the microtubules to rapidly explore
space, and reorganize if necessary (e.g., during mitosis). To capture
microtubule dynamics, various dynamic events in the system have to be
accounted for; the nucleation of microtubules, the polymerization of the
microtubule (which depends on tubulin and tau concentrations), the
depolymerization of microtubules and the hydrolysis of GTP-tubulin in the
microtubules to GDP-tubulin (drastically affecting polymerization kinetics).
Microtubules are nucleated at a given nucleation site with a rate of
$k_{new}=0.16\,\text{s}^{-1}$ [16]. The GTP-tubulin is converted in the wall
of microtubules into GDP-tubulin with a rate of
$k_{hydro}=0.04\,\text{s}^{-1}$. The polymerization kinetics of a microtubule
depends on the local tubulin and tau concentrations and, if the end of the
microtubule consists of a GTP-tubulin cap, is taken to be of the form

$k_{poly}=k_{poly,0}\phi_{tub}+k_{poly,\tau}\phi_{tub}\phi_{tau}$

where $k_{poly,0}=0.004\,\text{s}^{-1}\,\mu\text{M}^{-1}$ and
$k_{poly,\tau}=0.01\,\text{s}^{-1}\,\mu\text{M}^{-2}$ [5]. However, if the
GTP-tubulin at the end of the microtubule hydrolyzes to GDP-tubulin, then the
rate of polymerization is taken to be zero. Similarly, we assume the
depolymerization rate, $k_{depoly}$, is 0.7 $\,\text{s}^{-1}$ if the end of
the microtubule consists of GDP-tubulin, but is zero if it is a cap of GTP-
tubulin [6]. Once we have the rates associated with all of the possible events
within the system we define an average time over which at least one of these
events is likely to occur. In particular, the time step is taken to be the
inverse of the sum over all rates.

$\Delta t=\frac{1}{\sum_{i}^{N}k_{i}}$

This time step gives a timescale to the simulation. Assuming an event occurs
somewhere in the system, during this time interval, the probability of a
particular event occurring is the rate associated with this event relative to
the total rate (associated with all events)

$P_{i}=\frac{k_{i}}{\sum_{j}^{N}k_{j}}$

where, again, the summation is over all events. An event is chosen to occur
from the cumulative distribution function associated with these probabilities
and, therefore, the selection of an event takes into consideration the correct
probability weightings. The cumulative distribution function is defined as

$c_{i}=\frac{\sum_{j=0}^{i}k_{j}}{\sum_{j=0}^{N}k_{j}}$

In order to select an event which will occur during this time step, a random
number between 0 and 1 is chosen (RND [0,1]) and the ith event is found for
which $c_{i}<\text{RND [0,1]}<c_{i+1}$. In this manner, the polymerization
kinetics of the microtubule dynamics occur stochastically but with the correct
probability weightings.

Not only do the polymerization rates depend on the local concentrations, but
the concentrations are increased or decreased in response to the
depolymerization or polymerization of a microtubule. In particular, the
concentration of tubulin associated with $0.5\,\mu\text{m}$ of a microtubule
(which we add or subtract in the current model) is $11.6\,\mu\text{M}$ and the
amount of tau is taken to be $0.5\,\mu\text{M}$. Therefore, the average
concentration of tubulin and tau in the cell varies depending on the number
and average length of the microtubules. It is worth noting that the local
concentrations could also play a role in microtubule dynamics and, for
example, a depolymerizing microtubule could leave behind a trail of higher
tubulin and tau concentrations which might encourage polymerization along this
trail. However, capturing the diffusion of concentrations throughout this
three-dimensional system is beyond the constraints of an in-class
implementation.

One of the wonderful aspects of Glowscript exemplified in this model is the
graphics capabilities. In particular, we use curves in Glowscript. This is a
function that allows an array of points to characterize a three-dimensional
curve through space. Furthermore, each point along the curve can be assigned
different thicknesses or colors. For example, here we change the colors of the
curves that depict the microtubules in the simulations to show the hydrolysis
of GTP-tubulin in the microtubules to GDP-tubulin. Glowscript can include a
number of shapes, such as spheres, rectangles and cones (just to name a few).

Figure 3: Screenshot of the Glowscript implementation showing snapshots of the
morphology of the microtubules that are displayed as the simulation is
running. Systems with a) normal tau concentrations, and b) tau concentrations
reduced to 50%, are contrasted.

The output of simulations are depicted in Fig. 3. The microtubules are colored
blue if they consist of GDP-tubulin and red if they consist of GTP-tubulin.
These simulations took a few minutes to run within a Chrome browser on a
typical computer. Fig. 3a shows the morphology of the microtubules for a
system with normal levels of tau (around $1\,\mu\text{M}$). While Fig. 3b
shows the morphology of the microtubules for a system with reduced tau
concentrations ($0.5\,\mu\text{M}$). This might mimic a tauopathy where the
tau is deposited or aggregated as tau tangles (such as Alzheimer’s disease).
Clearly a reduction of tau concentrations results in a severe inhibition of
microtubule growth. Students can explore this model further by looking at the
statistics of the stochastic model. For example, the size distribution of
microtubules or the ratio of GDP-tubulin to GTP-tubulin within the
microtubules as a function of time. Students can also vary the parameters in
the model. For example, as shown here, the initial concentrations of the tau
available to the microtubules has a profound effect on the system. therefore,
students can investigate the onset of different neurodegenerative diseases
through the variations of the parameters in the model [15].

## Example: Drug Encapsulation in Polymer Nanoparticle

Premedical students, in particular, are often interested in numerical models
that capture therapeutic systems. As an example, we include a model of drug
delivery. Polymer nanoparticles can be used to deliver chemotherapeutic drugs
to a tumor, and selectively release the drug at the tumor site [3, 1]. To
introduce the students to these systems we capture the polymer nanoparticle
and the diffusion of the encapsulated drug from the nanoparticle. The polymer
can be constructed from acid-labile bonds that disintegrate in the more acidic
environment of the tumor. For computational simplicity (we want this to be
solved within a couple of minutes during class) the system is taken to be
spherically symmetric and modeled in one-dimension.

In spherical coordinates, the diffusion equation is

$\frac{\partial\phi}{\partial t}=\frac{1}{r^{2}}\frac{\partial}{\partial
r}\left(r^{2}D\frac{\partial\phi}{\partial r}\right)$

where $\phi$ is the concentration, $t$ is time, $r$ is the radial distance,
and $D$ is the diffusion coefficient. The finite difference approximation is
of the form

$\frac{d\phi}{dt}=\begin{cases}6\dfrac{D}{(\Delta
r)^{2}}(\phi(1)-\phi(0))&\text{if }i=0\vspace{5pt}\\\ 0,&\text{if
}i=N\vspace{5pt}\\\ \dfrac{D}{(\Delta
r)^{2}i}\left[(i+1)\phi(i+1)-2i\phi(i)+(i-1)\phi(i-1)\right],&\text{otherwise}\end{cases}$

The diffusion of the drug concentration is

$\frac{d\phi_{d}(i)}{dt}=\begin{cases}6\dfrac{D_{d}}{(\Delta
r)^{2}}(\phi_{d}(1)-\phi_{d}(0))&\text{if }i=0\vspace{5pt}\\\ 0,&\text{if
}i=N\vspace{5pt}\\\ \dfrac{D_{d}}{(\Delta
r)^{2}i}\left[(i+1)\phi_{d}(i+1)-2i\phi_{d}(i)+(i-1)\phi_{d}(i-1)\right],&\text{otherwise}\end{cases}$

where the subscript $d$ indicates quantities associated with the drug. The
diffusion coefficient is given by Darcy’s law which modifies the diffusion
coefficient to account for permeability (here taken to be the space that is
not solid polymer).

$D_{d}(i)=D_{d0}(1-\phi_{s}(i))$

In terms of the nanoparticle we have the solid polymer concentration,
$\phi_{ps}$, and the concentration of diffusing polymers that have
disintegrated from the nanoparticle, $\phi_{p}$. The rate at which the solid
particle is broken down will depend on the concentration of the blood, and the
exposure of the polymer to the more acidic environment near a tumor. The
concentration of blood is simply taken to be the volume of the system not
occupied by polymer or drug.

$\phi_{b}(i)=1-\phi_{d}(i)-\phi_{p}(i)-\phi_{s}(i)$

The degradation of the solid nanoparticle is then captured using

$\frac{d\phi_{s}(i)}{dt}=-k\phi_{s}(i)\phi_{b}(i)$

where $k$ is a constant, and the degradation depends on the concentration of
solid nanoparticle that is being degraded and the concentration of more acidic
blood that is responsible for the degradation. As the solid degrades we see an
increase in the concentration of the diffusing polymer, whose diffusion can
also be captured in spherical coordinates.

$\frac{d\phi_{p}(i)}{dt}=\begin{cases}6\dfrac{D_{p}}{(\Delta
r)^{2}}(\phi_{p}(1)-\phi_{p}(0))+k\phi_{s}(0)\phi_{b}(0)&\text{if
}i=0\vspace{5pt}\\\ 0,&\text{if }i=N\vspace{5pt}\\\ \dfrac{D_{p}}{(\Delta
r)^{2}i}\left[(i+1)\phi_{p}(i+1)-2i\phi_{p}(i)+(i-1)\phi_{p}(i-1)\right]+k\phi_{s}(i)\phi_{b}(i),&\text{otherwise}\end{cases}$

where the diffusion coefficient is again given by Darcy’s law

$D_{d}(i)=D_{d0}(1-\phi_{s}(i))$

Note that the degradation of the polymer nanoparticle could result in a
polydispersity in diffusing polymer fragments, and replacing this with a
single diffusing polymer concentration is an approximation.

The Glowscript implementation (along with all the code discussed here) is
provided as supplementary information. The built-in graphing capabilities
allow the concentrations to be plotted as a function of radial position, and
for these plots to animate as the simulation progresses. Figure 4 depicts an
output of a simulation. The solid nanoparticle dissolves and the polymer
fragments are created more at the particle interface and diffuse out. The
reduction in solid nanoparticle concentration increases permeability and the
chemotherapeutic drug is able to diffuse more readily from the nanoparticle.

Possible modifications and investigations that students might make to this
base model could be the calculation of drug release rate, the incorporation of
different sized polymer fragments (with different diffusion coefficients), the
investigation of the effects of the degradation rate on drug release rates,
the encapsulation of a drug in a hollow nanocapsule, or (in a more advanced
modification) the incorporation of hydrophobic drugs or polymers. Students,
therefore, not only implement and use a model of drug delivery from a polymer
nanoparticle, but are required to modify and explore this model further.
Students, therefore, exhibit creativity and implement scientific methodologies
that reflect research practices.

Figure 4: A screenshot of the Glowscript implementation. The concentrations of
solid polymer, diffusing polymer and chemotherapeutic drug from a
nanoparticle.

## Conclusions

Glowscript can be used to introduce computer programming, mathematical
modeling and numerical simulations to biology undergraduate students. Examples
of different models that have been implemented in a computational biology
course taught at Robert Morris University are highlighted; however, the
utility of Glowscript extends beyond the applications considered here. The
ability of Glowscript to display graphical objects, including a dynamic
plotting functionality, makes this an ideal programming environment for
introducing computer simulations to science students. One aspect of this
course that we have observed over the years is that a number of biology
students have been motivated to implement computer models outside of this
class, especially in senior research projects. We hope that others will read
this article and be inspired to implement Glowscript within their classrooms,
especially in the biological sciences where students are traditionally less
exposed to numerical modeling.

## References

  * [1] Gavin A Buxton.  Simulating the co-encapsulation of drugs in a “smart” core-shell-shell polymer nanoparticle.  The European Physical Journal E, 37(3):1–7, 2014.
  * [2] Gavin A Buxton.  Mathematical modelling and computer simulations in undergraduate biology education.  Spreadsheets in Education, 10(3):4673, 2018.
  * [3] Gavin A Buxton and Nigel Clarke.  Acid-labile core-layer–shell polymer nanoparticles: modeling particle disintegration and subsequent drug release.  Soft Matter, 8(34):8891–8895, 2012.
  * [4] Gavin A Buxton, Sandra L Siedlak, George Perry, and Mark A Smith.  Mathematical modeling of microtubule dynamics: insights into physiology and disease.  Progress in neurobiology, 92(4):478–483, 2010.
  * [5] David N Drechsel, AA Hyman, Melanie H Cobb, and MW Kirschner.  Modulation of the dynamic instability of tubulin assembly by the microtubule-associated protein tau.  Molecular biology of the cell, 3(10):1141–1154, 1992.
  * [6] Harold P Erickson and E Timothy O’Brien.  Microtubule dynamic instability and gtp hydrolysis.  Annual review of biophysics and biomolecular structure, 21(1):145–166, 1992.
  * [7] Tim W Fawcett and Andrew D Higginson.  Heavy use of equations impedes communication among biologists.  Proceedings of the National Academy of Sciences, 109(29):11735–11739, 2012.
  * [8] Jason Feser, Helen Vasaly, and Jose Herrera.  On the edge of mathematics and biology integration: improving quantitative skills in undergraduate biology education.  CBE—Life Sciences Education, 12(2):124–128, 2013.
  * [9] Linda Grandell, Mia Peltomäki, Ralph-Johan Back, and Tapio Salakoski.  Why complicate things? introducing programming in high school using python.  In Proceedings of the 8th Australasian Conference on Computing Education-Volume 52, pages 71–80, 2006.
  * [10] Louis J Gross.  Education for a biocomplex future.  Science, 288(5467):807–807, 2000.
  * [11] Fionnuala Johnson, Stephen McQuistin, and John O’Donnell.  Analysis of student misconceptions using python as an introductory programming language.  In Proceedings of the 4th Conference on Computing Education Practice 2020, pages 1–4, 2020.
  * [12] Eva Mészárosová.  Is python an appropriate programming language for teaching programming in secondary schools.  International Journal of Information and Communication Technologies in Education, 4(2):5–14, 2015.
  * [13] Bradley N Miller and David L Ranum.  Teaching an introductory computer science sequence with python.  In Proc. 38th Midwest Instructional and Computing Symposium, Eau Claire, Wisconsin, USA, 2005.
  * [14] Philip Munz, Ioan Hudea, Joe Imad, and Robert J Smith?  When zombies attack!: mathematical modelling of an outbreak of zombie infection.  Infectious disease modelling research progress, 4:133–150, 2009\.
  * [15] Dulal Panda, Jonathan C Samuel, Michelle Massie, Stuart C Feinstein, and Leslie Wilson.  Differential regulation of microtubule dynamics by three-and four-repeat tau: implications for the onset of neurodegenerative disease.  Proceedings of the National Academy of Sciences, 100(16):9548–9553, 2003.
  * [16] Michelle Piehl, U Serdar Tulu, Pat Wadsworth, and Lynne Cassimeris.  Centrosome maturation: measurement of microtubule nucleation throughout the cell cycle by using gfp-tagged eb1.  Proceedings of the National Academy of Sciences, 101(6):1584–1588, 2004.
  * [17] J Vandermeer.  How populations grow: the exponential and logistic equations.  Nature Education Knowledge, 3(10):15, 2010.