1. Introduction

A successful drug discovery is an extremely time-consuming procedure, including initial target identification and validation, pre-clinical trials on animals, regulatory approval to start trials in humans, clinical trials, submission of marketing and manufacturing authorization, licensing review, product on sale, and post-marketing surveillance [11, 13, 14, 16, 21].

Figure 1 illustrates a flow chart of a typical drug discovery process. After target identification and validation, researchers in biology, biochemistry, and pharmaceutical science test diverse chemical compounds in order to find hits, which are compounds with desired effects on targets. With the development of robotics and high-speed computing technology, it is feasible to develop automatic systems that can screen a large number of biochemical compounds in a short period of time. Such an automatic compound screening and analyzing process is called high-throughput screening (HTS) in drug discovery of pharmaceutical industries. HTS provides a practical and efficient method to test a large number of synthetic compounds in miniaturized in vitro assays to identify hit targets of interest. Then, the chemical compounds that have therapeutic and useful pharmacological or biological activities, called leads, are evaluated and undergo lead optimization to identify promising lead compounds. Followed by the initial synthesis and animal testing in preclinical trials and three phases of clinical trials on humans, a drug can be put on market after the Food and Drug Administration (FDA) approval.

Figure 1. Drug discovery process flowchart

HTS systems in drug discovery run multiple and repetitive experiments and tests on various chemical and biological compounds. Some of these trials are repetitions for statistical analysis. Other trials have various small modifications in the process with the goal of finding the optimum composition of chemicals and biological agents. All these trials may easily cover hundreds of experiments and consume tremendous amount of time, even though most of the work is done by robotic systems. From a business perspective, it is essential to work quickly so that the results can be produced before the competition has a chance to catch up. There are several possible ways to schedule the experiments run more efficiently. Most importantly, all the small steps have to be performed in the right order and with the right timing. This is similar to cooking multiple dishes in a large kitchen. You have to consider availability of equipment as multiple process steps call for heating, washing, and cutting, and there are only so many pieces of equipment to perform such tasks. On one hand, you do not need to invest into purchasing another piece of equipment if you know that at times the equipment you already have is sitting idle because processes are not ready for it. At other times many processes are just sitting and waiting for their turns on that same piece of equipment.

The main objective of this paper is to provide a much more efficient, robust and accurate scheduling and control paradigm based on the just-in-time control criterion in order to minimize stock holding and costs, eliminate waste without further delays, and improve HTS system efficiency and throughput. This proposed scheduling method provides an algebraic description of the HTS processes and converts HTS processes into equations, which will generate the optimal scheduling of the testing prior to the computer simulations. This mathematical model is called “discrete-event system”, i.e. a discrete-state event-driven system of which the system evolution completely depends on the asynchronous occurrence of discrete events over time. The discrete-event system model allows us to model the HTS process as algebraic equations that can be easily entered into computers for finding the optimal scheduling to get things done in the fastest and most efficient way. Drug discovery processes inherently are not numerical quantities. Rather, they are so called discrete events. They still can be quantified with numbers, but the essential descriptions here are in a sequence of events and timing. Overall, the description must combine numerical values, such as time delays, sample sizes or weights, arrival times and so forth with abstract structures, which group them in records, branches, pages, as we have on the charts. However, computers prefer grouping numbers into matrices and vectors, because those algebraic structures can be handled much faster by today's computer hardware. In contrast, charts are handled by multiple if-than statements, which require programming in which most HTS lab users are not specialized. An intuitive explanation of this proposal is to perform process optimization by entering equations in a computer rather than tweaking the parameters of the scheduling programs on a computer, as it was done most of the time. Especially HTS lab users often are not computer scientists, so they would prefer to have an optimal solution in hand rather than programming the computer to find one feasible solution, which is usually not optimal. The successful completion of this project is expected to provide a better scheduling and control approach that decreases the HTS screening time, reduces the cost of automation, and ultimately improves the quality of HTS testing by improving the lead-finding ratio and reducing false positives and negatives. Furthermore, the scheduling and control methodology developed in this project can be potentially applied other industrial systems, VLSI manufacturing systems, and transportation networks.

