1 Department of Econometrics and Business
Statistics, Monash University, Melbourne,
Victoria, Australia
2 Institute of
Statistics, University of Natural Resources and Life
Sciences, Vienna, Austria
3 Guangdong Provincial Key
Laboratory of Interdisciplinary Research and Application for Data
Science, BNU-HKBU United International College, Zhuhai, Guangdong, China

Indexes are useful for summarizing multivariate information into single
metrics for monitoring, communicating, and decision-making. While most
work has focused on defining new indexes for specific purposes, more
attention needs to be directed towards making it possible to understand
index behavior in different data conditions, and to determine how their
structure affects their values and variation in values. Here we discuss
a modular data pipeline recommendation to assemble indexes. It is
universally applicable to index computation and allows investigation of
index behavior as part of the development procedure. One can compute
indexes with different parameter choices, adjust steps in the index
definition by adding, removing, and swapping them to experiment with
various index designs, calculate uncertainty measures, and assess
indexes' robustness. The paper presents three examples to illustrate the
pipeline framework usage: comparison of two different indexes designed
to monitor the spatio-temporal distribution of drought in Queensland,
Australia; the effect of dimension reduction choices on the Global
Gender Gap Index (GGGI) on countries' ranking; and how to calculate
bootstrap confidence intervals for the Standardized Precipitation Index
(SPI). The methods are supported by a new R package, called

indexesdata pipelinesoftware designuncertaintydecision-making

Shaded[frame hidden, breakable, enhanced, boxrule=0pt, interior hidden, borderline west=3pt0ptshadecolor, sharp corners]

§ INTRODUCTION

Indexes are commonly used to combine and summarize different sources of
information into a single number for monitoring, communicating, and
decision-making. They serve as critical tools across the natural and
social sciences. Examples include the Air Quality Index, El
Niño-Southern Oscillation Index, Consumer Price Index, QS University
Rankings, and the Human Development Index. In environmental science,
climate indexes are produced by major monitoring centers, like the
United States Drought Monitor and National Oceanic and Atmospheric
Administration, to facilitate agricultural planning and early detection
of natural disasters. In economics, indexes provide insight into market
trends through combining prices of a basket of goods and services. In
social sciences, indexes are used to monitor human development, gender
equity, or university quality. The problem is that every index is
developed in its own unique way, by different researchers or
organizations, and often indexes designed for the same purpose cannot
easily be compared.

To construct an index, experts typically start by defining a concept of
interest that requires measurement. This concept often lacks a direct
measurable attribute or can only be measured as a composite of various
processes, yet it holds social and public significance. To create an
index, once the underlying processes involved are identified, relevant
and available variables are then defined, collected, and combined using
statistical methods into an index that aims to measure the process of
interest. The construction process is often not straightforward, and
decisions need to be made, such as the selection of variables to be
included, which might depend on data availability and the statistical
definition of the index to be used, among others. For instance, the
indexes constructed from a linear combination of variables require to
decide the weight assigned to each variable. Some indexes have a spatial
and/or temporal component, and variables can be aggregated to different
spatial resolutions and temporal scales, leading to various indexes for
different monitoring purposes. Hence, all these decisions can result in
different index values and have different practical implications.

To be able to test different decision choices systematically for an
index, the index needs to be broken down into its fundamental building
blocks to analyze the contribution and effect of each component. We call
this process of breaking the index construction into different steps the
index pipeline. Such decomposition of index components provides the
means to standardize index construction via a pipeline and offers
benefits for comparing among indexes, calculating index uncertainty, and
assessing index robustness.

In this work, we provide statistical and computational methods for
developing a data pipeline framework to construct and customize indexes
using data. The proposed pipeline comprises various modules, including
temporal and spatial aggregation, variable transformation and
combination, distribution fitting, benchmark setting, and index
communication. Given the decisions analysts need to make when combining
multivariate data into indexes, the proposed pipeline enables the
evaluation of how the specific choice can affect the index, as well as
how the index may appear under alternative options. Furthermore,
uncertainty calculation can also flow through the pipeline, providing
the index with confidence measures.

The rest of the paper is structured as follows.
Section <ref> provides background about the development of
indexes. Section <ref> reviews the tidy framework in R and how
index construction can benefit from such a framework. The details of the
pipeline modules are presented in Section <ref>.
Section <ref> explains the design of the
package that implements the modules. Examples are given in
Section <ref> to illustrate three use cases of the

§ BACKGROUND TO INDEX DEVELOPMENT

There are many documents providing advice on how to construct indexes
for different fields, and review articles describing the range of
available indexes for specific purposes. The OECD handbook (OECD,
European Union, and Joint Research Centre - European Commission 2008)
provides a comprehensive guide for computing socio-economic composite
indexes, with detailed steps and recommendations. The drought index
handbook (Svoboda, Fuchs, et al. 2016) provides details of various
drought indexes and recommendations from the World Meteorology
Organization. Zargar et al. (2011), Hao and Singh (2015) and Alahacoon
and Edirisinghe (2022) are review papers describing the range of
possible drought indexes.

There is also some attention being given to the diagnosis of indexes,
and incorporation of uncertainty. Jones and Andrey (2007) investigates
the methodological choices made in the development of indexes for
assessing vulnerable neighborhoods. Saisana, Saltelli, and Tarantola
(2005) describes incorporating uncertainty estimates and conducting
sensitivity analysis on composite indexes. Tate (2012) and Tate (2013),
similarly, make a comparative assessment of social vulnerability indexes
based on uncertainty estimation and sensitivity analysis. Laimighofer
and Laaha (2022) studies five uncertainty sources (record length,
observation period, distribution choice, parameter estimation method,
and GOF-test) of drought indexes.

There are also a few R packages supporting index calculation. The
package (Vicente-Serrano, Beguería, and López-Moreno 2010)
computes two drought indexes. The  package (Martin 2023)
computes price indexes, and the  package (Grenié
and Gruson 2023) computes functional diversity indexes for ecological
study. The package  (Becker et al. 2022) is more
ambitious, making a start on following the broader guidelines in the
OECD handbook to construct, analyze, and visualize composite indexes.

From reviewing this literature, and in the process of developing methods
for making it easier to work with multivariate spatio-temporal data, it
seems possible to think about indexes in a more organised, cohesive and
standard manner. Actually, the area could benefit from a tidy

§ TIDY FRAMEWORK

The tidy framework consists of two key components: tidy data and tidy
tools. The concept of tidy data (Wickham 2014) prescribes specific rules
for organizing data in an analysis, with observations as rows, variables
as columns, and types of observational units as tables. Tidy tools, on
the other hand, are concatenated in a sequence through which the tidy
data flows, creating a pipeline for data processing and modeling. These
pipelines are data-centric, meaning all the tidy tools or functions take
a tidy data object as input and return a processed tidy data object,
directly ready for the next operations to be applied. Also, the pipeline
approach corresponds to the modular programming practice, which breaks
down complex problems into smaller and more manageable pieces, as
opposed to a monolithic design, where all the steps are predetermined
and integrated into a single piece. The flexibility provided by the
modularity makes it easier to modify certain steps in the pipeline and
to maintain and extend the code base.

Examples of using a pipeline approach for data analysis can be traced
back to the interactive graphics literature, including A. Buja et al.
(1988); Sutherland et al. (2000); Wickham et al. (2009); Xie, Hofmann,
and Cheng (2014). Wickham et al. (2009) argue that whether made explicit
or not, a pipeline has to be presented in every graphics program, and
making them explicit is beneficial for understanding the implementation
and comparing between different graphic systems. While this comment is
made in the context of interactive graphics programs, it is also
applicable generally to any data analysis workflow. More recently, the
tidyverse suite (Wickham et al. 2019) takes the pipeline approach for
general-purpose data wrangling and has gained popularity within the R
community. The pipeline-style code can be directly read as a series of
operations applied successively on tidy data objects, offering a method
to document the data wrangling process with all the computational
details for reproducibility.

Since the success of tidyverse, more packages have been developed to
analyze data using the tidy framework for domains specific applications,
a noticeable example of which is  for building
machine learning models (Kuhn and Wickham 2020). To create a tidy
workflow tailored to a specific domain, developers first need to
identify the fundamental building blocks to create a workflow. These
components are then implemented as modules, which can be combined to
form the pipeline. For example, in supervised machine learning models,
steps such as data splitting, model training, and model evaluation are
commonly used in most workflow. In the , these steps
are correspondingly implemented as package ,
, and , agnostic to the specific model
chosen. The uniform interface in tidymodels frees analysts from
recalling model-specific syntax for performing the same operation across
different models, increasing the efficiency to work with different
models simultaneously.

For constructing indexes, the pipeline approach adopts explicit and
standalone modules that can be assembled in different ways. Index
developers can choose the appropriate modules and arrange them
accordingly to generate the data pipeline that is needed for their
purpose. The pipeline approach provides many advantages:

* makes the computation more transparent, and thus more easily debugged.

* allows for rapidly processing new data to check how different
features, like outliers, might affect the index value.

* provides the capacity to measure uncertainty by computing confidence
intervals from multiple samples as generated by bootstrapping to
original data.

* enables systematic comparison of surrogate indexes designed to measure
the same phenomenon.

* it may even be possible to automate diagrammatic explanations and
documentation of the index.

Adoption of this pipeline approach would provide uniformity to the field
of index development, research, and application to improve
comparability, reproducibility, and communication.

§ DETAILS OF THE INDEX PIPELINE

In constructing various indexes, the primary aim is to transform the
data, often multivariate, into a univariate index. Spatial and temporal
considerations are also factored into the process when observational
units and time periods are not independent. However, despite the
variations in contextual information for indexes in different fields,
the underlying statistical methodology remains consistent across diverse
domains. Each index can be represented as a series of modular
statistical operations on the data. This allows us to decompose the
index construction process into a unified pipeline workflow with a
standardized set of data processing steps to be applied across different

An overview of the pipeline is presented in
Figure <ref>, illustrating the nine available modules
designed to obtain the index from the data. These modules include
operations for temporal and spatial aggregation, variable transformation
and combination, distribution fitting, benchmark setting, and index
communication. Analysts have the flexibility to construct indexes by
connecting modules according to their preferences.

Diagram of pipeline modules for index
construction. The highlighted path illustrates one possible construction
using the dimension reduction and simplification module.

Now, we introduce the notation used for describing pipeline modules.
Consider a multivariate spatio-temporal process,

