\section{Introduction}
\label{sec:intro}
Robotic computing is on the rise. A myriad of robots such as drones, legged robots, and self-driving cars are on the verge of becoming an integral part of our life~\cite{wan2021survey,liu2021robotic}. Robotics is typically an art of system integration both in software and hardware (Fig.~\ref{fig:overview}). The continuous proliferation of robots, however, face computing challenges, raised from the higher performance requirements, resource constraints, miniaturization of machine form factors, dynamic operating scenarios, and cybersecurity considerations. Therefore, it is essential to choose a proper computing substrate for robotic system that can meet real-time and power requirements and adapt to changing workloads.

CPUs and GPUs are two widely-used computing platforms, however, their performance and efficiency are still incompetent in real-time computation for complex robots.
Take the motion planning task as an example, CPU typically takes a few seconds to find the collision-free trajectory~\cite{hauser2015lazy}, making it too slow for complex navigation tasks. GPUs can finish planning tasks in hundreds of milliseconds, still insufficient for many scenarios while at hundreds of watts cost~\cite{pan2012gpu}. ASICs are recently developed for specific robotic workloads with low power and high performance~\cite{suleiman2019navion,yoon2020neuroslam,wan2022circuit}, but their fixed architecture has difficulty in adapting to rapid-evolving robotic algorithms and dynamic scenarios, and is vulnerable to cybersecurity threats.

As an alternative, we believe FPGA is the promising compute substrate for robotic applications. First, FPGA increases the performance with massive parallelism and deeply pipelined datapath, making it capable of meeting real-time requirements with high energy efficiency compared to CPUs and GPUs. Second, FPGA can adaptively generate custom architectures and update with the fast-evolving of robotic algorithms without going through re-fabrication as ASIC~\cite{mayoral2021adaptive}. Third, FPGA is flexible in dealing with highly diverse robotic workloads, especially with partial reconfiguration allowing modification part of the operating board. Fourth, FPGA provides reliable design by leveraging reconfiguration to patch flows, compared to potential vulnerabilities detected in fixed architectures~\cite{kocher2019spectre}, which is especially essential in safety-critical scenarios~\cite{wan2021analyzing}. Overall, FPGA has the potential to deliver high-performance, low-power, reconfigurable, adaptive, and secure features in robotic computing, and is booming in autonomous applications. However, several challenges, such as tedious development procedures, inefficient system support, and huge design space, remain in the FPGA-based robotic computing and impede the way ahead.

In this paper, we will discuss the current progress, challenges, and opportunities for FPGA-based robotic computing. Section~\ref{sec:algo} introduces the cross-layer stack of robotic system. Section~\ref{sec:FPGA} presents current FPGA accelerators and systems for robotic computing, with an emphasis on design techniques. Section~\ref{sec:challenges} discusses challenges and opportunities for FPGA-based robotic computing, and our view of the road ahead.

\section{Cross-Layer Robotic Computing Systems}
\label{sec:algo}
This section introduces the abstraction layers of the robotic computing stack. We traverse down Fig.~\ref{fig:overview} to explain robotic-specific algorithms and systems building blocks.

\begin{figure}[t!]
        \centering\includegraphics[width=\columnwidth]{figs/overview_v3.pdf}
        \caption{Cross-layer stack of the robotic computing system, requirements, research challenges, and open opportunities.}
        \label{fig:overview}
        \vspace{-10pt}
\end{figure}

\subsection{Robotic-Computing Algorithm Layer}
Fig.~\ref{fig:algorithm} illustrates the representative algorithm building blocks in robotic computing, including sense-plan-act (perception, localization, planning, control) and end-to-end learning.

\textbf{Perception.}
The goal of perception is to sense the dynamic surroundings and build a reliable and detailed representation based on sensory data (e.g., camera, IMU, GPS, LiDAR). Perception usually includes feature extraction, stereo vision, object detection, scene understanding, etc. In feature extraction, key points are usually detected using FAST feature and ORB descriptor~\cite{rublee2011orb}. Compared with all image pixels, operating on feature points can improve the robustness and compute efficiency. Stereo vision is to obtain 3D structure information of the scene through disparity calculation. Local, semi-global, and global stereo matching algorithms are proposed based on operational scenarios~\cite{lu2021resource}. Recently, advances in deep learning have exposed robotic perception systems to more tasks.

\textbf{Localization.}
The goal of localization is to calculate the position and orientation of a robot itself in a given frame of reference. Knowing the position fundamentally enables robots to plan the trajectory and navigate, and knowing the orientation further helps robots stabilize. Simultaneous localization and mapping (SLAM) is a commonly-used algorithm where the robot simultaneously constructs a map of the environment while localizing itself~\cite{mur2017orb}, and one principled mathematical approach to solving SLAM is maximum a posteriori estimation. The filtering-based approach has recently been developed with Multi-State Kalman Filter-based algorithms such as MSCKF VIO~\cite{sun2018robust} and OpenVINS~\cite{geneva2020openvins}.

\textbf{Motion planning and control.}
The goal of motion planning is to find the optimal collision-free trajectory from the start position to the goal position, which is invoked during a robot movement to adapt to environmental changes. Motion planning is usually followed by a control module continuously tracking the differences between actual poses and poses on the pre-defined trajectory. Sampling-based solutions are widely used for motion planning, such as Probabilistic Roadmap (PRM)~\cite{ichter2020learned}, Rapidly-exploring Random Tree (RRT)~\cite{lavalle2001rapidly} and their variants, which generally contain three steps: roadmap construction, collision detection, and graph search.

\textbf{End-to-End learning system.}
End-to-end algorithms enable skill learning directly from sensor input and perform all the following cognitive robotic tasks using a single neural network model. Maps or separate planning stages are not required in end-to-end learning. The neural network model can be trained using reinforcement learning~\cite{anwar2020autonomous} or supervised learning~\cite{loquercio2018dronet}. The challenges of end-to-end learning include alleviating the model simulation-to-reality performance gap, designing optimal reward functions, and improving model explainability and robustness, which are actively explored.