2. Literature Review

The optimal scheduling of testing a large number of compounds requires innovative modeling, advanced scheduling algorithms, effective communication tools, and friendly user interfaces. However, the HTS scheduling is far more challenging than traditional scheduling in chemical or manufacturing engineering in the following aspects: (1) each batch may pass the same machine more than once; (2) more than one batch may be present in an HTS system at the same time; (3) there are no buffers between machines; (4) a batch may occupy multiple machines simultaneously; (5) there are upper time bounds, i.e. due dates for processes, or time window constraints because of process requirements; (6) operating schemes are strictly cyclic, meaning that the distance between two consecutive chemical batches is always constant. Due to these challenging factors, the optimal scheduling of HTS systems has been known as an “np-hard” problem and has not been thoroughly studied. The origin, evolution, current state of arts, challenges, and future opportunities for HTS systems can be found in [11, 13, 14, 16, 21]. From 1980 to 2011, the speed of HTS has increased from 10-100 compounds per week to 8000 tests per day. Drug discovery requires HTS processes have higher quality, higher cost-efficiency, and higher throughput in compound testing in order to reduce the number of false diagnostics, such as false positives and false negatives, and increase the number of compounds tested within a certain time. In this proposal, I investigate the highly automated HTS processes in drug discovery to identify the optimal scheduling for testing compounds. The expected outcome of this project will reduce compound waste, help produce consistent compound mixtures, and increase the number and quality of hits in the testing. A better scheduling algorithm and a better automation control paradigm proposed in this project target three key challenges in HTS systems: time, costs, and quality [13]. Time can be reduced by a more efficient, robust, and accurate scheduling algorithm resulting in less delay, fewer compound waste, and better assay adaptations. The optimal scheduling can reduce the cost of automation, personnel, consumables, and reagents in HTS systems. Moreover, the quality of the screening will be improved as a consequence of the improvement of the process efficiency.

Scheduling of testing almost an unlimited number of compounds in an HTS system belongs to a general class of cyclic scheduling problems [9]. Except for special cases, cyclic scheduling problems are np-hard, in other words, that the computational complexity of the scheduling problems is exponentially increasing with respect to the size of the systems. Previous research developed globally optimal solutions for scheduling in HTS systems using mixed integer linear programming (MIP) models [13]. There are several drawbacks of these MIP models. First, they are inefficient and require a substantial amount of time to compute an optimal solution for even a small size scheduling problem. This has limited their applications to only trivial problems and available MIP models cannot be used for real-time scheduling of HTS systems. Secondly, existing HTS scheduling methods do not employ any algorithms that take advantage the properties of HTS scheduling and improve computation efficiency. For example, available MIP models for HTS systems use optimization software packages (e.g., CPLEX) to find the optimal scheduling. Thirdly, none of the existing HTS scheduling methods has been tested on real-world HTS systems.

This research paper provides an efficient HTS scheduling model using discrete-event systems and validate it using real-world HTS systems. The scheduling in HTS systems is typically implemented by computer software embedded in the HTS robots [15]. A globally optimal scheduling method based on discrete-event system modeling was developed for cyclic HTS systems [12]. Discrete-event systems are event-driven systems where the state evolutions depend on the occurrences of asynchronous discrete events over time. Such cyclic scheduling is called overtaking [17] or interleaving [15]. Using the predetermined static schedule, an algebraic model of the HTS system's activities was established [4, 5]. Because static schedules do not perform well when deviations occur during run-time [15], a feedback controller is designed to generate an optimal just-in-time schedule at the presence of unexpected deviations from the cyclic scheme. To summarize, the existing online scheduling algorithms [15] do not perform well where malfunctions are present. The predetermined static schedules [4, 5, 12] require the schedules to be known in advance. This project integrates and adapts the algebraic modeling for HTS systems [4, 5] and advanced scheduling algorithms in manufacturing, and develops an innovative optimal scheduling and control strategy for HTS systems.

3. Max-plus Linear System Modeling