\begin{equation}
\mathbf{x}(s;t) = \{x_1(s;t), x_2(s;t), \cdots, x_p(s;t)\} \qquad s \in D_s \subseteq \mathbb{R}^m, t \in D_t \subseteq \mathbb{R}^n
\end{equation}

* \(x_j(s, t)\) represents a variable of interest for example
precipitation, \(j = 1, \cdots, p\) and

* \(s\) represents the geographic locations in the space
\(D_s \subseteq \mathbb{R}^m\). Examples of geographic locations
include a collection of countries, longitude and latitude coordinates
or regions of interest and,

* \(t\) denotes the temporal order in \(D_t \subseteq \mathbb{R}^n\).
For instance, time measurements could be recorded hourly, yearly,
monthly, quarterly, or by season.

In what follows when geographic or temporal components of the
\(x_j(s,t)\) process are fixed we will be using suffix notation. For
example, \(x_{sj}(t)\) represents the data for a fixed location \(s\) as
a function of time \(t\). While \(x_{tj}(s)\) denotes the spatial
varying process for a fixed \(t\). An overview of the notation for
pipeline input, operation, and output is present in
Table <ref>.

>p(- 6) * 0.2188
>p(- 6) * 0.1406
>p(- 6) * 0.2344
>p(- 6) * 0.4062@
Summary of the notation for input,
operation, and output of each pipeline module.

Temporal processing     \(x_{sj}(t)\)     \(f[x_{sj}(t)]\)
\(x^{\text{Temp}}_{sj}(t^\prime) \quad t^\prime \in D_{t^\prime}\)

Spatial processing     \(x_{tj}(s)\)     \(g[x_{tj}(s)]\)
\(x^{\text{Spat}}_{tj}(s^\prime) \quad s^\prime \in D_{s^\prime}\)

Variable transformation     \(x_{j}(s; t)\)     \(T[x_j(s;t)]\)
\(x^{\text{Trans}}_j(s;t)\)

Scaling     \(x_j(s; t)\)     \([x_j(s;t) - \alpha]/\gamma\)
\(x^{\text{Scale}}_j(s;t)\)

Dimension reduction     \(\mathbf{x}(s;t)\)     \(h[\mathbf{x}(s;t)]\)
\(\mathbf{y}(s;t) \quad \mathbf{y} \subseteq \mathbb{R}^d, d < p\)

Distribution fit     \(x_j(s; t)\)     \(F[x_j(s;t)]\)
\(P_j(s;t) \quad P(.) \in [0, 1]\)

Normalising     \(x_j(s; t)\)     \(\Phi^{-1}[x_j(s; t)]\)
\(z_j(s; t)\)

Benchmarking     \(x_j(s; t)\)     \(u[x_j(s;t)]\)     \(b_j(s;t)\)

Simplification     \(x_j(s; t)\)     \(v[x_j(s;t)]\)
\(A_j(s;t) \in \{a_1, a_2, \cdots, a_z\}\)

§.§ Temporal processing

The temporal processing module takes as input argument a single variable
\(x_{sj}(t)\) at location \(s\) as a function of time. In this step the
original time series can be transformed or summarized into a new one via
time aggregation, the transformation is represented by the function
\(f\), \(x^{\text{Temp}}_{sj}(t^\prime) = f[x_{sj}(t)]\) where
\(t^\prime\) refers to the new temporal resolution after aggregation. An
example of temporal processing done in the computation of the
Standardized Precipitation Index (SPI) (McKee et al. 1993), consists of
summing the monthly precipitation series over a rolling time window of
size \(k\). That is also known as the time scale. For SPI, the choice of
the time scale \(k\) is used to control the accumulation period for the
water deficit, enabling the assessment of drought severity across
various types (meteorological, agricultural, and hydrological).

§.§ Spatial processing

The spatial processing module takes a single variable with a fixed
temporal dimension, \(x_{tj}(s)\), as input. This step transforms the
variable from the original spatial dimension \(s\) into the new
dimension \(s^\prime \in D_{s^\prime}\) through
\(x^{\text{Spat}}_{tj}(s^\prime) = g[x_{tj}(s)]\). The change of spatial
dimension allows for the alignment of variables collected from different
measurements, such as in-situ stations and satellite imagery, or
originating from different resolutions. This also includes the
aggregation of variables into different levels, such as city, state, and
country scales.

§.§ Variable transformation

Variable transformation takes the input of a single variable
\(x_j(s;t)\) and reshapes its distribution using the function \(T[.]\)
to produce \(x^{\text{Trans}}_{j}(s;t)\). When a variable has a skewed
distribution, transformations such as log, square root, or cubic root
can adjust the distribution towards normality. For example, in the Human
Development Index (HDI), a logarithmic transformation is applied to the
variable Gross National Income per capita (GNI), to reduce its impact on
HDI, particularly for countries with high GNI values.

§.§ Scaling

Unlike variable transformation, scaling maintains the distributional
shape of the variable. It includes techniques such as centering, z-score
standardization, and min-max standardization and can be expressed as
\([x_{j}(s;t) - \alpha]/\gamma\). In the Human Development Index (HDI),
the three dimensions (health, education, and economy) are converted into
the same scale (0-1) using min-max standardization.

Although the scaling might be considered to be a transformation, we have
elected to make it a separate module because it is neater. Scaling
simply changes the numbers in the data not the shape of a variable.
Transformation will most likely change the shape, and is usually

Comparison of the module
scaling (green) and variable transformation (orange). While both modules
change the variable range, scaling maintains the same distributional
shape, which is not the case with variable transformation.

§.§ Dimension reduction

Dimension reduction takes the multivariate information
\(\mathbf{x}(s;t)\), where \(\mathbf{x} \subseteq \mathbb{R}^p\), or a
subset of variables in \(\mathbf{x}(s;t)\), as the input. It summarises
the high-dimensional information into a lower-dimension representation
\(\mathbf{y}(s;t)\), where \(y \subseteq \mathbb{R}^d\) and \(d < p\),
as the output. The transformation can be based on domain-specific
knowledge, originating from theories describing the underlying physical
processes, or guided by statistical methods. For example, the
Standardized Precipitation-Evapotranspiration Index (SPEI) (Beguería and
Vicente-Serrano 2017) calculates the difference \(D\) between
precipitation (\(P\)) and potential evapotranspiration (\(\text{PET}\)),
using a water balance model (\(D = P - \text{PET}\)). This is the only
step that differs from the Standardized Precipitation Index (SPI), and
can be considered to be a dimension reduction using a particular linear

Linear combinations of variables are commonly used to reduce the
dimension in statistical methodology, and chosen using a method like
principal component analysis (PCA) (Hotelling 1933) or linear
discriminant analysis (Ronald A. Fisher 1936), preparing contrasts to
test particular elements in analysis of variance (Ronald Aylmer Fisher
1970), or hand-crafted by a content-area expert. Linear combinations
also form the basis for visualizing multivariate data, in methods such
as tours (Wickham et al. 2011). This dimension reduction method can
accommodate linear combinations as provided by any method, and hence is
linear by design. The transformation module provides variable-wise
non-linear transformation.

§.§ Distribution fit

Distribution fit applies the Cumulative Distribution Function (CDF)
\(F\) of a distribution on the variable \(x_j(s; t)\) to obtain the
probability values \(P_j(s;t) \in [0, 1]\). In SPEI, many distributions,
including log-logistic, Pearson III, lognormal, and general extreme
distribution, are candidates for the aggregated series. Different
fitting methods and different goodness of fit tests may be used to
compare the distribution choice on the index value. This could be
considered to be a variable transformation because it is usually
conducted separately for each variable. However, very occasionally a fit
is conducted on two or more variables simultaneously. For this reason,
and because it usually is applied later in the pipeline it is neater to
make this a separate module.

§.§ Normalising

Normalizing applies the inverse normal CDF \(\Phi^{-1}\) on the input
data to obtain the normal density \(z_{j}(s;t)\). Normalizing can
sometimes be confused with the scaling or variable transformation
module, which does not involve using a normal distribution to transform
the variable. It is arguably whether normalizing and distribution fit
should be combined or separated into two modules. A decision has been
made to separate them into two modules given the different types of
output each module presents (probability values for distribution fit and
normal density values for normalizing).

§.§ Benchmarking

Benchmark sets a value \(b_j(s,t)\) for comparing against the original
variable \(x_j(s;t)\). This benchmark can be a fixed value consistently
across space and time, perhaps extracted from expert knowledge or
determined by the data through the function \(u[x_j(s;t)]\). Once a
benchmark is set, observations can be highlighted for adjustments in
other modules or can serve as targets for monitoring and planning.

§.§ Simplification

Simplification takes a continuous variable \(x_j(s;t)\) and categorises
it into a discrete set \(A_j(s;t) \in \{a_1, a_2, \cdots, a_z\}\)
through a piecewise constant function,

\begin{equation}
v[x_i(s;t)] =
\begin{cases}
a_0 & C_1 \leq x^i(s; t) < C_0 \\
a_1 & C_2 \leq x^i(s; t) < C_1 \\
a_2 & C_3 \leq x^i(s; t) < C_2 \\
\cdots \\
a_z & C_z \leq x^i(s; t)
\end{cases}
\end{equation}

This is typically used at the end of the index pipeline to simplify the
index to communicate to the public the severity of the concept of
interest measured by the index. An example of simplification is to map
the calculated SPI to four categories: mild, moderate, severe, and
extreme drought.

§ SOFTWARE DESIGN

The R package  implements a proof-of-concept of the
index pipeline modules described in Section <ref>. These
modules compute an index in a sequential manner, as shown below:

DATA |
module1(...) |
module2(...) |
module3(...) |

Each module offers a variety of alternatives, each represented by a
distinct function. For example, within the
module, three methods are available:
, , and
and they can be used as:

dimension_reduction(V1 = aggregate_linear(...))
dimension_reduction(V2 = aggregate_geometrical(...))
dimension_reduction(V3 = manual_input(...))

Each method can be independently evaluated as a recipe, for example,

manual_input(~x1 + x2)

takes a formula to combine the variables  and  and

This recipe will then be evaluated in the pipeline module with data to
obtain numerical results. The package also offers wrapper functions that
combine multiple steps for specific indexes. For instance, the
function bundles three steps (temporal aggregation,
distribution fit, and normalizing) into a single command, simplifying
the syntax for computation. Analysts are also encouraged to create
customized indexes from existing modules.