\subsection{Robotic-Computing System Layer}
\textbf{Robot Operating System (ROS).}
ROS is a commonly used operating system to provide tools, libraries, and package management for robotics development.
It is a distributed framework of processes that enables executables to be individually designed and loosely coupled at runtime.
Conceptually, the peer-to-peer network of ROS processes is called computation graph. The basic ROS computation graph includes nodes, topics, services, and masters, all of which provide data to the graph in different ways (Fig.~\ref{fig:overview}). Each ROS node is a process used to perform a task. ROS nodes communicate with each other via topics or services. Topics allow one node to publish messages that multiple other nodes can subscribe. Services allow for creating a one-to-one communication between a service node and a client node.
The ROS master is responsible for storing operating parameters and managing other nodes.

\begin{figure}[t!]
        \centering\includegraphics[width=\columnwidth]{figs/algorithm_v2.pdf}
        \caption{Applications and algorithm building blocks in robotic systems.}
        \label{fig:algorithm}
        \vspace{-10pt}
\end{figure}

\section{Current Progress and Design Techniques}
\label{sec:FPGA}
This section presents our designs and current progress for FPGA-based robotic computing, with an emphasis on the design traits and techniques.

\textbf{Perception on FPGAs.}
The perception typically contributes significantly to the end-to-end latency of robotic applications. Take the ORB perception module as an example, it usually accounts for 50\%-80\% compute latency of the whole localization scenario. To alleviate that, \cite{wan2021energy,gan2021eudoxus} accelerate ORB-based perception on FPGA for both aerial and ground robots. The key design principles are to exploit task-level parallelisms by frame-multiplexing feature extraction, and customize on-chip memories to suit different types of data reuse. Another series of illustrative works is to accelerate stereo matching, which is the bottleneck of the stereo vision system. \cite{perri2020stereo} implements local stereo matching algorithms on FPGA with characterized hardware-software partitions. \cite{kamasaka2018fpga} proposes a parallel 3D graph cut algorithm for accelerating global stereo matching, achieving 166$\times$ speedup to CPU. FP-Stereo~\cite{zhao2020fp} present streaming architecture and sampling-insensitive disparity algorithm on FPGA to accelerate semi-global stereo matching. Recently, the Bayesian approach with generative probabilistic models facilitates efficient dense matching. An appealing example is iELAS~\cite{gao2021ielas}, a hardware-friendly large-scale stereo algorithm implemented on FPGA. iELAS reforms the computational-intensive and irregular triangulation modules in a regular manner with intelligent points interpolation. Additionally, FPGA has been widely used in accelerating neural networks for robotic perception and end-to-end learning. Several techniques, including quantization, loop optimization, array partitioning, data reuse, and memory optimization, are proposed. Interested readers are pointed to~\cite{abdelouahab2018accelerating} for details.

\textbf{Localization on FPGAs.}
The backbone of SLAM is a complex non-linear optimization problem, bundle adjustment (BA), which consumes a significant amount of time and power. $\pi$-BA~\cite{liu2020pi} designs a co-observation optimization technique to accelerate BA based on the key inspection that not all 3D points appear on all images in a BA problem. Pisces~\cite{asgari2020pisces} co-optimizes SLAM power consumption and latency by exploiting inherent SLAM sparsity. By orchestrating sparse data, Pisces aligns correlated data and enables direct, parallel, and deterministic memory access. Going beyond point solutions, Archytas~\cite{liu2021archytas} presents a template hardware synthesis solution that automatically generates a SLAM accelerator given the hardware template and algorithm data-flow graph. To make the design adaptable to various environments, \cite{liu2022energy} dynamically optimizes the SLAM accelerator with an offline constructed lookup table and clock gating without sending new bitstreams to FPGA at run time. Typically, SLAM is suitable for unknown indoor environments, while Registration is used for known indoor environments and visual-inertial-odometry (VIO) functions well for outdoor environments. No single localization algorithm fits all scenarios. Interestingly, these algorithms share fundamental computation kernels amenable to matrix blocking. Eudoxus~\cite{gan2021eudoxus} implements SLAM, Registration, and VIO on FPGA by accelerating common matrix operations, with a lightweight runtime scheduler reducing variation.

\textbf{Motion planning and control on FPGAs.}
Among the motion planning pipeline, the computation of collision detection is usually the bottleneck. Take RRT as an example, when it runs on CPU, 99\% of the instructions are executed for collision detection, taking up 90\% of total computation time. Recent efforts have proposed to accelerate motion planning kernels through algorithm-hardware co-design on FPGAs. \cite{murray2016microarchitecture} constructs robot-specific circuitry and architecture with roadmap pre-computation and massive path search parallelism, which is able to solve a motion planning query in 16~$\mu$s. \cite{murray2019programmable} further presents a programmable dataflow architecture with a low-cost interconnection network, reducing the latency to 2.3~$\mu$s.

Several key design and optimization techniques are leveraged in FPGA-accelerated perception, localization, planning, and control. From the software aspect, hardware-friendly algorithms and data structures are proposed to promote parallelism with reduced intrinsic recursions and algorithm complexity. From the hardware aspect, robot-specific architecture, data sparsity, locality, parallelism, optimized interconnection networks, and reduced data movement contribute to high-performant and flexible motion planning design. These techniques can be generalized to other implementations, serving as a guide for future works.

\textbf{Multi-robot collaboration on FPGAs.}
Going beyond single-robot applications, swarm robotics has been increasingly deployed in real-life scenarios where a team of robots collaboratively finish a task. Multi-robot workload typically demonstrates unique compute challenges. Several algorithm kernels may need to process the data at the same time, leading to hardware resources conflicts. Therefore, the FPGA accelerator should support multi-thread and dynamic scheduling. An intriguing example is INCAME~\cite{yu2021incame}, a single-core multi-robot exploration framework that supports dynamic multi-task scheduling with a virtual-instruction-based interrupt method. The perception and control tasks are assigned high priorities, while long-term decisions and optimization have low priorities. We envision that the multi-core multi-tasking FPGA accelerator will further improve the performance of the multi-robot system.