Discrete-event systems, which have been used to model communication networks, manufacturing systems, and transportation infrastructures, are used to mode HTS systems [12]. A special algebraic structure, Max-plus algebra [4, 5], is used for offline optimal scheduling and control synthesis. This proposal adopts the discrete-event system modeling approach because it provides a systematic method for modeling and controller synthesis for HTS systems. A very simple HTS system is shown Figure 2, where the robot deposits the chemical compounds in a microplate, which takes 5 seconds. Then the robot takes the microplate and drops off the microplate in the incubator, which takes another 2 seconds.

Figure 2. An example of transporting a microplate to the incubator for an HTS system (top) and its corresponding discrete-event system model (bottom)

This is a simple process that can be modeled by a discrete-event system model in Figure 2, in which the bars are called transitions and the circles are called places. The dark circles inside places are called tokens, which indicates the readiness of the system transitions. The numbers next to the places indicate the time duration needed for the process. The input variable u represents the time that the compound is ready. The state variable x₁ represents the time starting to drop off compound. After 5 time units, the state variable x₂ represents the microplate is ready to be picked up. The output variable y represents the finishing time after the robot drops off the microplate in the incubator, which takes 2 time units.

The discrete-event system model shown in Figure 2 is called timed-event graphs (TEGs) in the existing literature, in which each place (circle) has only one upstream and one downstream transition (bar). For instance, the process time of each transition in the TEG model shown in Figure 2 can be described by the following linear equations by comparing the time/event variables:

(1)

Moreover, the robot has to wait for 5 time units to finish dropping the compounds, so The output equation implies that the robot takes 2 time units to move the microplate to the incubator after the k-th compound was dropped (x₂(k)).

Typically, discrete-event systems are very complex, nonlinear, and hierarchical and often have synchronization behaviors in the occurrences of the events. The max-plus linear systems [1] are used to describe timed discrete-event systems by incorporating the traditional linear system theory in modeling and analysis of nonlinear synchronization behaviors in event-driven systems. This simple HTS example is used to illustrate the main advantage of max-plus linear system modeling for discrete-event systems, which is the adoption of the traditional linear system schematics (Eq. (1)), where denotes the max-operation, denotes the addition:

(2)

The states x, the controls u, the disturbances q, and the output y record the processing time in a manufacturing process. The main difference between the traditional linear difference equations and the max-plus linear systems is replacing the addition and multiplication operations by the max and addition operations, respectively. I was able to establish an elegant control system theory for max-plus linear systems in optimal scheduling [18], state feedback and output feedback controller design [19, 20], and state estimation theory for such systems [8, 10]. The benefit of discrete-event systems is the computational reduction independent of the number of transitions and places in the TEGs, which can be achieved using the toolbox MinMaxGD, a C++ library allowing to handle periodic series [7], interfaced with Scilab and MATLAB. Moreover, Eq.(2) can be transformed to the so-called γ-domain analysis, i.e. analogy of the z-transform in discrete-time linear system. Eq. (2) will become

(3)

where

4. Optimal Scheduling and Controller Design for HTS Systems using Max-Plus Linear Systems

If we have an enzyme compound screening, the microplates contain compounds, enzyme solution gets manually added during the barcode reader 1 min incubation, and then enzyme substrate is added by a multidrop dispenser, which is an automatic, programmable dispenser for micro volume dispensing. Then, this is followed by an incubation period or the reaction and then reading the plate on the envision, as a measure of the enzyme reaction. Time constraint would be most critical for the substrate addition, incubation and readout process to ensure that the readout is obtained in the linear range and in a consistent manner for all plates (9 minutes of incubation all the time). Time constraints may also be put in place to avoid the plate sitting without a lid for extended periods of time (so any evaporation is limited). The HTS system can be modelled as a special discrete-event systems, called a “timed-event graph” (TEG), which is a special discrete-event model with on one upstream and one downstream transition, respectively. Figure 3 illustrates the HTS workflow for the enzyme compound screening.