idx_spi - function(...){
DATA |
temporal_aggregate(...) |

The accompanied package, tidyindex, is not intended to offer an
exhaustive implementation for all indexes across every domains. Instead,
it provides a realization of the pipeline framework proposed in the
paper. When adopting the pipeline approach to construct indexes,
analysts may consider developing software that can be readily deployed
in the cloud for production purposes.

§ EXAMPLES

This section uses the example of drought and social indexes to show the
analysis made possible with the index pipeline. The drought index
example computes two indexes (SPI and SPEI) with various time scales and
distributions simultaneously using the pipeline framework to understand
the flood and drought events in Queensland. The second example focuses
on the dimension reduction step in the Global Gender Gap Index to
explore how the changes in linear combination weights affect the index
values and country rankings.

§.§ Every distribution, every scale, every index all at

The state of Queensland in Australia frequently experiences natural
disaster events such as flood and drought, which can significantly
impact its agricultural industry. This example uses daily data from
Global Historical Climatology Network Daily (GHCND), aggregated into
monthly precipitation, to compute two drought indexes – SPI and SPEI –
at various time scales and fitted distributions, for 29 stations in the
state of Queensland in Australia, spanning from January 1990 to April
2022. This example showcases the basic calculation of indexes with
different parameter specifications within the pipeline framework. The
dataset used in this example is available in the
package as  and blow prints the first few rows of the

Index pipeline for two drought indexes: the
Standardized Precipitation Index (SPI) and the Standardized
Precipitation-Evapotranspiration Index (SPEI). Both indexes share
similar construction steps with SPEI having two additional steps
(variable transformation and dimension reduction) to convert temperature
into evapotranspiration and combine it with the precipitation series.

Figure <ref> illustrates the pipeline steps of the two indexes.
The two indexes are similar with the difference that SPEI involves two
additional steps – variable transformation and dimension reduction –
prior to temporal processing. As introduced in
Section <ref>, wrapper functions are available for both
indexes as  and , which allows
for the specification of different time scales and distributions for
fitting the aggregated series. In , multiple indexes
can be calculated collectively using the function
. Both SPI and SPEI are calculated across
four time scales (6, 12, 24, and 36 months). The SPEI is fitted with two
distributions (log-logistic and general extreme value distribution) and
the gamma distribution is used for SPI:

.scale - c(6, 12, 24, 36)
idx - queensland %%
mutate(month = lubridate::month(ym)) |
init(id = id, time = ym, group = month) |
spei = idx_spei(
.tavg = tavg, .lat = lat,
.scale = .scale, .dist = list(dist_gev(), dist_glo())),
spi = idx_spi(.scale = .scale)

The output contains the original data, index values (),
parameters used (, , and ),
and all the intermediate variables (, , and
). This data can be visualized to investigate the
spatio-temporal distribution of the drought or flood events, as well as
the response of index values to different time scales and distribution
parameters at specific single locations.
Figure <ref> and Figure <ref>
exemplify two possibilities. Figure <ref> presents
the spatial distribution of SPI during two periods: October 2010 to
March 2011 for the 2010/11 Queensland flood and October 2019 to March
2020 for the 2019 Australia drought, which contributes to the notorious
2019/20 bushfire season. Figure <ref> displays the
sensitivity of the SPEI series at the Texas post office to different
time scales and fitted distributions. Larger time scales produce a
smoother index across time, however, all time scales indicate an extreme
drought (corresponding to -2 in SPEI) in 2020, confirming the severity
of the drought across different time horizons. Moreover, the chosen
distribution has less influence on the index, with general extreme value
distribution tending to produce more extreme outcomes than log-logistic
distribution for the extreme events (index  2 or

Spatial distribution of Standardized
Precipitation Index (SPI-12) in Queensland, Australia during two major
flood and drought events: 2010/11 and 2019/20. The map shows a
continuous wet period during the 2010/11 flood period and a mitigated
drought situation, after its worst in 2019 December and 2020 January,
likely due to the increased rainfall in February from the meteorological

Time series plot of Standardized
Precipitation-Evapotranspiration Index (SPEI) at the Texas post office
station (highlighted by a diamond shape in panel a). The SPEI is
calculated at four time scales (6, 12, 24, and 36 months) and fitted
with two distributions (Log Logistic and GEV). The dashed line at -2
represents the class “extreme drought” by the SPEI. A larger time
scale gives a smoother index series, while also taking longer to recover
from an extreme situation as seen in the 2019/20 drought period. The
SPEI values from the two distributional fits mostly agree, while GEV can
result in more extreme values, i.e. in 1998 and 2020.

§.§ Does a puff of change in variable weights cause a tornado in

The Global Gender Gap Index (GGGI), published annually by the World
Economic Forum, measures gender parity by assessing relative gaps
between men and women in four key areas: Economic Participation and
Opportunity, Educational Attainment, Health and Survival, and Political
Empowerment (World Economic Forum 2023). The index, defined on 14
variables measuring female-to-male ratios, first aggregates these
variables into four dimensions (using the linear combination given by
in Table <ref>). The weights are the
inverse of the standard deviation of each variable, scaled to sum to 1,
thus ensuring equal relative contribution of each variable to each of
the four new variables. These new variables are then combined through
another linear combination ( in
Table <ref>) to form the final index value.
Figure <ref> illustrates that the pipeline is constructed by
applying the dimension reduction module twice on the data. The data for
GGGI does not needs to be transformed or scaled so these steps are not
included, but they might still need to be used for other similar

Index pipeline for the Global Gender Gap
Index (GGGI). The index is constructed as applying the module dimension
reduction twice on the data.

Weights for the two applications of dimension reduction to compute the
Global Gender Gap Index. V-wgt is used to compute four new variables
from the original 14. These are then equally combined to get the final
index value.

Variable     V-wgt     Dimension     D-wgt     wgt

Labour force participation     0.199     Economy     0.25     0.050

Wage equality for similar work     0.310               0.078

Estimated earned income     0.221               0.055

Legislators senior officials and managers     0.149               0.037

Professional and technical workers     0.121               0.030

Literacy rate     0.191     Education     0.25     0.048

Enrolment in primary education     0.459               0.115

Enrolment in secondary education     0.230               0.058

Enrolment in tertiary education     0.121               0.030

Sex ratio at birth     0.693     Health     0.25     0.173

Healthy life expectancy     0.307               0.077

Women in parliament     0.310     Politics     0.25     0.078

Women in ministerial positions     0.247               0.062

Years with female head of state     0.443               0.111

The 2023 GGGI data is available from the Global Gender Gap Report 2023
in the country's economy profile and can be accessed in the
package as  with
Table <ref> as . The index can be
reproduced with:

gggi %%
init(id = country) %%
add_meta(gggi_weights, var_col = variable) %%
index_new = aggregate_linear(
weight = weight))

After initializing the  object and attaching the
as meta-data, a single linear combination within
the dimension reduction module is applied to the 14 variables (from
column  to
), using the weight
specified in the  column of the attached metadata. While
computing the index from the original 14 variables, it remains unclear
how the missing values are handled within the index, which impacts 68
out of the total 146 countries. However, after aggregating variables
into the four dimensions, where no missing values exist, the index is
reproducible for all the countries.

Figure <ref> illustrates doing sensitivity analysis for
GGGI, for a subset of 16 countries. Frame 12 shows the dot plot of the
original index values sorted from highest to lowest. Leading the
rankings are the Nordic countries (Iceland, Norway, and Finland) and New
Zealand. The index values are between 0 and 1, and indicate proportional
difference between men and women, with a value of 0.8 indicating women
are 80% of the way to equality of these measures. There is a gap in
values between these countries and the middle group (Brazil, Panama,
Poland, Bangladesh, Kazakhstan, Armenia, and Slovakia), and another big
drop to the next group (Pakistan, Iran, Algeria, and Chad). Afghanistan
lags much further behind.

Exploring the sensitivity of the GGGI, by
varying the politics component's contribution, for a subset of
countries. Each panel shows a dotplot of the index values, computed for
the linear combination represented by the segment plots below. Frame 12
shows the actual GGGI values, and countries are sorted from highest to
lowest on this. Frames 1 and 6 show the GGGI if the politics component
is reduced. Frames 18, 24, 29 show the GGGI when the politics component
is increased. The most notable feature is that Bangladesh's GGGI drops
substantially when politics is removed, indicating that this component
plays a large role in its relatively high value. Also, politics plays a
substantial role in the GGGI's for the top ranked countries, because
each of them drops, to the state of being similar to the middle ranked
countries when the politics component's contribution is reduced. The
animation can be viewed at https://vimeo.com/847874016.

To make a simple illustration of sensitivity analysis, we slightly vary
the weight for politics, between 0.07 and 0.52, while maintain equal
weights among other dimensions. This can be viewed as an animation to
examine change in relative index values as a response to the changing
weights. This visualization technique, which presents a sequence of data
projections, is referred to as a “tour” and the specific kind of tour
used here to move between nearby projections is known as a “radial
tour” (see Andreas Buja et al. (2005), Wickham et al. (2011) and
Spyrison and Cook (2020) for more details).

Frames 1 and 6 show linear combinations where politics contributes less
than the original. It is interesting to note that the gap between the
Nordic countries and the middle countries dissipates, indicating that
this component was one reason for the relatively higher GGGI values of
these countries. Also interesting is the large drop in value for

Frames 18, 24, 29 show linear combinations where politics contributes
more than the original. The most notable feature is that Bangladesh
retains it's high index value whereas the other middle group countries
decline, indicating that the politics score is a major component for
Bangladesh's index value.

Ideally, an index should be robust against minor changes in its
construction components. This is not the case with GGGI, where small
changes to one component lead to fairly large change in the index. The
modular pipeline framework for computing the index makes it easy to
conduct this type of sensitivity analysis, where one or more components
are perturbed and the index re-calculated. One aspect of the GGGI not
well-described in the Global Gender Gap Report is handling of missing
values that are present in the initial variables for many countries,
something that is common for this type of data. This could also be made
more transparent with the dimension reduction module, by specifying an
imputation method or providing warnings about missing values.

§.§ Decoding uncertainty through the wisdom of the

Errors in measurement, variability and sampling error, may arise at
various stages of the pipeline calculation, including from different
parameterization choices, as illustrated from
Section <ref>, or from the statistical summarization
procedures applied in the pipeline. Although it may not be possible to
perfectly measure these errors, it is important that they are recognised
and estimated for an index, so that it is possible to compute confidence
intervals. In this example, the Texas post office station highlighted in
Figure <ref> is used to illustrate one possibility
to compute a confidence interval for the Standardized Precipitation
Index (SPI). Bootstrapping is used to account for the sampling
uncertainty in the distribution fit step of the index pipeline and to
assess its impact on the SPI series.

In SPI, the distribution fit step fits the gamma distribution to the
aggregated precipitation series separately for each month. This results
in 32 or 33 points, from January 1990 to April 2022, for estimating each
set of distribution parameters. To account for this sampling uncertainty
with these samples, bootstrapping is used to generate replicates of the
aggregated series. In the  package, this bootstrap
sampling is activated when the argument  is set to a
value other than the default of 1. In the following code, the
Standardized Precipitation Index (SPI) is calculated using a time scale
of 24. The bootstrap procedure samples the aggregated precipitation
() for 100 iterations () and then
fits the gamma distribution. The resulting gamma probabilities are then
transformed into normal densities in the normalizing step with

DATA |
temporal_aggregate(.agg = temporal_rolling_window(prcp, scale = 24)) |
distribution_fit(.fit = dist_gamma(var = ".agg", method = "lmoms",
.n_boot = 100)) |
normalise(.index = norm_quantile(.fit))

The confidence interval can then be calculated using the quantile method
from the bootstrap samples. Figure <ref> presents the
80% and 95% confidence interval for the Texas post office station, in
Queensland, Australia. From the start of 2019 to 2020, the majority of
the confidence intervals lie below the extreme drought line (SPI = -2),
suggesting a high level of certainty that the Texas post office is
suffering from a drastic drought. The relatively wide confidence
interval, as well as during the excessive precipitation events in
1996-1998 and 1999-2000, suggests a high variation of the gamma
parameters estimated from the bootstrap samples and its difficulty to
accurately quantify the drought and flood severity in extreme events.

80% and 95% confidence interval of
the Standardized Precipitation Index (SPI-24) for the Texas post office
station, in Queensland, Australia. A bootstrap sample of 100 is taken
from the aggregated precipitation series to estimate gamma parameters
and to calculate the index. The dashed line at SPI = -2 represents an
extreme drought as defined by the SPI. Most parts of the confidence
intervals from 2019 to 2020 sit below the extreme drought line and are
relatively wide compared to other time periods. This suggests that while
it is certain that the Texas post office is suffering from a drastic
drought, there is considerable uncertainty in quantifying its severity,
given the extremity of the event.

§ CONCLUSION

The paper introduces a data pipeline comprising nine modules designed
for the construction and analysis of indexes within the tidy framework.
The pipeline offers a modular workflow to compute indexes with different
parameterizations, to test minor changes to the original index
definition, and to quantify uncertainties. The framework proposed in the
paper is universal to index across diverse domains. Examples are
provided, including the drought indexes (SPI and SPEI) and the Global
Gender Gap Index (GGGI), to demonstrate the index calculation with
different time scales and distributions, to illustrate how slight
adjustment of linear combination weights impact the index, and to
calculate confidence intervals on the index.

§ ACKNOWLEDGEMENT

This work is funded by a Commonwealth Scientific and Industrial Research
Organisation (CSIRO) Data61 Scholarship. Nicolas Langrené acknowledges
the partial support of the Guangdong Provincial Key Laboratory IRADS
(2022B1212010006, R0400001-22) and the UIC Start-up Research Fund
UICR0700041-22. The article is created using Quarto (Allaire et al.
2022) in R (R Core Team 2021). The source code for reproducing this
paper can be found at:

§ REFERENCE

Alahacoon, Niranga, and Mahesh Edirisinghe. 2022. “A Comprehensive
Assessment of Remote Sensing and Traditional Based Drought Monitoring
Indices at Global and Regional Scale.” Geomatics, Natural
Hazards and Risk 13 (December): 762–99.

Allaire, J. J., Charles Teague, Carlos Scheidegger, Yihui Xie, and
Christophe Dervieux. 2022. Quarto (version 1.2).

Becker, William, Giulio Caperna, Maria Del Sorbo, Hedvig Norlen, Eleni
Papadimitriou, and Michaela Saisana. 2022. “COINr: An R Package for
Developing Composite Indicators.” Journal of Open Source
Software 7 (78): 4567. <https://doi.org/10.21105/joss.04567>.

Beguería, Santiago, and Sergio M. Vicente-Serrano. 2017. SPEI:
Calculation of the Standardised Precipitation-Evapotranspiration Index.

Buja, A, D Asimov, C Hurley, and JA McDonald. 1988. “Elements of a
Viewing Pipeline for Data Analysis.” In Dynamic Graphics for
Statistics, 277–308. Wadsworth, Belmont.

Buja, Andreas, Dianne Cook, Daniel Asimov, and Catherine Hurley. 2005.
“Computational Methods for High-Dimensional Rotations in Data
Visualization.” Handbook of Statistics 24: 391–413.

Fisher, Ronald A. 1936. “The Use of Multiple Measurements in Taxonomic
Problems.” Annals of Eugenics 7 (2): 179–88.

Fisher, Ronald Aylmer. 1970. “Statistical Methods for Research
Workers.” In Breakthroughs in Statistics: Methodology and
Distribution, 66–70. Springer.

Grenié, Matthias, and Hugo Gruson. 2023. fundiversity: Easy
Computation of Functional Diversity Indices.

Hao, Zengchao, and Vijay P. Singh. 2015. “Drought Characterization
from a Multivariate Perspective: A Review.” Journal of
Hydrology 527 (August): 668–78.

Hotelling, Harold. 1933. “Analysis of a Complex of Statistical
Variables into Principal Components.” Journal of Educational
Psychology 24 (6): 417.

Jones, Brenda, and Jean Andrey. 2007. “Vulnerability Index
Construction: Methodological Choices and Their Influence on Identifying
Vulnerable Neighbourhoods.” International Journal of Emergency
Management 4 (2): 269–95.

Kuhn, Max, and Hadley Wickham. 2020. Tidymodels: A Collection of
Packages for Modeling and Machine Learning Using Tidyverse Principles.

Laimighofer, Johannes, and Gregor Laaha. 2022. “How Standard Are
Standardized Drought Indices? Uncertainty Components for the SPI &
SPEI Case.” Journal of Hydrology 613 (October): 128385.

Martin, Steve. 2023. Gpindex: Generalized Price and Quantity
Indexes. <https://CRAN.R-project.org/package=gpindex>.

McKee, Thomas B, Nolan J Doesken, John Kleist, et al. 1993. “The
Relationship of Drought Frequency and Duration to Time Scales.” In
Proceedings of the 8th Conference on Applied Climatology,
17:179–83. 22. Boston, MA, USA.

OECD, European Union, and Joint Research Centre - European Commission.
2008. Handbook on Constructing Composite Indicators:
Methodology and User Guide. OECD.

R Core Team. 2021. R: A Language and Environment for Statistical
Computing. Vienna, Austria: R Foundation for Statistical Computing.

Saisana, M., A. Saltelli, and S. Tarantola. 2005. “Uncertainty and
Sensitivity Analysis Techniques as Tools for the Quality Assessment of
Composite Indicators.” Journal of the Royal Statistical Society
Series A: Statistics in Society 168 (2): 307–23.

Spyrison, Nicholas, and Dianne Cook. 2020. “Spinifex: An R Package
for Creating a Manual Tour of Low-Dimensional Projections of
Multivariate Data.” The R Journal 12: 243–57.

Sutherland, Peter, Anthony Rossini, Thomas Lumley, Nicholas Lewin-Koh,
Julie Dickerson, Zach Cox, and Dianne Cook. 2000. “Orca: A
Visualization Toolkit for High-Dimensional Data.”
Journal of Computational and Graphical Statistics 9 (3): 509–29.

Svoboda, Mark, Brian Fuchs, et al. 2016. “Handbook of Drought
Indicators and Indices.” Drought and Water Crises: Integrating
Science, Management, and Policy, 155–208.

Tate, Eric. 2012. “Social Vulnerability Indices: A Comparative
Assessment Using Uncertainty and Sensitivity Analysis.” Natural
Hazards 63 (2): 325–47.

———. 2013. “Uncertainty Analysis for a Social
Vulnerability Index.” Annals of the Association of American
Geographers 103 (3): 526–43.

Vicente-Serrano, Sergio M., Santiago Beguería, and Juan I. López-Moreno.
2010. “A Multiscalar Drought Index Sensitive to Global
Warming: The Standardized Precipitation Evapotranspiration
Index.” Journal of Climate 23 (7): 1696–1718.

Wickham, Hadley. 2014. “Tidy Data.” Journal of Statistical
Software 59 (September): 1–23.

Wickham, Hadley, Mara Averick, Jennifer Bryan, Winston Chang, Lucy
D'Agostino McGowan, Romain François, Garrett Grolemund, et al. 2019.
“Welcome to the Tidyverse.” Journal of Open Source Software
4 (43): 1686. <https://doi.org/10.21105/joss.01686>.

Wickham, Hadley, Dianne Cook, Heike Hofmann, and Andreas Buja. 2011.
“Tourr: An R Package for Exploring Multivariate Data with
Projections.” Journal of Statistical Software 40 (2).

Wickham, Hadley, Michael Lawrence, Dianne Cook, Andreas Buja, Heike
Hofmann, and Deborah F. Swayne. 2009. “The Plumbing of Interactive
Graphics.” Computational Statistics 24 (2): 207–15.

World Economic Forum. 2023. “The Global Gender Gap Report 2023.”

Xie, Yihui, Heike Hofmann, and Xiaoyue Cheng. 2014. “Reactive
Programming for Interactive Graphics.” Statistical
Science 29 (2): 201–13.

Zargar, Amin, Rehan Sadiq, Bahman Naser, and Faisal I Khan. 2011. “A
Review of Drought Indices.” Environmental Reviews 19 (NA):
333–49. <https://www.jstor.org/stable/envirevi.19.333>.


# In Situ 3D Spatiotemporal Measurement of Soluble Biomarkers in Organoid
Culture

Alexander McGhee  Eric McGhee  Jack E Famiglietti  W Gregory Sawyer
Biomolecular Science & Engineering Division Naval Research Laboratory
Washington, DC 20375  Research and Development Aurita Bioscience Gainesville,
FL 32601  Department of Mechanical Engineering, University of Wisconsin-
Madison, Madison, WI 53706, USA  Department of Mechanical & Aerospace
Engineering Herbert Wertheim College of Engineering University of Florida
Gainesville, FL 32601

###### Abstract

Advanced cell culture techniques such as 3D bio-printing and hydrogel-based
cell embedding techniques harbor many new and exciting opportunities to study
cells in environments that closely recapitulate in-vivo conditions.
Researchers often study these environments using fluorescence microscopy to
visualize the protein association with objects such as cells within the 3D
environment, yet quantification of concentration profiles in the
microenvironment has remained elusive. Here, we present a method to
continuously measure the time-dependent concentration gradient of various
biomarkers within a 3D cell culture assay using bead-based immunoassays to
sequester and concentrate the fluorescence intensity of these tagged proteins.
This assay allows for near real-time in situ biomarker detection and enables
spatiotemporal quantification of biomarker concentration. Snapshots of
concentration profiles can be taken, or time series analysis can be performed
enabling time-varying biomarker production estimation. Example assays utilize
an osteosacroma tumoroid as a case study for a quantitative single-plexed gel
encapsulated assay, and a qualitative multi-plexed 3D bioprinted assay. In
both cases, a time-varying cytokine concentration gradient is measured. An
estimation for the production rate of the IL-8 cytokine per second per
osteosarcoma cell results from fitting an analytical function for continuous
point source diffusion to the measured concentration gradient and reveals that
each cell produces approximately two IL-8 cytokines per second. Proper
calibration and use of this assay is exhaustively explored for the case of
diffusion-limited Langmuir kinetics of a spherical adsorber.

###### keywords:

Biomarker; Bead-based ELISA; 3D cell culture; Immunoassay; Osteosarcoma

††journal: arxiv

##  1 Introduction

The tumor microenvironment creates a protective immunosuppressive barrier and
represents a significant obstacle to therapies [1, 2, 3, 4]. Local
measurements of the secreted cytokines and signaling molecules involved in
immunoregulation within this interphase between tumor tissue and the non-
cancerous surrounding tissues is particularly challenging due to difficulties
in probing this exquisitely small volume. On cellular dimensions ( O(10)
$\upmu$m in thickness) this total volume of the interphase region can be as
small as a microliter (depending on the tumor size); in vitro interphase
volumes further exacerbate the difficulty of measurement since the tumor
explants and tumoroids are often much smaller than in vivo with interphase
volumes of O(1-10)nL. The second major challenge with local measurements is
that the in vitro tumor models are increasingly becoming three-dimensional,
which introduces significant sampling and imaging challenges [5, 6]. Accurate
measurements of local concentrations over nano-liter volumes in 3D pose
significant challenges in both collection and detection. To overcome these
challenges we have suspended ELISA beads within traditional matrigel-based 3D
culture as well as 3D bioprinted Liquid Like Solid culture. These suspended
ELISA beads are then able to participate in capture dynamics and when imaged
in a standard confocal microscope, allow for quantifiable in-situ detection of
cytokines in the immediate tumor vicinity. In this manuscript, we will discuss
the approach of using fields of ELISA beads to measure biometric data such as
spatiotemporal concentration gradients and cellular production rates.

Standard benchtop in vitro experiments enjoy a plethora of techniques to
quantify biomarker concentrations[7]. Often centered around media samples
captured during cell feeding, the majority of immunoassays offer robust
measurement of biomarkers for time-point analysis. Although these benchtop
immunoassays are the gold standard for many biomarkers, they typically do not
have a sensitivity high enough for biomarkers that are scarcely produced and
can only measure the average biomarker concentration in the entire assay
volume. To overcome this issue, researchers have developed an in-situ
immunoassay called FluoroSpot which enables the simultaneous measurement of
multiple biomarkers within 2D cell culture at high sensitivity [8]. In this
assay, the cells act as a point source for the diffusion of biomarkers, thus,
a gradient of captured biomarkers surrounds each cell; this gradient is
detected through the use of a fluorophore-conjugated antibody which binds to
the captured biomarkers in a fashion commonly referred to as the sandwich
ELISA [9]. A fluorescence measurement device can then be used to measure the
resulting intensity gradient and further analyzed to reveal the concentration
gradient profile. Although the FluoroSpot method works well for 2D cell
culture, there is currently no available analog for 3D culture.

Recent advancements in immunoassay development have enabled the detection of
multiple biomarkers from a single sample by utilizing antibody-conjugated
microbeads, commonly known as bead-based multiplexed ELISA assays. These
assays are becoming increasingly popular with many off-the-shelf biomarker
detection kits available from numerous manufacturers. Most of these assays
utilize the sandwich ELISA method to capture and detect biomarkers in a
solution through the use of proprietary flow cytometers such as Magpix®,
LUNARIS™, and Quanterix® [10]. Manufacturers of bead-based assays typically
require the use of custom-designed flow cytometers for the detection of
fluorescence intensity with their bead-based immunoassay to maximize the
fluorescence intensity signal and achieve high sensitivity and linearity of
the assay. Although these detection systems are custom-built, they all utilize
the same basic principles outlined in Fulton et al[11]. Briefly, microbeads
with one or more unique fluorescence signatures are used to identify the
biomarker specificity of each bead while other fluorescence signatures are
used to detect biomarker concentration through a calibration curve transfer
function[12]. Great care is taken in the design of the bead adsorption
kinetics to optimize and further enhance the detectability of the
biomarkers[13].

Here, we utilize these bead-based immunoassays to develop an in-situ biomarker
detection method for 3D cell culture. This technique uses fluorescence
microscopy to detect fluorescence intensities corresponding to local biomarker
concentration, as overviewed in (Fig. 1(B)). Using this method, we are able to
quantify the biomarker concentration in a singleplex assay and show that this
technique is able to capture in situ spatiotemporal information of biomarker
concentration in near real-time. Qualitative detection of biomarker production
is also described in a multiplex assay using a much-simplified version of the
quantitative method. Adaptation of the quantitative method for use in such a
multiplex assay is straightforward but not presented. Furthermore, protocols
for incorporating capture beads into standard 3D culture systems are developed
for liquid-like solid 3D bioprinting [14] and typical gel-based 3D culture
[15]. These results show the first successful implementation of a commercially
available real-time in-situ immunoassay for 3D cell culture and 3D bioprinting
assays in both quantitative singleplex, and qualitative multiplex fashions.

Figure 1: Overview of traditional ELISA methods versus our approach A) typical
bench-top methods for measuring biomarkers with immunoassays from (a) 2D cell
culture and (b) 3D cell culture rely on the diffusion of biomarkers into a
bath of growth media which is collected at various time points. These methods
dilute the biomarker concentration by the total volume of growth media in the
bath. Furthermore, the sample collection and testing can often take multiple
days to obtain results. B) our method pictured in (c) utilizes immunosorbent
beads scattered throughout the volume of a 3D cell culture assay. The bead
fluorescence intensity corresponds to the local concentration of biomarker
within a local bead volume of approximately 20 nL and these local
concentrations may decrease with increasing distance from a source following a
concentration gradient of a continuous point source (d). For a continuous
source of biomarker from a collection of cells (e), each biomarker (f)
released from a cell will bind to a detection antibody (g) to form a pair (h),
this antibody biomarker pair then adsorbs to the immunosorbant bead (i) via
conjugated capture antibodies (j). The antibody biomarker pair will desorb (k)
from the capture antibody at some reverse reaction rate $k_{r}$. Through this
process, the capture beads will produce a time-varying signal proportional to
the local concentration of biomarker at any time.