\textbf{ROS on FPGAs.}
ROS-compliant FPGAs have recently been developed with ROS becoming increasingly common in robotics. ROS-compliant FPGAs must consider four functions: encapsulation of FPGA circuits, interface between ROS software and FPGA hardware, subscribe interface, and publish interface to ROS topic. Typically, large communication latency between ROS components is the bottleneck of offloading computing to FPGAs. \cite{sugata2017acceleration} reduces the latency by implementing publish and subscribe messaging of ROS as hardware circuits, making the direct ROS-FPGA communication possible and efficient. Recently, \cite{leal2020automated,lienen2020reconros} propose tools and frameworks to offload and accelerate ROS computational graph on FPGA. However, the ecosystem of ROS on FPGAs is still in its infancy, better interface, automated tools, and whole ROS acceleration are to be developed.

\section{Research Challenges and Future Directions}
\label{sec:challenges}
This section discusses the research opportunities for FPGA-based robotic computing, and our view for road ahead (Fig.~\ref{fig:overview}).

\circled{1} \textbf{Reconfiguring robotic computing at run time.}
Robots usually operate in highly dynamic environments, thus designing runtime-reconfigurable compute platforms is critical and can enable robots to be adaptive in various scenarios.
Partial reconfiguration (PR) is a key feature of FPGA. Using PR, part of FPGA can be reconfigured at runtime without compromising the integrity of the applications running on those parts of the device that are not being reconfigured.
Therefore, PR can allow various robotic computing kernels to time-share part of an FPGA, leading to high performance and energy efficiency, and making FPGA a more suitable computing platform for dynamic and complex robotic workloads.

\circled{2} \textbf{Modularizing robotic computing kernels design.}
The number of robotic algorithms is booming, but many algorithm variants share similar key computation blocks. It is thereby imperative to modularize the robotic computing kernel design.
We can build optimized hardware acceleration blocks for these kernels as libraries or packages, while exploring their inherent task-specific features such as sparsity, data flow, and memory access patterns. During the design phase, robotics practitioners can directly import these robotics-specific libraries and building blocks to build their FPGA design without delving into hardware engineering, which will greatly ease the design process.
Modularizing the robotic algorithm design can help roboticists create custom accelerators for a kernel without hardware expertise.

\circled{3} \textbf{Mapping robotic computing on heterogeneous platforms.}
One of the key technical challenges of designing robotic compute systems is to develop a suitable computer architecture, along with a software stack that allows computational flexibility.
To improve the overall performance, FPGA-based System-on-Chip (SoC) solutions for robotic computing would be of the essence~\cite{mayoral2021adaptive}, which holistically integrates various computing technologies, including CPU, GPU, FPGA, and accelerators.
The OpenCL framework can be used for programming and executing programs across heterogeneous platforms, and accelerator-level parallelism is expected to be explored~\cite{hill2021accelerator}.
By doing so, the SoCs are equipped with both software and hardware programmability, having the capability to deliver high performance, low power, adaptive, and reliable robotic computing.

\circled{4} \textbf{Connecting FPGA to ROS ecosystem.}
With ROS increasingly utilized in robotics applications of all scales, robotic FPGA platforms need to be able to efficiently map ROS computational graphs on silicon.
Going beyond the current work on accelerating specific ROS libraries, the inter-process and intra-process between ROS nodes also need to be accelerated~\cite{lienen2022reconros}. It is worth noting that the hardware acceleration must be directly integrated into the ROS ecosystem to provide a seamless user experience for roboticists. A better interface between ROS and FPGA is expected to be delivered.
Furthermore, through dynamically and efficiently mapping ROS to heterogeneous compute platforms, holistic hardware acceleration for robotic computing on ROS applications is expected to be achieved.

\circled{5} \textbf{Benchmarking robotic computing kernels.}
Given the proliferation of robotic kernels and the rapid advances of hardware platforms, benchmarking these robotic algorithms and systems in a comparable, quantitative, and validatable manner is imperative.
Such benchmarking comes into two folds, benchmarking a robotic algorithm across various hardware platforms, and benchmarking various robotic algorithms within the same hardware~\cite{neuman2019benchmarking}. Particularly, benchmarks should consider the interactions of ROS and its computational graph.
Benchmarking robotic computing will guide the robotics and hardware researchers to investigate the trade-offs in accuracy, performance, and energy efficiency of various robotic algorithms, and implement (or select) algorithms on FPGAs and other platforms in a performance-portable way.

\circled{6} \textbf{Automating robotic computing design flow.}
Given the increasing complexity of robotic algorithms and the cross-stack nature, the development of robotic computing systems is becoming slow and tedious.
Thus, building a push-button flow with robotic task requirements as input to automatically generate robotic accelerator design is critical~\cite{krishnan2021autopilot,neuman2021robomorphic,krishnan2021autosoc}.
We envision the agile framework will intelligently search the huge design space and automatically choose the optimal algorithm-hardware parameters with the help of modular kernels, benchmarking, and machine learning-assist methods. New robotic-centric electronic design automation (EDA) tool needs to be developed to convert the design to FPGA implementation.
Automating the design follow will greatly facilitate the FPGA-based robotic computing development, and make FPGAs an ideal platform for fast prototyping and commercialization.

\circled{7} \textbf{Building customized robotic computing with the open-source framework.}
The field of robotic computing is still in its infancy and fast-changing, and numerous opportunities still exist in task-specific acceleration.
The open-source design framework with iteratively deployment, profiling, and optimization has recently been developed for machine learning applications~\cite{prakash2022cfu}, but it is still to be explored for robotic computing applications. Designers can build their custom specialized and optimized processors based on the RISC-V instruction set architecture (ISA).
Defining and building an open-source FPGA-based RISC-V robotics-on-chip processor with open-source frameworks would considerably facilitate the design process and allow us to adapt to the rapidly changing landscape of robotic computing algorithms and accelerators.