The corresponding variables are defined in Table 1. Figure 3 is a simplified time-event graph from a 35-transition state machine initially for easier understanding of the readers. The control signal is the time instant that the user starts to the k-th screening of the compound. Between robot picking up the microplate in the incubator, x₁ and robot dropping off the microplate in the incubator, x₂, the time duration of 221 seconds include the picking up the microplate, removing the lid, changing the grip position, move to the barcode reader, etc. Once the microplate is in the incubator already, it is required to stay in the incubator for 9 minutes, no more and no less ideally for consistency in order to avoid contamination and compound reactions. The k+1-th microplate can be picked up after the k-th microplate is in the incubator, which is represented by the transition from x₂ to x₁, with one place with one token. The k-th microplate cannot be picked up after 9 minute (540 seconds) of incubation before before the k-1-th microplate is back in the incubator, which is represented by the transition from x₄ to x₃, with one place with one token. The disturbances marked in red transitions represent all the possible malfunctions and delays during this process. The transition from x₃ to x₁ with -761 seconds are the time constraints to guarantee the incubation time is no longer than 9 minutes (540 seconds), i.e. x₁+221+540<x₃.

Figure 3. HTS workflow for the enzyme screening

Table 1. Variable Definitions

The TEG model has close connection with max-plus linear systems in Eq. (2), therefore, Figure 3 can be written in terms of the max-plus linear systems with the state space matrices given as

The system equation in Eq. (3) can be rewritten as follows by solving the implicit equation (3):

(4)

If we calculate the transfer functions Hyu and Hyq using the toolbox MinMaxGD [7] interfaced with Scilab and MATLAB. We can have:

The state and output series represent the time instants of the k-th screening. If the screening needs to meet certain standards with respect to the speed, cost, and throughput, for instant, denoted as the reference transfer function G_ref, see Figure 4. We can design an open-loop profiler controller such that G_ref=HP, where H can be H_yu or G_ref can be H_yq, if one wants to minimize the cycle time or maximize the throughput despites of disturbances, assuming some degrees of measurement with disturbances.

Figure 4. Model reference control in max-plus systems [19]

4.1. Disturbance Decoupling

The optimal solution for the prefilter can be calculated by the left and right residuation defined as follows, respectively:

Therefore, the optimal prefilter is obtained by the following argument for any disturbances:

Then, the optimal prefilter solving the model reference control is

Notice that P is not causal, we can obtain the causal projection by remov ing the negative time instants, i.e.

The causal prefilter can be realized through the time-event graph shown in Figure 5. The cyclic place with one token and 227 second time delay presents (227γ)^*, each element represent how to generate the control u from the disturbances. For instance, from q₁ to v, there is an empty place with no delay and no token, which is P_causal(1,1); from q₃ to v, there is a place with one token and 6 seconds time delay, which is P_causal(1,2); from q₂ to v, there is a place with four tokens and 6 second time delay, which is P_causal(1,3); and from q₄ to v, there is a place with five tokens and 147 seconds time delay, which is P_causal(1,4). The optimal scheduling generated by the residuation is given as below:

Where the disturbance Q(γ) is 5 second delay of the operation due to malfunction on x₁. If the system runs autonomously without disturbances, the scheduling is the following:

You may observe for every transition, it is exactly 5 second delay of the operation. Therefore, the prefilter P_causal provides the optimal scheduling such that the system will start running right away after the malfunction is cleared without further delay based on the just-in-time criterion.

Figure 5. Timed-event graph model for the HTS system with prefilter

5. Conclusions

Discrete-event systems are event-driven systems where the state evolutions depend on the occurrences of asynchronous discrete events over time. Max-plus linear systems are used to model a special discrete event systems with up to one upward transition and up to one downward transition. The advantage of the max-plus linear systems is modelling the synchronization and nonlinearity of discrete event systems in a linear manner, especially the optimal scheduling problem in discrete event systems. This paper used the max-plus linear system modelling to find a prefilter solution to the model reference control in high throughput screening in drug discovery. Such a prefilter will generate an optimal scheduling despites of the system malfunctions and delays. The main results are illustrated using enzyme screening compound screening applications. In fact, state feedback and output feedback can also be designed for such system to achieve optimal scheduling. Moreover, industrial standard scheduling algorithms can be used to high-throughput screening systems as well.