##  2 Assay Creation

###  2.1 Building a Qualitative Biomarker Detection Assay

The simplest form of this assay can be implemented in a qualitative manner
with relative ease resulting in spatial information of biomarker presence
within detectable ranges. In practice, a combination of cytometric beads and
detection antibody within a 3D assay can provide a boolean interpretation of
biomarker presence. In short, if a mean fluorescence intensity ($MFI$) exists
in the proper channel after background subtraction, then the biomarker is
present at detectable levels, and vice versa. However, the authors suggest
careful consideration of Appendix B before use of a qualitative assay due to
various issues related to the kinetics of singlepot-ELISA. For example, the
use of custom cytometric beads or overly dense bead concentrations may lead to
poor or misleading results. Nevertheless, this level of the assay can allow a
user to rapidly probe an experiment to provide real-time spatial information
about the state of the organoid in a 3D culture system.

###  2.2 Building a Quantitative Biomarker Detection Assay

In this manuscript, we utilize cytometric beads in a manner that is far from
the manufacturer’s intended use. For example, measurement of the fluorescence
intensity is achieved using a fluorescence microscope which has many
associated issues with capturing quantitative information[16, 17, 18]. These
issues ultimately lead to relatively large uncertainties of measured biomarker
concentration, however, the ability to capture the spatiotemporal
concentration of biomarker in real-time remains an attractive feature. The
adaptation of these commercial assays to capture in-situ quantitative
information for 3D culture systems demands attention to two critical
considerations: the conversion of $MFI$ to biomarker concentration and
determination of acceptable bead density. Typical bead-based assays assume a
pure convection boundary at the bead surface coupled with an abundance of
biomarker much beyond the total capture capacity of the beads in solution.
These assumptions cannot be made upon adaptation to 3D culture systems which
typically have minimal or no convection[19]. This leads to diffusion-limited
reactions in the 3D environment and thus special considerations must be taken
to ensure that measured bead $MFI$ is calibrated correctly. For this
calibration, a standard curve is generated under pure diffusion, and we
introduce a dynamic measurement based on the concentration gradient of the
biomarker with known boundary conditions. The dynamic measurement is based on
a flat wall diffusion gradient of biomarker from a fully mixed convective
chamber into a 3D culture analog. The purpose of the diffusion chamber is to
ensure the density of capture beads used in the assay is appropriate to the
assay kinetics. Since each capture bead sequesters some amount of biomarker
from the local environment, capture beads may compete with one another if they
are too densely packed, effectively lowering the detected local concentration
as shown in Fig. 2(D). Thus, the capture bead concentration must be
experimentally verified to ensure an appropriate density exists.