\circled{8} \textbf{Integrating robotic computing hardware in a simulation loop.}
The FPGA-accelerated kernels are usually part of the whole autonomy computing pipeline. The correlation among compute stages and other robotic cyber-physical components will impact the final robotic system performance and lead to inaccurate hardware evaluation~\cite{krishnan2020sky,krishnan2022roofline}. Thus, instead of isolated hardware development, adopting the hardware-in-the-loop (HIL) method is critical~\cite{boroujerdian2018mavbench}. HIL requires plugging the hardware platforms into the simulation to understand how robots respond to stimuli on FPGA or other compute substrates. HIL can help designers quantify the FPGA real-time performance within the whole system and enable robust evaluation without risking real robots. Particularly, HIL can alleviate the FPGA hardware-induced gaps between training and deployment in learning-based systems. To perform faster performance evaluation at an earlier design stage, a closed-loop co-simulation framework of both FPGA architectural behavior (e.g., FireSim~\cite{karandikar2018firesim}, SystemModeler~\cite{acevedo2016fpga}) and robotic environment simulator (e.g., AirSim~\cite{shah2018airsim}) is necessary.

The abundance of challenges raised above provides plentiful opportunities for research development at all levels. Endeavoring to solve these problems requires interdisciplinary approaches across all layers of computing stack, from algorithm and system to architecture, micro-architecture, and circuits.

\section{Conclusion}
\label{sec:conclusion}
Robotic computing is a rising area and critically depends on efficient, adaptive, and reliable compute substrates. This paper presents the cross-layer robotic computing stack and illustrates the current progress, along with FPGA design techniques.
We conclude the paper by discussing the challenges, research opportunities, and roadmap for the next-generation FPGA-based robotic computing systems.

\section*{Acknowledgements}
This work was supported in part by C-BRIC, one of six centers in JUMP, a Semiconductor Research Corporation (SRC) program sponsored by DARPA.
\section{Introduction}

\bibliographystyle{ieeetr}


\section{Introduction}

Among most therapeutic areas, failures in clinical trials are common and costly. Failure rates of drugs entering Phase I trials have hit 90\% across most therapeutic areas \citep{mullard2016parsing}. Oncology in particular has one of the highest rates of clinical trial failures, with only 4\% \citep{mullard2016parsing} of therapies entering phase 1 FDA clinical trials ultimately being approved, despite having the most active clinical trials, with approximately 32\% of Phase I FDA clinical trials being in oncology \citep{thomas2016clinical}.
In part due to the desire to accelerate new medicines into the clinic to address unmet medical need, but also driven by competition in the industry, drug development organizations do not devote adequate time and resources to improve preclinical model systems that might be more predictive of clinical results \citep{honkala2021harnessing}.  Instead, the standard approach adopted by most of the industry relies on \textit{in vitro} and \textit{in vivo} tumor model systems that are poorly predictive of activity in patients due to the reductionist nature of the systems and inadequate attention devoted to understanding the molecular similarities and differences between preclinical and clinical data.

In this work, we describe an end-to-end machine learning pipeline, SystemMatch, that optimizes a preclinical model to best approximate the behavior of a target \textit{in sapiens} population to enable drug developers to quickly optimize preclinical models and identify those with the greatest predictive power for therapeutic priorities. SystemMatch uses single-cell genomic data from many preclinical models and evaluates them against a multi-study atlas of single-cell genomics data from the tumor or tissue of interest, helping decide which of these models is most likely to provide clinically meaningful predictions.
Furthermore, SystemMatch utilizes a Compositional Perturbational Autoencoder (CPA) \citep{lotfollahi2021compositional} to predict the behavior of single cells in previously untested combinations of experimental conditions. We use these predictions to recommend experimental changes to the preclinical models to enhance their similarity to the target population.

We demonstrate SystemMatch on a large multi-condition perturbational dataset of \textit{in vitro} differentiating macrophages at the single-cell level, and we compare these model systems' proximity to a target population of human tumor suppressive macrophages obtained from a multi-study single-cell macrophage atlas we collated, integrated, and annotated for this purpose. We demonstrate that SystemMatch produces systems that contradict the standard dogma for generating \textit{in vitro} suppressive macrophages, and we recommend further optimizations of our model system to generate systems with more predictive power for this purpose. This is, to our knowledge, the first computational pipeline for assessing and optimizing the fitness of preclinical models to an \textit{in sapiens} target population.

\section{Background}

\subsection{Single-cell omics}

Single cell omics refers to the quantitative characterization of biological phenotype at the cellular level. Early work in single-cell omics focused on the development of single-cell RNA sequencing technology \citep{kolodziejczyk2015technology}, in which the quantity of gene transcripts (or mRNA) in each cell is counted through complex microfluidic assays. Further work has expanded this technology to the measurement of many modalities, including genes, proteins, metabolites, transcripts, lipids and more \citep{chen2020single}, as well as multimodal data such as CITE-seq, which measures both RNA expression and protein abundance in the same assay \citep{citeseq}. As the quality and quantity of single-cell omics data rise, it is increasingly straightforward to precisely define rare cell types and hitherto poorly understood cellular heterogeneity (see, e.g., \citet{jaitin2014massively,zeisel2015cell}). Bulk RNA sequencing, on the other hand, refers to the sequencing of transcripts present in a large number of cells at the resolution of a cell type or tissue. Most large bulk datasets are not cell type specific, which makes it impractical to understand the phenotypic profile of specific cell subpopulations. For example, with bulk RNAseq data obtained from solid tumors, it is typically not possible to accurately understand the transcriptomic profile of specific immune subpopulations, like CD8 memory T-cells or immunosuppressive macrophages, despite recent work to deconvolve bulk data to cell type resolution \citep{newman2015robust,finotello2019molecular}.

\subsection{Preclinical model systems}

Due to both the cost and ethical implications of testing novel drugs in humans, most or all drugs are first tested in preclinical models, which range from microorganisms, to cell- and tissue-based models, to animal models including mice and non-human primates. Preclinical data are typically required for FDA approval to test a drug in humans~\citep{mcelvany2009fda} and are additionally used to prioritize selection of drugs to advance to clinical trial \citep{denayer2014animal}.  However, failures to translate success in preclinical models to humans have cast doubt on the predictive power of these models, prompting some to question their utility in drug prioritization (see, e.g., \citet{Schnabel2008,Suckling2008}), with some even going as far as to recommend forgoing animal models altogether and testing drugs directly in humans \citep{shanks2009animal}. On the other hand, reverse translational approaches seek to inform the development of preclinical models through the study of clinical success, creating an iterative process between preclinical and clinical studies to optimize later generations of drugs, the success of which can be seen, for example, in the development of EGFR tyrosine kinase inhibitors for the treatment of non-small cell lung carcinoma \citep{honkala2021harnessing}.

\subsection{Perturbation prediction}

Predicting cellular responses to perturbations is an important goal in computational biology.  \citet{ji2021machine} detail several uses for modeling perturbational single-cell data, including perturbation response prediction, target and mechanism prediction, perturbation interaction prediction, and chemical property prediction. Here, we introduce a new use case for perturbational modeling, which is to predict perturbations that will generate an optimal model system.  We generate \textit{in silico} predictions for a wide variety of possible perturbations, then select those closest to our target model system.  In this way, we can more rapidly converge on an ideal model system.

\section{Related work}

\subsection{Preclinical model system development}

Classical protocols of preclinical model system design use direct measurement of the preclinical system's phenotype to evaluate the quality of the system. For example, \citet{mia2014optimized} select among a set of possible protocols to generate \textit{in vitro} immunosuppressive macrophages by measuring a) secreted proteins known to be markers of immunosuppression and b) \textit{in vitro} suppression of T cells, finding that \textit{M-CSF} + \textit{IL-4} + \textit{IL-10} + \textit{TGF$\beta$} generates the most suppressive macrophages. \citet{fogg2020ovarian} take a different approach, culturing monocytes with ovarian cancer cell lines in order to understand the pathways activated by the cancer cells, finding an alternative pathway to macrophage polarization through \textit{TGF$\alpha$}. Reverse translational medicine, on the other hand, uses deep characterization of clinical response to existing drugs to understand the mechanisms of action of these drugs in order to design preclinical models that replicate the drug resistance mechanisms in humans \citep{honkala2021harnessing}. However, to our knowledge, no prior work has been done to optimize a preclinical model's phenotype in the high-dimensional space made visible through single-cell genomics.

\subsection{Perturbation response prediction}
Machine learning models have been used in a number of different ways to predict cellular perturbation response \cite{ji2021machine}.  A common setting is \emph{causal imputation} \cite{squires2020causal}, where a model is required to predict a response to an intervention in a particular context after having been trained on related interventions and contexts.  For example, \citet{squires2020causal} predict genomic response data by training on perturbations in some cell types and predicting response in other cell types.  \citet{lotfollahi2019scgen} introduces \textit{scGen} for causal imputation on single-cell gene expression data, which uses a variational autoencoder to represent interventions in latent space and then adds the interventions in latent space to an unperturbed representation to obtain the perturbational response.  In our setting, where we wish to predict perturbations that generate gene expression corresponding to an \textit{in vivo} model system, combinations of treatments are also important to consider.  Recently, the Compositional Perturbational Autoencoder \cite{lotfollahi2021compositional} was demonstrated to be able to predict combinations of gene knockouts from being trained on each knockout individually.

\section{SystemMatch}

Here, we introduce our iterative preclinical model optimization pipeline, SystemMatch. Fig.~\ref{fig:schematic} shows a schematic representation of the SystemMatch process. First, in order to define the \textit{target} population, we collect an atlas of single-cell data from relevant clinical cohorts giving a robust universal representation of the cell type of interest across multiple disease states. We apply a simple reference-based cell type annotation (see, e.g., \citet{seurat,scvi}) to only retain the cells associated with the target population (e.g., all macrophages). Next, we integrate this single-cell atlas using a fixed gene list derived from expert domain knowledge to capture only the biology relevant to the system at hand and use this as input to a batch correction algorithm (see, e.g., \citet{harmony,seurat,haghverdi2018mnn}). Within this integrated single-cell cell type--specific atlas, we annotate subtypes using a combination of graph-based clustering \citep{leiden} and expert domain knowledge. We choose the target population from these subtypes as the subtype that most closely recapitulates the phenotype of the target population as described in the literature (e.g., macrophages that express known markers of immunosuppression).

\begin{figure}[tb]
\centering
\includegraphics[width=\textwidth]{figures/2022-03-15_schematic_small.png}
\caption{Schematic of the iterative SystemMatch pipeline.}
\label{fig:schematic}
\end{figure}

Second, we collect the \textit{query} model systems. We generate this data by developing many different model systems in parallel experiments (here, \textit{in vitro} macrophage differentiation and polarization) with multiple experimental conditions (here, stimulation with different cytokines and combinations thereof) selected from prior domain knowledge to generate a heterogeneous set of model systems. We perform single-cell sequencing on the cells of interest from each model system and then perform quality control on these datasets to remove low quality and outlier cells in order to obtain a clean single-cell resolution representation of each model system, and project these datasets to a subspace of the full gene expression matrix by retaining only those genes known to be related to the biological function of interest.

Next, we compute the distances in this subspace between each query dataset and the target population using some distance metric $m$. In the simplest case, $m$ is simply the L2 distance between the average (or "pseudobulked") gene expression of the query ${\mathbb{X}} := \{{x}_i \in \mathbb{R}^d\}$ and the target ${\mathbb{Y}} := \{{y}_i \in \mathbb{R}^d\}$

$$m_{L2}({\mathbb{X}}, {\mathbb{Y}}) := \left\Vert \frac{\sum_i^{n_{\mathbb{X}}} {x}_i}{n_{\mathbb{X}}} - \frac{\sum_i^{n_{\mathbb{Y}}} {y}_i}{n_{\mathbb{Y}}} \right\Vert_2$$