##  3 Methodology in Demonstration

###  3.1 Culture of osteosarcoma spheroids

OS521 was cultured in a T-75 flask (Corning) using a complete culture medium
with 0.3 mg/mL G418 (Mediatech) at $37^{\circ}$C, 5% $CO_{2}$ [20]. After
reaching confluence, cells were released from the flask using a trypsin
dissociation reagent (TrypLE, Gibco) and resuspended in culture media to a
density of approximately 106 cells/mL. This cell suspension was then
distributed to an ultra-low attachment round-bottom 96-well plate (Corning )
and incubated at $37^{\circ}$C, 5% $CO_{2}$ for 72 hours with daily media
changes until tight spheroids were formed. Immediately before use in assays,
spheroids were dyed using CellTracker Green CMFDA dye using standard
manufacturer-specified protocols prior to the experiment.

###  3.2 preparation of gel-based distributed bead assay

A human IL-8 cytometric bead array (CBA) kit (BD biosciences) consisting of
three components, IL8 capture beads, fluorophore-conjugated detection
antibody, and lyophilized human IL8 standard was used in a matrigel based
culture of osteosarcoma spheroids. First, the detection antibody was diluted
in a solution of culture media as well as a solution of matrigel at a ratio of
1:200. This ratio may vary between assays and manufacturers since the
detection antibody concentration is not consistent. For this assay, the target
detection antibody concentration allows for 10x the number concentration of
detection antibodies to the highest local number concentration of biomarkers
expected. The detection antibody concentration can be determined using the
protocol outlined in Appendix C.4. Next, the capture beads were added to the
matrigel solution at a concentration of 50 beads/$\upmu$L (the concentration
of capture beads given in BD CBA kits is approximately 3,000 beads/$\upmu$L),
gently pipette mixed, and placed on ice. An osteosarcoma spheroid, chilled on
ice, was then added to the matrigel-bead solution and gently pipette mixed.
The solution was then pipetted onto a glass-bottom well to form a dome, the
spheroid was gently moved towards the center, and the sample was placed in an
incubation chamber at $37^{\circ}$C, 5% $CO_{2}$ until the matrigel was set
(approximately 2 minutes). Finally, the media + detection antibody solution
was added to the culture dish ensuring the matrigel was fully submerged. The
sample was then placed on a fluorescent microscope and imaged every 4 hours
using laser channels of: GFP, to capture the OS521 spheroid; cy5, to capture
bead positions; and TRITC, to capture detection antibody intensity. The
resulting flourescent images can be seen in Fig. 3(A).