where $n_{\mathbb{X}}:=|X|$ and $n_{\mathbb{Y}}:=|Y|$. However, other more complex distance metrics between single-cell datasets could also be used, e.g. the Earth Mover's Distance (EMD) or Wasserstein metric \citep{kantorovich1960mathematical}

$$m_{EMD}({\mathbb{X}}, {\mathbb{Y}}) := \min\limits_{\bm{F}} {\sum_{i=1}^{n_{\mathbb{Y}}}\sum_{j=1}^{n_{\mathbb{X}}} {F}_{i,j}d(x_i,x_j)}$$

where ${\bm{F}}$ is the $n_{\mathbb{X}} \times n_{\mathbb{Y}}$ flow matrix with ${F}_{i,j} \ge 0$ and $d$ is a distance metric between cells, e.g. the Euclidean distance. Then, using our choice of $m$, we rank all query datasets to produce an ordered set of queries, with the query least distant from the target denoted the most representative model system of those tested.

However, this still leaves an open question: can we further improve the models beyond the original set of tested systems? To assist with this question, we employ generative deep neural networks to generate all possible combinations of conditions tested in the experimental queries to generate combinatorially many \textit{in silico} model systems. Then, in order to avoid over-reliance on the accuracy of these predictions, we leverage the \textit{in silico} queries to generate a search space of possible experimental conditions that are substantially different from the systems already tested. We use a modified form of the $k$-medoids algorithm \citep{kaufman2009finding} applied to the pseudo-bulked \textit{in silico} queries. Briefly, $k$-medoids selects $k$ equidistant model centroids as in $k$-means but requires that these centroids be selected from the existing data points. We extend this algorithm to enforce that all selected medoids are also equidistant from all existing experimentally tested queries, which leaves us with the smallest possible subset of experiments to run in the next iteration while ensuring that we do not leave any region of the search space untested. We can further combine these recommendations with expert knowledge by examining the genes driving the difference between our best-ranked query and the target population to propose additional experimental conditions.

Finally, we use the outputs of this pipeline to re-run the model system generation experiment, and we iterate upon this process until the generated model system is either sufficiently similar to the target population, or further iterations fail to improve upon the existing queries.

\begin{figure}[tb]
     \centering
     \begin{subfigure}[t]{0.35\textwidth}
         \vskip 0pt
         \centering
         \includegraphics[width=\textwidth]{figures/2022-03-06_in-vivo-umap.png} \\
         (a) \\
         \includegraphics[width=\textwidth]{figures/2022-03-07_in-vitro-umap.png} \\
         (b)
     \end{subfigure}
     \hfill
     \begin{subfigure}[t]{0.17\textwidth}
         \vskip 0pt
         \centering
         \includegraphics[width=\textwidth]{figures/2022-03-06_in-vivo-heatmap.png}
     \end{subfigure}
     \hspace*{-0.9cm}
     \begin{subfigure}[t]{0.3925\textwidth}
         \vskip 0pt
         \centering
         \includegraphics[width=\textwidth]{figures/2022-03-07_in-vitro-heatmap.png}
     \end{subfigure} \\
     \vspace*{-1em}
     \hfill
     \begin{minipage}{0.57\textwidth}
     \centering
     (c)
     \end{minipage}
    \caption{Characterization of human tumor macrophage atlas and \textit{in vitro} stimulated macrophages. a) Human tumor macrophage atlas UMAP shows seven major clusters; b) \textit{in vitro} macrophage UMAP shows \textit{in vitro} conditions separate into three major groups; c) heatmap of major macrophage function genes in human macrophages (left) and \textit{in vitro} macrophages (right) shows the
    \textit{gene36-gene46+} cluster is most strongly associated with known suppressive response genes \textit{gene20} and \textit{gene46}, while this combination is absent in vitro.}
    \label{fig:dataset}
\end{figure}

\subsection{Generating \protect\textit{in silico} queries}

Since the space of possible model systems is extremely large and diverse, we are unable to test all possible experimental conditions \textit{in vitro}. Hence, to better explore the space of untested systems that could best emulate the target population, we predict the gene expression of combinations of the cytokine stimulations in our \textit{in vitro} experiment. We use CPA \cite{lotfollahi2021compositional} to generate \textit{in silico} perturbed cells for 66 additional conditions. Briefly, CPA is an autoencoder trained to decompose the data into disentangled latent representations for three key attributes of each cell: the cell's "basal state", the perturbation effect, and the covariate effect. These latent representations are combined in the decoder and trained via a reconstruction loss. To enforce independence of the latent representations, CPA is trained using a discriminator network and adversarial loss such that no signal from the observed perturbation or covariates is captured in the cellular basal state embedding. As the latent representations are combined by summation, CPA can combine multiple perturbations, allowing us to generate \textit{in silico} samples for combinations of cytokine stimulations. We further modify CPA to allow the incorporation of multiple covariates in order to account for experimental batch covariates included in our \textit{in vitro} dataset.

\subsection{Selecting new queries for experimental iteration}

To select a subset of the generated \textit{in silico} perturbations for experimental validation, we devise a scheme based on $k$-medoids \citep{kaufman2009finding} which ensures that the conditions selected for experimental validation are a) sufficiently different from the conditions already tested; and b) sufficiently heterogeneous to cover the space of possible perturbations. Briefly, we run the regular k-medoids algorithm, but at the medoid-update step of the algorithm, we consider each of the existing \textit{in vitro} perturbations as fixed medoids; this way, the $k$-medoids algorithm selects $k$ new perturbations for testing which are both heterogeneous and far from any of the previously tested conditions.

\section{Results}

\subsection{Characterizing \protect\textit{in vitro} and \protect\textit{in vivo} tumor associated macrophages}