###  3.3 preparation of 3D printed bead assay

Four BD CBA kits of human IL-11, human IL-8, human IL-6, and human VEGF are
checked for cross-reactivity by consulting the BD CBA manual. A solution of
LLS is created using culture media according to the protocol outlined in
Bhattacharjee et al.[21]. Detection antibodies from each kit are diluted into
culture media as well as a solution of LLS at a ratio of 1:200 v/v. Capture
beads from each kit are combined with the LLS solution in separate aliquots at
a concentration of 50 beads/$\upmu$L. In a separate aliquot, an osteosarcoma
spheroid is combined with a small amount of the LLS solution and printed into
a glass-bottom petri dish filled with LLS following a method outlined in
Bhattacharjee et al.[21]. Each capture bead solution is then printed in a
4-arm configuration centered about the osteosarcoma spheroid as shown in Fig.
3(C). Finally, the media + detection antibody solution is carefully added atop
the printed structures being careful not to disturb the sample. The sample is
then placed on a fluorescent microscope and imaged every 4 hours using both
the bead detection (cy5) channel as well as the detection antibody (TRITC)
channels.

##  4 Results

###  4.1 Calibration

We developed a PDMS device for the calibration of $MFI$ versus the biomarker
concentration using a serial dilution of biomarker standard as well as the
validation of bead kinetics. For this study, we used the human IL-8 CBA kit to
produce a calibration curve and used the results to produce a concentration
gradient profile for two bead concentrations - 50 and 500 beads/$\upmu$L. The
results presented in Fig. 2(D) show that the bead concentration of 50
beads/$\upmu$L produces a concentration profile similar to the analytical
model described by eq.2 and the bead concentration of 500 beads/$\upmu$L
greatly underpredicts the analytical model. Kinematic variables of the assay
were determined using the protocols described in Appendix C finding the
antibody + biomarker diffusion coefficient $D$ to be 55 $\upmu$ $m^{2}$/s, the
reverse reaction rate $k_{r}$ to be approximately 0.001 $1/s$, and the number
of capture sites per bead N to be 2.6 E6 sites. Using these kinematic
variables with a biomarker concentration of the top standard $C$ = 2500 pg/mL,
the equilibrium bead coverage $\theta_{eq}$ was found to be 0.12. Using these
variables in eq.4 reveals the sequestration ratio $\Gamma$ for bead
concentrations of 50 beads/$\upmu$L and 500 beads/$\upmu$L to be 0.1 and 1
respectively. This result demonstrates that eq.4 can accurately predict the
appropriate bead density when the sequestration ratio is $\Gamma<<1$.

Figure 2: Calibration device for standard curve and dynamic diffusion curve
analysis. A) a cross-sectional view of the calibration device showing the
standard curve wells, reservoir, diffusion chamber, and membrane filter. The
colorscale from light pink to dark orange signifies the concentration of
biomarkers in from low to high respectively. B) an isometric view of the
calibration device shows the overall construction. C) The standard curve is
generated using a linear spaced concentration of biomarker from 500 to 2500
pg/mL which corresponds to this assays linear range. The $MFI$ of each bead is
measured and then normalized ($nMFI$) by the top of the dynamic range. D)
Beads in the diffusion chamber are measured after 10 hours and the $nMFI$ is
transferred into cytokine concentration using the standard curve. An
analytical model for flat wall diffusion (black dashed line) is compared to
two bead densities $\Gamma=0.1$ (low density, red circles) and $\Gamma=1$
(high density, blue circles). These results show that as $\Gamma$ becomes much
less than one, the beads measurement more closely follows the expected value.
Figure 3: Quantitative single-plex, and qualitative multi-plexed assay results
for an osteosarcoma spheroid. A) a max intensity z-projection shows a 500
$\upmu$m diameter osteosarcoma spheroid in green on the left with a white
dashed line corresponding to the area where no detection beads are found.
Breakouts of six beads at various radial distances from the spheroid show a
decreasing intensity with increasing radial distance. B) The resulting
transformation of the measured bead intensities to IL8 concentration for
timepoints 100, 300, and 500 minutes shows good agreement with the analytical
model for diffusion from a continuous point source with a production rate $q$
of approximately 2 IL-8 proteins per cell per second. C) A 3D printed bead
array surrounding an osteosarcoma tumoroid shows the relative position of four
bead types IL6, IL8, IL11, and VEGF. D) The detectability of each biomarker
can be compared directly with IL8 and VEGF showing relatively high readings
for each, and relatively low readings for IL-11 and IL-6. These readings would
indicate, using the qualitative assay, that IL-8 and VEGF are present while
IL-11 and IL-6 are not being produced in sufficient quantities to be measured.

###  4.2 Singleplex matrigel-based assay

OS521 is known to produce the IL-8 cytokine in abundance, so we chose to focus
on this cytokine for our trial study[22]. After the IL-8 capture bead density
and $MFI$ vs concentration calibration curve were found, we analyzed the
resulting intensity gradient of the beads surrounding the spheroid to recover
the concentration gradient for three time points 100, 300, and 500 minutes as
shown in Fig. 3(A,B). An analytical solution for a continuous point source
given by Crank[23] in eq.1

$C\left(r,t\right)=\frac{q}{4\pi Dr}erfc\left[\frac{r}{2\sqrt{Dt}}\right]$ (1)

is fitted to the experimental data using the rate of production $q$ as a free
variable with $r$ the radial distance from the point source as the independent
variable. The diffusion coefficient for the detection antibody $+$ biomarker
$D$, and the assay time in seconds $t$ are known constants for each curve
analyzed. The rate of cytokine production $q$ for the 500$\upmu$m diameter
spheroid (estimated to be $\sim$10k cells) was found to be nearly 2 IL-8
proteins per cell per second by fitting eq.1 to the concentration gradient at
multiple timepoints.

###  4.3 Multiplex LLS-based assay

Bead-based immunoassays are often used to measure multiple biomarkers from the
same sample, and thus the detection and capture antibodies have been optimized
to ensure little to no cross-reactivity. We utilize this capability in a 3D
bioprinted assay by separating the bead types by location in a 4-arm
configuration centered around the spheroid as shown in Fig. 3(C). The bead
intensity gradient for the IL-11, IL-8, IL-6, and VEGF was measured for
multiple time points and presented in Fig. 3(D). These results show that the
four cytokines measured have different concentration gradient time profiles
and thus four different production rates. No calibration was completed for
this multiplexed assay so results remain qualitative.

##  5 Discussion

The rise of 3D culture provides exciting opportunities to study in vitro
biological systems in their closest analogous configuration to in vivo[24, 25,
26, 19]. These methods offer potentially enhanced biological complexity
through mass transport modification, 3D intercellular contact, cell-matrix
signaling, customizable geometric configuration of tissues, tissue shielding,
and allow complex immune interaction. The addition of perfusion into these 3D
culture systems allows for long-duration cultures providing opportunity for
emergent biological complexity to flourish in vitro[27, 28, 29, 30]. Such
methods are of keen interest in the pharmaceutical industry as the development
and deployment of in vivo-like high throughput drug screening may offer better
drug development outcomes[31, 32]. Unfortunately, interrogating biochemical
compositions of 3D in vitro systems using standard bioassays can be extremely
challenging to near impossible. Therefore, new developments in bioassay
technologies are needed to provide the insights researchers have come to know
as standard for 2D in vitro systems.

In this study, we have examined the potential for the use of a commercially
available bead-based immunosorbent assay as a viable method for the in-situ
measurement of cytokine concentration within 3D culture assays. This assay
method not only represents the first of its kind measurement for 3D culture
systems but is also easily implemented using this manuscript as a guide.
Additionally, we have shown that this assay can be implemented in simple gel-
based culture systems as well as more advanced 3D bioprinting systems at low
cost. Typical bead-based assays supply approximately 300,000 beads, and more
than enough detection antibody to coat each bead at the highest standard. The
3D assays used in this study used a bead density of 50 beads/$\upmu$L with a
total gel volume of approximately 100$\upmu$L per sample. At this ratio, this
experiment can be repeated 60 times at a total cost of a single bead-based
assay kit. At the time of this writing, an IL-8 CBA kit costs approximately
$278.00 making the cost of each experiment less than five dollars. If a 3D
printing technique were implemented, as we have done in this manuscript, this
cost can be reduced significantly.