We collated public and proprietary data from tumour-infiltrating immune cells from one study we generated internally and four public studies~\citep{zhang2021single,bassez2021single,yost2019clonal,qian2020pan}, comprising of 125 patients across five tumor types. We filtered these data to just retain the macrophage populations using reference-based mapping \citep{seurat}. We then filtered low-quality cells (defined as having fewer than 1500 genes with non-zero counts) and run batch integration across all 125 patients using Harmony~\citep{harmony} in a curated subspace of 319 macrophage function genes. We then ran Leiden clustering~\citep{leiden} and characterized each population independently according to expression of important marker genes. Finally, we selected the \textit{gene36-gene46+} population as our target population, as it expressed known markers associated with suppressive macrophages and correlated with poor clinical response to common immunotherapies.

\begin{figure}[tb]
     \centering
     \includegraphics[width=0.75\textwidth]{figures/2022-03-14_condition-ranking.png}
    \caption{SystemMatch-produced ranking of \textit{in vitro} conditions based on their computed similarity to the target condition.}
    \label{fig:ranking}
\end{figure}

Fig.~\ref{fig:dataset}a shows the UMAP~\citep{umap} dimensionality reduction of the integrated macrophage atlas annotated by its marker genes. We identify inflammatory, suppressive, proliferating, and co-suppressive/inflammatory cell subsets, the last of which is not typically described in reviews of macrophage biology. For the purposes of this study, we target the suppressive population for \textit{in vitro} system optimization.

To study the heterogeneity of macrophage development systems, we differentiated monocytes \textit{in vitro} with either \textit{M-CSF} or \textit{GM-CSF} for 6 days and then stimulated these macrophages for an additional 4 days with one of nine combinations of cytokines typically understood to play a role in regulatory macrophage phenotype. Fig.~\ref{fig:dataset}b shows the UMAP dimensionality reduction of the \textit{in vitro} systems. We see that these systems primarily separate into three distinct groups, which can be described as pro-inflammatory, pro-suppressive, and basal.

Finally, to compare the \textit{in vitro} cells to our human tumor macrophage atlas, we show a heatmap of selected marker genes for both datasets (Fig.~\ref{fig:dataset}c). We observe qualitatively that not one of the \textit{in vitro} systems fully recapitulates the gene expression of the suppressive \textit{gene36-gene46+} population. This is further recapitulated by the similarity ranking of the \textit{in vitro} systems to the target population (Fig.~\ref{fig:ranking}), in which a) the best-ranked system using M-CSF, which serves as the basis for all systems most commonly used to generate suppressive macrophages \citep{mia2014optimized,fogg2020ovarian} is ranked 7th; and b) the closest system only reduces the distance to the target by 20\% from the least common system, which is well known to produce a pro-inflammatory response.

\begin{figure}
     \centering
     \begin{subfigure}[t]{0.3\textwidth}
         \vskip 0pt
         \centering
         \includegraphics[width=\textwidth]{figures/2022-03-13_groundtruth-umap-big-outline.png}
     \end{subfigure}
     \hfill
     \begin{subfigure}[t]{0.34\textwidth}
         \vskip 0pt
         \centering
         \includegraphics[width=\textwidth]{figures/2022-04-13_ranks_as_bar.png}
     \end{subfigure}
     \hfill
     \begin{subfigure}[t]{0.34\textwidth}
         \vskip 0pt
         \centering
         \includegraphics[width=\textwidth]{figures/2022-04-14_scores_by_genes_more_samples_200_iterations.png}
     \end{subfigure}
     \hfill
    \caption{Quantitative evaluation of distance metrics for ranking data with known ground truth. UMAP (left) shows the annotated held-out, ground truth macrophage data; computed distances (middle) between the held-out macrophage populations and the target population for each method; evaluation scores (right) of each distance metric relative to the known ground truth.}
    \label{fig:evaluation}
\end{figure}

\subsection{Evaluating similarity ranking methods}

To validate the effectiveness of SystemMatch, we compare held-out ground truth macrophage data from \citet{tang2018human} to our human tumor macrophage atlas. We preprocessed the held-out data by filtering out low-quality cells and annotated cells in the dataset using the same gene set and protocol as for the atlas, marking the macrophages as highly suppressive, moderately suppressive, or nonsuppressive. Fig.\ref{fig:evaluation} (left) shows the UMAP of the final annotated clusters of the ground truth macrophage dataset.

Because our target cells are suppressive macrophages, we expected the ranking of the ground truth cells regarding their proximity to the target suppressive cells to be from most suppressive (\textit{gene36-gene46+}) to least suppressive (\textit{gene36+gene46-}). The following distance metrics were considered (see Methods): 1) L2 distance on log-normalized pseudobulked gene expression; 2) L2 distance on z-scored, log-normalized pseudobulked gene expression; 3) Earth Mover's Distance on log-normalized single cell gene expression; 4) Earth Mover's Distance on z-scored, log-normalized single cell gene expression. On the full dataset, all four distance metrics correctly ranked the three datasets (Fig.~\ref{fig:evaluation} (middle)), so we define the following score to quantitatively evaluate a) how well it correctly ranked the held-out conditions and b) the separation of the most suppressive condition from the least suppressive condition.