Bead-based immunoassays are not typically designed to be used in the manner we
prescribe here. Some assays may not be viable due to incompatible assay
kinetics, background saturation, cross-reactivity, etc. For example, the
immunoassay must be capable of single-pot capture and detection which requires
the antigen to have unique binding epitopes for the paratopes of the capture
and detection antibodies. Many assays are capable of single-pot capture and
detection but do not recommend this incubation method due to a reduction in
assay sensitivity via the inclusion of detection antibodies within the sample
which may obscure signals at the low end of the detection range. In our
experience, these factors typically result in a 10 to 20% reduction of the
assays reported range when used in our method. However, since our assay
assumes the pure diffusion of biomarkers from one or more point sources, the
local concentration of biomarkers accumulates over time and may eventually
become detectable. Finally, the active implementation of this assay may limit
biological interaction. The inclusion of detection antibody and subsequent
conjugation of the detection antibody in solution may affect intercellular
chemokine signaling and as a result, change cellular behavior. This challenge
could be circumvented by cycling of detection antibody for periods before
imaging, especially in thin sections of 3D gel-based assays or through
perfusion pulsation in LLS-based experimental platforms.

##  6 Conclusions

We developed a novel method for the utilization of commercially available
bead-based immunosorbent assays for use in 3D cell culture platforms. Two 3D
culture systems, a simple matrigel-based assay, and an advanced 3D bioprinting
assay were studied using osteosarcoma spheroids to generate and measure
multiple cytokines of interest. Protocols to back out the kinematic variables
from commercial cytometric assays have been developed and allow users to
easily calculate the ideal beads density and the detection antibody
concentration appropriate for an individualized assay. Finally, we show that
this assay can be easily implemented and adapted for use in multiple 3D
culture systems in both qualitative and quantitative capacities.

##  Appendix A Bead Measurements

Commercial bead-based multiplex immunoassays use one or more fluorescent dyes
incorporated into the bead to identify the antigen specificity within a
mixture of many bead types[33, 34]. This internal fluorescence can be used to
locate the bead centers within a volumetric scan of the 3D assay. Each bead in
the volumetric scan should be cropped into individual max intensity
projections using the centroid of the bead as measured by the internal
fluorescence channel. Once the bead has been located and cropped, each bead
image should then be analyzed in the fluorescence channel of the respective
detection antibody using the $MFI$ of a circular masked region with dimensions
equal to the diameter of the bead in real-world units. The uncropped inverse
space of the full volume can then be used to develop a background fluorescence
intensity used for background subtraction. Positive fluorescence intensity for
any bead can be interpreted as detection of the respective biomarker in that
location.

##  Appendix B Calibration of Quantitative Biomarker Detection Assay

###  Appendix B.1 Calibration curve generation

This assay utilizes a bead fluorescence intensity measurement using a
fluorescence microscope and assumes the intensity is proportional to the local
concentration of biomarkers due to one or more point sources in a 3D cell
culture assay. This relationship between bead intensity and local biomarker
concentration must be experimentally calibrated. To accomplish this, we have
created a PDMS chip featuring six wells, and a reservoir $+$ diffusion chamber
which is a rectangular well separated by a filter membrane (as shown in Fig.
2(A,B)). The six wells are used to create a standard curve by creating a
serial dilution of the biomarker standard. Detection antibodies and beads are
added to the standards and the beads are allowed to sink to the bottom of each
well. The detection antibody associates with the biomarker in solution and the
beads are allowed to reach an adsorption equilibrium for their respective
concentrations in a diffusion-limited reaction.

To measure the assay kinetics within a 3D culture system, the diffusion
chamber is first filled with the 3D culture medium mixed with capture beads
and detection antibody distributed in a manner described in sections 3.2 and
3.3. In the reservoir, a solution of the highest concentration standard mixed
with detection antibodies is filled to the same level as the 3D culture
medium. The sample is placed onto a fluorescence microscope and the beads
within the standard wells are used to calibrate the microscope intensity
settings such that the dynamic range is maximized for the lowest and highest
bead intensities used in the standards. Once a dynamic range is established,
the bead $MFI$ should be measured according to the protocol outlined above to
establish the standard curve. The detection settings (laser power, gain, lens,
zoom, etc.) should be recorded and kept constant for all future measurements
referring to this standard calibration. The standard concentration versus
measured bead $MFI$ should be analyzed according to the manufacturer’s
instructions to create a calibration curve. This calibration curve should then
be used to transform the intensity measurements from the assay into biomarker
concentration.

Beads within the diffusion chamber are measured every two hours until an
intensity gradient fills the view. The resulting concentration profile is then
compared to the analytical function for flat wall diffusion given by Crank[23]
in eq.2.

$C\left(x,t\right)=C_{o}\;erfc\left[\frac{x}{2\sqrt{Dt}}\right]$ (2)

Where $C(x,t)$ is the time-varying concentration at some distance $x$ from the
flat wall of constant concentration $C_{0}$ and $D$ is the diffusion
coefficient for the biomarker + detection antibody. This diffusion coefficient
can be measured using the protocol outlined in section Appendix C.1. If the
measured concentration profile does not match the analytical function, the
bead density must be decreased. Refer to section Appendix B.2 for more
information related to bead density.

###  Appendix B.2 Determination of ideal bead density

Langmuir kinetics can be used to approximate the sequestration of biomarkers
by each bead within the 3D assay with the number of sequestered biomarkers
reaching an equilibrium surface concentration ($\theta_{eq}$ ) for some given
local concentration ($C$) according to eq.3.

$\theta_{eq}=\frac{K_{eq}C}{1+K_{eq}C}$ (3)

The equilibrium constant $K_{eq}$ describes the binding affinity of the
capture antibody on the bead surface to the biomarker. This equilibrium
constant is unique to each assay and can be determined using eq.8 after the
kinetic variables have been found (refer to Appendix C for more information).
To find the maximum bead density for each 3D assay, we consider the ratio of
the total number of surface sites occupied at equilibrium ($q_{eq}$) versus
the number of biomarkers ($q_{V}$) within some local volume ($V$) surrounding
the bead for a given local concentration ($C$). This sequestration ratio can
be found using eq.4.

$\Gamma=\frac{q_{eq}}{q_{v}}=\frac{N\theta_{eq}}{CV}$ (4)

Where $N$ is the number of capture sites on the bead, and $\theta_{eq}$ is the
surface coverage at equilibrium. When the ratio given in eq.4 is much less
than one, the sequestration of biomarkers from a bead does not significantly
affect the local concentration of neighboring beads. In this case, the maximum
bead density can be calculated using eq.5.

$\rho=\frac{1}{V}$ (5)

##  Appendix C Determination of Bead Assay Kinetic Variables

Many of the variables used in Appendix B are considered proprietary
information by the manufacturer of the bead assay and thus are not readily
available. However, these variables can be found through experimentation. The
following sections outline the procedures for estimating the kinetic variables
for any commercially available assay in a diffusion-limited system.

###  Appendix C.1 Determination of detection antibody + biomarker diffusion
coefficient

A solution of detection antibody mixed with a 1:1 ratio of biomarker is added
to the gel solution used in the 3D culture platform then set within a glass-
bottom culture dish. Fluorescence recovery after photobleaching (FRAPS) is
then used to determine the diffusion coefficient according to the protocol
outlined in Deschout et. al[35].

###  Appendix C.2 Determination of the reaction rates

The detection antibody and biomarker both attach and detach from the capture
bead at a forward reaction rate $k_{f}$ and reverse reaction rate $k_{r}$. The
forward reaction rate can be approximated using eq.6 adapted from the
Einstein-Smoluchowski equation describing the diffusional association of two
reactive species in solution[36, 37]. This forward reaction rate describes the
total number of analyte impinging on the surface of the entire bead normalized
by the number of capture sites and describes the rate the analyte will impinge
on any given capture site.

$k_{f}=\frac{4\pi DR}{N}$ (6)

where $R$ is the radius of the capture bead, $N$ is the number of capture
sites on the bead, and $D$ is the detection antibody + biomarker diffusion
coefficient. The reverse reaction rate should be determined experimentally
since this rate is determined by the antibody-antigen binding affinity. Begin
by transferring a bead that has reached equilibrium surface coverage of
biomarker and detection antibody 111Use beads from the highest concentration
standard when measuring off-rate kinetics to a buffer solution with a zero
concentration of biomarker and detection antibody. The fluorescence intensity
of the bead should then be measured every 5 minutes for 20 minutes. The
resulting $MFI$ measured at some time $t$ can be used to find the reverse
reaction rate $k_{r}$ using eq.7.

$k_{r}=-\frac{1}{t}\;ln\left(\frac{MFI-b}{MFI_{0}}\right)$ (7)

Where $MFI_{0}$ is the initial mean fluorescence intensity and $b$ is the
background fluorescence intensity. This exponential decay is predicted by the
Langmuir adsorption model[38]. The equilibrium constant $K_{eq}$ is simply the
ratio of the forward and reverse reaction rate as seen in eq.8

$K_{eq}=\frac{k_{f}}{k_{r}}=\frac{4\pi DR}{Nk_{r}}$ (8)

###  Appendix C.3 Estimation of the number of capture sites $N$ and
equilibrium surface coverage $\theta_{eq}$

Kinematic variables for the Langmuir equilibrium surface coverage
$\theta_{eq}$ versus local biomarker concentration $C$ can be found using the
$MFI$ versus $C$ data provided by the standard curve. First, the biomarker
concentration versus $MFI$ profile is cropped in the linear range. This allows
for the assumption that the $MFI$ normalized by the maximum reported $MFI$
will be proportional to $\theta_{eq}$ normalized by the maximum equilibrium
surface coverage at the highest biomarker concentration in the linear region
max according to eq.9

$\frac{MFI}{MFI_{max}}=\frac{\theta_{eq}}{\theta_{max}}$ (9)

Substituting eq.8 into eq.3 reveals the relationship between $\theta_{eq}$ and
the assay kinematic variables for diffusion-limited reactions.

$\theta_{eq}=\frac{4\pi DRC}{Nk_{r}+4\pi DRC}$ (10)

Finally, substituting eq.10 into eq.9 allows the reported $MFI$ vs
concentration profile within its linear range to be used to fit eq.11 using
the number of capture sites on each bead $N$ and $\theta_{max}$ as free
variables.

$\frac{MFI}{MFI_{max}}=\frac{4\pi DRC}{\theta_{max}\left(Nk_{r}+4\pi
DRC\right)}$ (11)

###  Appendix C.4 Estimation of detection antibody concentration

The concentration of the detection antibody should be high enough to ensure
all biomarkers generated within the 3D culture assay are captured by the
detection antibody before reaching the embedded capture bead. However, the
background fluorescence signal increases with increasing detection antibody
concentration. We, therefore, utilize the calibration device to determine the
detection antibody concentration in the following manner. The diffusion
chamber from the calibration device is filled with a PBS 1X buffer solution
and the detection antibody solution is added to one side of the chamber
separated by the membrane filter. The fluorescence intensity of the second
chamber is continuously measured every 5 minutes until an intensity gradient
fills the image of a fluorescent microscope. Assuming the $MFI$ of each pixel
is directly proportional to the antibody concentration, the flat wall
diffusion equation from eq.2 can be used to fit the intensity versus distance
profile from the filter membrane using the concentration as a free variable.
The diffusion coefficient of the detection antibody can be found using the
FRAP method described in Appendix C.1.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

## Acknowledgements

The authors would like to thank Dr. Parker Gibbs for generously donating
OS521Oct4 cell lines to this work.

## References

  * [1] Dennis E Discher, David J Mooney, and Peter W Zandstra.  Growth factors, matrices, and forces combine and control stem cells.  Science, 324(5935):1673–1677, 2009.
  * [2] Douglas Hanahan and Robert A Weinberg.  Hallmarks of cancer: the next generation.  Cell, 144(5):646–674, March 2011.
  * [3] Doris P Tabassum and Kornelia Polyak.  Tumorigenesis: it takes a village.  Nature reviews. Cancer, 15(8):473–483, July 2015.
  * [4] Pierre-Alexandre Vidi, Mina J Bissell, and Sophie A Lelièvre.  Three-Dimensional culture of human breast epithelial cells: The how and the why.  In Scott H Randell and M Leslie Fulcher, editors, Epithelial Cell Culture Protocols: Second Edition, pages 193–219. Humana Press, Totowa, NJ, 2013.
  * [5] Jarno Drost and Hans Clevers.  Organoids in cancer research.  Nature reviews. Cancer, 18(7):407–418, April 2018.
  * [6] Bauer L LeSavage, Riley A Suhar, Nicolas Broguiere, Matthias P Lutolf, and Sarah C Heilshorn.  Next-generation cancer organoids.  Nature materials, 21(2):143–159, February 2022.
  * [7] Karen L Cox, Viswanath Devanarayan, Aidas Kriauciunas, Joseph Manetta, Chahrzad Montrose, and Sitta Sittampalam.  Immunoassay Methods.  Eli Lilly & Company and the National Center for Advancing Translational Sciences, July 2019.
  * [8] Sylvia Janetzki, Markus Rueger, and Tomas Dillenbeck.  Stepping up ELISpot: Multi-Level analysis in FluoroSpot assays.  Cells, 3(4):1102–1115, November 2014.
  * [9] Suleyman Aydin.  A short history, principles, and types of ELISA, and our laboratory experience with peptide/protein analyses using ELISA.  Peptides, 72:4–15, October 2015.
  * [10] Joel Fleury Djoba Siawaya, Teri Roberts, Chantal Babb, Gillian Black, Hawa Jande Golakai, Kim Stanley, Nchinya Bennedict Bapela, Eileen Hoal, Shreemanta Parida, Paul van Helden, and Gerhard Walzl.  An evaluation of commercial fluorescent bead-based luminex cytokine assays.  PloS one, 3(7):e2535, July 2008.
  * [11] R J Fulton, R L McDade, P L Smith, L J Kienker, and J R Kettman, Jr.  Advanced multiplexed analysis with the FlowMetrix system.  Clinical chemistry, 43(9):1749–1756, September 1997.
  * [12] Edward Morgan, Rudi Varro, Homero Sepulveda, Julia A Ember, John Apgar, Jerry Wilson, Larry Lowe, Roy Chen, Lalita Shivraj, Anissa Agadir, Roberto Campos, David Ernst, and Amitabh Gaur.  Cytometric bead array: a multiplexed assay platform with applications in various areas of biology.  Clinical immunology, 110(3):252–266, March 2004.
  * [13] Tuomas Huovinen, Laurens Lindenburg, Ralph Minter, and Florian Hollfelder.  Multiplexed affinity characterization of protein binders directly from a crude cell lysate by covalent capture on suspension bead arrays.  Analytical chemistry, 93(4):2166–2173, February 2021.
  * [14] Tapomoy Bhattacharjee, Steven M Zehnder, Kyle G Rowe, Suhani Jain, Ryan M Nixon, W Gregory Sawyer, and Thomas E Angelini.  Writing in the granular gel medium.  Science advances, 1(8):e1500655, September 2015.
  * [15] Maddaly Ravi, V Paramesh, S R Kaviya, E Anuradha, and F D Paul Solomon.  3D cell culture systems: advantages and applications.  Journal of cellular physiology, 230(1):16–26, January 2015.
  * [16] Mark Fricker, John Runions, and Ian Moore.  Quantitative fluorescence microscopy: from art to science.  Annual review of plant biology, 57:79–107, 2006.
  * [17] James Jonkman, Claire M Brown, Graham D Wright, Kurt I Anderson, and Alison J North.  Tutorial: guidance for quantitative confocal microscopy.  Nature protocols, 15(5):1585–1611, May 2020.
  * [18] Catherine A Royer.  Characterizing proteins in their cellular environment: Examples of recent advances in quantitative fluorescence microscopy.  Protein science: a publication of the Protein Society, 28(7):1210–1221, July 2019.
  * [19] Duy T Nguyen, Jack E Famiglietti, Ryan A Smolchek, Zadia Dupee, Nickolas Diodati, Diego I Pedro, Juan M Urueña, Matthew A Schaller, and W Gregory Sawyer.  3D in vitro platform for cell and explant culture in liquid-like solids.  Cells, 11(6):967, March 2022.
  * [20] Padraic P Levings, Sean V McGarry, Thomas P Currie, David M Nickerson, Steven McClellan, Steven C Ghivizzani, Dennis A Steindler, and C Parker Gibbs.  Expression of an exogenous human oct-4 promoter identifies tumor-initiating cells in osteosarcoma.  Cancer research, 69(14):5648–5655, July 2009.
  * [21] Tapomoy Bhattacharjee, Carmen J Gil, Samantha L Marshall, Juan M Urueña, Christopher S O’Bryan, Matt Carstens, Benjamin Keselowsky, Glyn D Palmer, Steve Ghivizzani, C Parker Gibbs, W Gregory Sawyer, and Thomas E Angelini.  Liquid-like solids support cells in 3D.  ACS biomaterials science & engineering, 2(10):1787–1795, October 2016.
  * [22] Amy C Gross, Hakan Cam, Doris A Phelps, Amanda J Saraf, Hemant K Bid, Maren Cam, Cheryl A London, Sarah A Winget, Michael A Arnold, Laura Brandolini, Xiaokui Mo, John M Hinckley, Peter J Houghton, and Ryan D Roberts.  IL-6 and CXCL8 mediate osteosarcoma-lung interactions critical to metastasis.  JCI insight, 3(16), August 2018.
  * [23] J Crank.  The mathematics of diffusion.  Oxford science publications. Oxford University Press, London, England, 2 edition, March 1979.
  * [24] Sharmin Alhaque, Michael Themis, and Hassan Rashidi.  Three-dimensional cell culture: from evolution to revolution.  Philosophical transactions of the Royal Society of London. Series B, Biological sciences, 373(1750), July 2018.
  * [25] Tyler Jacks and Robert A Weinberg.  Taking the study of cancer cell survival to a new dimension.  Cell, 111(7):923–925, December 2002.
  * [26] Donglai Lv, Zongtao Hu, Lin Lu, Husheng Lu, and Xiuli Xu.  Three‑dimensional cell culture: A powerful tool in tumor research and drug discovery (review).  Oncology letters, 14(6):6999–7010, December 2017.
  * [27] Srikkanth Balasubramanian, Kui Yu, Diana Vasquez Cardenas, Marie-Eve Aubin-Tam, and Anne S Meyer.  Emergent biological endurance depends on extracellular matrix composition of Three-Dimensionally printed escherichia coli biofilms.  ACS synthetic biology, 10(11):2997–3008, November 2021.
  * [28] Mototsugu Eiraku, Nozomu Takata, Hiroki Ishibashi, Masako Kawada, Eriko Sakakura, Satoru Okuda, Kiyotoshi Sekiguchi, Taiji Adachi, and Yoshiki Sasai.  Self-organizing optic-cup morphogenesis in three-dimensional culture.  Nature, 472(7341):51–56, April 2011.
  * [29] E L Jackson and H Lu.  Three-dimensional models for studying development and disease: moving on from organisms to organs-on-a-chip and organoids.  Integrative Biology, 8(6):672–683, May 2016.
  * [30] Masamitsu Oshima, Miho Ogawa, and Takashi Tsuji.  Regeneration of complex oral organs using 3D cell organization technology.  Current opinion in cell biology, 49:84–90, December 2017.
  * [31] J P Hughes, S Rees, S B Kalindjian, and K L Philpott.  Principles of early drug discovery.  British journal of pharmacology, 162(6):1239–1249, March 2011.
  * [32] Carrie J Lovitt, Todd B Shelper, and Vicky M Avery.  Advanced cell culture techniques for cancer drug discovery.  Biology, 3(2):345–367, May 2014.
  * [33] Edmond J Breen, Veronika Polaskova, and Alamgir Khan.  Bead-based multiplex immuno-assays for cytokines, chemokines, growth factors and other analytes: median fluorescence intensities versus their derived absolute concentration values for statistical analysis.  Cytokine, 71(2):188–198, February 2015.
  * [34] David Holmes, Joseph K She, Peter L Roach, and Hywel Morgan.  Bead-based immunoassays using a micro-chip flow cytometer.  Lab on a chip, 7(8):1048–1056, August 2007.
  * [35] Hendrik Deschout, Joel Hagman, Sophia Fransson, Jenny Jonasson, Mats Rudemo, Niklas Lorén, and Kevin Braeckmans.  Straightforward FRAP for quantitative diffusion measurements with a laser scanning microscope.  Optics express, 18(22):22886–22905, October 2010.
  * [36] Otto G Berg and Peter H von Hippel.  Diffusion-Controlled macromolecular interactions.  Annual review of biophysics and biophysical chemistry, 14(1):131–158, June 1985.
  * [37] M v Smoluchowski.  Grundriß der koagulationskinetik kolloider lösungen.  Kolloid-Zeitschrift, 21(3):98–104, September 1917.
  * [38] Irving Langmuir.  The adsorption of gases on plane surfaces of glass, mica, and platinum.  Journal of the American Chemical Society, 40(9):1361–1403, September 1918.