For a given distance metric $m$, query subsets $\{{\mathbb{X}}_i ~|~ i \in \{\textit{gene36-gene46+}, \textit{gene10}, \textit{gene36+gene46-}\}$ and target population ${\mathbb{Y}}$, we compute distances between the conditions in the held-out dataset to the target condition giving a distance vector
$${d}^m = \left( m({\mathbb{X}}_{\textit{gene36-gene46+}}, {\mathbb{Y}}),  m({\mathbb{X}}_{\textit{gene10}}, {\mathbb{Y}}),  m({\mathbb{X}}_{\textit{gene36+gene46-}}, {\mathbb{Y}}) \right)$$
with corresponding rank vector ${r}^m$, where the expected rank vector ${r}^{
\textit{true}} = (1, 2, 3)$. For a ranking $r_m$ of length $n$, we define the score

$$\mathrm{Score(m)} = \left( \frac{\sum_i^n I({r}^m_i = r^{\textit{true}}_i)}{n} + \frac{{d}^m_{gene36+gene46-} - {d}^m_{gene36-gene46+}}{{d}^m_{gene36+gene46-}}\right) / 2$$

We compute this score for different levels of dataset corruption, which we implement by subsampling genes to remove information. We repeat this subsampling up to 200 times, or until the score converges. Fig.~\ref{fig:evaluation} (right) shows the scores for the four metrics; EMD with log-normalized genes performs best with the full gene set, but pseudo-bulked L2 distance is most robust to corruption.

\subsection{Recommending future experiments}

In order to leverage the recommendations of SystemMatch to further improve the \textit{in vitro} systems beyond those initially tested, we combine compositional perturbational autoencoders with a modified $k$-medoids algorithm to select a subset of possible combinations of the initially tested cytokines to test in the next experimental iterations. First, we use a modified form of CPA~\citep{lotfollahi2021compositional} to generate all possible double- and triple-combinations of cytokines (Fig.~\ref{fig:recommender} (left)). We see that \textit{M/GM-CSF} and \textit{CKD} drive the most significant differences in the predicted cells, which is consistent with our understanding of these cytokines. Further, we see as expected that the double-combinations of cytokines are generally more similar to the \textit{in vitro} single cytokine perturbations.

To validate the accuracy of our \textit{in silico} predictions, we compared the similarity of the held-out \textit{in vitro} triplet perturbations to the corresponding \textit{in silico} perturbations. Fig.~\ref{fig:recommender} (middle) shows the distance (using the L2 pseudobulk distance normalized to a range of $[0, 1]$) from each of the held-out triplets to the corresponding predicted triplets. All of the four are closest to the corresponding prediction. Equivalent heatmaps for the other distance metrics are shown in Fig.~\ref{fig:triplets}.

Next, we run our modified $k$-medoids algorithm on the average expression on each condition to select a subset of these double- and triple-combinations to test experimentally in future work to maximize the heterogeneity of all tested perturbations. Alternatively, we can use SystemMatch to test the proximity of the generated conditions; this will produce a less heterogeneous set of combinations for further testing but will optimize more directly towards the target condition. The result of this comparison is shown in Tab.~\ref{fig:ranking_in_silico}; interestingly, the top predicted condition uses M-CSF, even though none of the top six \textit{in vitro} conditions did. The \textit{in silico} preference (however mild) for \textit{M-CSF} over \textit{GM-CSF} is consistent with commonly used models of suppressive macrophages \citep{mia2014optimized,fogg2020ovarian}, giving confidence in our recommendations.

\begin{figure}[tb]
\begin{center}

     \centering
     \begin{subfigure}[b]{0.30\textwidth}
         \centering
         \includegraphics[width=\textwidth]{figures/2022-03-15_in-silico-umap-labelled.png}
     \end{subfigure}
     \hfill
     \begin{subfigure}[b]{0.36\textwidth}
         \centering
         \includegraphics[width=\textwidth]{figures/2022-04-13_l2_lg_heatmap_main.png}
     \end{subfigure}
     \hfill
     \begin{subfigure}[b]{0.25\textwidth}
         \centering
         \includegraphics[width=\textwidth]{figures/2022-04-20_experiment-recommender.png}
     \end{subfigure}
\end{center}
\caption{Sample output of the experiment recommender system. UMAP embedding of \textit{in vitro} and \textit{in silico} single cells (left), with \textit{in silico} predictions shown in grey; heatmap of normalized distances from predicted triplets to held-out triplets (middle); UMAP embedding of average expression per condition (right), with $k$-medoids--recommended perturbations highlighted in orange.}
\label{fig:recommender}
\end{figure}

\section{Discussion}

In this work we provide an end-to-end ML pipeline for assessing the fitness of and optimizing preclinical models to maximize their predictive power of clinical results. We do this through an iterative process of comparing model systems to the target \textit{in sapiens} population, applying perturbation prediction to suggest experimental changes, and iterating upon this process to produce preclinical models that are more representative of the relevant tumor or tissue data. We demonstrate our pipeline through a use case of developing an \textit{in vitro} system optimized to produce macrophages most similar to human tumor-derived suppressive macrophages, and we recommend a set of further experiments to fully explore the space of possible perturbations and ultimately generate a system more suited to developing drugs targeting these cells than those systems currently used in the literature.

We validate our pipeline's performance through a series of tests using ground truth data with a held-out \textit{in sapiens} public dataset in which we manually annotate cells known to be close and far from the target population. We show that a) our pipeline is largely robust to the choice of distance metric; b) distance metrics computed on average expression are generally most robust to data corruption; and c) in the limit of low data corruption, single-cell optimal transport metrics may outperform these "pseudobulk" methods. We then validate our perturbation prediction using experimentally generated \textit{in vitro} controls of cytokine combinations, through which we show that the \textit{in silico} generated combinations are most similar to the corresponding matching \textit{in vitro} populations.

We note that the robustness to the use of average expression (compared to single-cell) indicates that bulk RNA-seq could suffice in some cases for the query datasets; however, this would preclude the use of CPA, and the resulting lower-resolution data would provide less predictive power to generate \textit{in silico} samples. Additionally, in the case of heterogeneous model systems (e.g. in systems with developmental trajectories or large numbers of cycling cells), single-cell metrics may score more highly. We also note that the use of \textit{in silico} generated combinations of experimental conditions explores only a subset of possible model systems; systems generated with conditions not included in the initial experiment cannot be discovered in this way, and must instead be added to the iterative process by manual expert review.

In recent years, the utility of preclinical model data to evaluate the clinical relevance of a drug has been questioned. We show here, through a combination of single-cell genomic data and machine learning, a method by which these preclinical model systems can be optimized to enhance their potential predictive power. While ML methods have generated a great deal of impact in target discovery and validation, this is to our knowledge the first ML pipeline for actually establishing the preclinical model systems in which those targets can be evaluated.