Analytical method for calculation of deviations from intended dosages during multi-infusion

In this paper, a new method is presented that combines mechanical compliance effects with Poiseuille flow and push-out effects (“dead volume”) in one single mathematical framework for calculating dosing errors in multi-infusion set-ups. In contrast to existing numerical methods, our method produces explicit expressions that illustrate the mathematical dependencies of the dosing errors on hardware parameters and pump flow rate settings. Our new approach uses the Z-transform to model the contents of the catheter, and after implementation in Mathematica (Wolfram), explicit expressions are produced automatically. Consistency of the resulting analytical expressions has been examined for limiting cases, and three types of in-vitro measurements have been performed to obtain a first experimental test of the validity of the theoretical results. The relative contribution of various factors affecting the dosing errors, such as the Poiseuille flow profile, resistance and internal volume of the catheter, mechanical compliance of the syringes and the various pump flow rate settings, can now be discerned clearly in the structure of the expressions generated by our method. The in-vitro experiments showed a standard deviation between theory and experiment of 14% for the delay time in the catheter, and of 13% for the time duration of the dosing error bolus. Our method provides insight and predictability in a large range of possible situations involving many variables and dependencies, which is potentially very useful for e.g. the development of a fast, bed-side tool (“calculator”) that provides the clinician with a precise prediction of dosing errors and delay times interactively for many scenario’s. The interactive nature of such a device has now been made feasible by the fact that, using our method, explicit expressions are available for these situations, as opposed to conventional time-consuming numerical simulations.

in vivo studies [5][6][7][8][9], indicating that physical effects related to the infusion hardware are equally important to understand [10,11]. Especially in critical care, on the ICU and the OR, where multiple medications are typically delivered through one single lumen of a thin catheter, these physical effects may cause ambiguous and counter-intuitive discrepancies between the intended dose and the dose that has been actually delivered. This is particularly true if the dosing rate is adjusted ad hoc [12]; for example, fast-acting and critical inotropic drugs are often titrated, based on the mean arterial blood pressure (MABP).
There are three major factors that can produce significant deviations from the intended medication dosing rate scheme: (i) the length of the catheter causes a delay in the administration of the medication into the patient, and therefore, a mixture, corresponding to previous medication dosing rates, may still be present inside the catheter, i.e. the contents inside the catheter constitutes a "memory" in which the effects of previous medication dosing rates may be stored. This has been called "dead volume effect" in the literature [7,13]. (ii) The syringes that are used in clinical practice are far from ideal, because these syringes have a significant mechanical compliance primarily due to the compressibility of the rubber plunger inside the syringe [14]. Therefore, whenever the flow rate setting of one of the infusion pumps in a multi-infusion set-up is altered, pressure changes within the entire system cause a change in the deformation of compressible and expandable parts (i.e., the other syringes) within the system, and hence may cause a deviation from the intended flow rates. (iii) The low velocities of the fluids in catheters ensure that the flow inside these lines will be laminar (low Reynolds number), and hence exhibit a Poiseuille flow profile, in which the fluid particles near the central longitudinal axis of the catheter travel faster than fluid particles near the wall of the catheter, thus giving rise to a "mixing effect" of its own. Several infusion simulation studies have shown the importance of these problems by demonstrating the influence of the physical effects of dead volume, compliance and the Poiseuille effect on drug delivery in isolation: There are a number of studies that describe the "dead volume effect" in isolation from the syringe compressibility ("compliance") effect [7]. In simple cases, calculation of the "dead volume effect" is straightforward: If the actual volumetric flow rate u cath inside a single lumen catheter is assumed to be constant, and the internal volume V cath of the catheter (i.e., the internal volume as measured starting from the mixing point, at which the medications from all syringes come together, up to the distal catheter tip inside the vasculature of a patient) is known, then calculation of the "delay time" t delay , i.e. the time needed for a droplet of medication to travel through the dead volume before it reaches the blood stream of the patient, is simply t delay = V cath /u cath .
In clinical practice, however, the situation is often more complex: the flow rate may vary during the delay time due to changes in the pump flow rate settings, and the actual flow rate also typically differs from the pump flow rate setting value temporarily, due to the effects of mechanical compliance as mentioned above. These compliance effects have been studied in isolation as well [2,3,15]. Its basic mechanism is a "capacitor" effect, which has been quantified using the electric analog of a system containing capacitors, and the calculations of the dosing errors have been performed using the Laplace transform.
The third major factor, the mixing effect due to the Poiseuille profile [16], can be described as a convolution, as will be explained in more detail in "Appendix". It causes a spreading out of the dosing error in time, in which the first arrival of the dosing error occurs sooner then it would have without Poiseuille mixing effect.
As a result, due to the combination of the effects (i), (ii), and (iii) mentioned above, non-trivial deviations from the intended medication dosages scheme may occur, as will be explained below.
In an earlier paper, we have shown that the "dead volume effect" on one hand, and syringe compliance ("capacitor") effect on the other hand, produce opposite deviations from the pump flow rate settings in the actual drug output concentrations, making the net result hard to predict and often counterintuitive [4].
In this paper, however, we focus on the mathematical method of calculating dosing errors in multi-infusion setups in complex situations. The aim of the new method described in this paper is to obtain analytical expressions for the deviations from the intended dosages during multi-infusion that result from the combination of the effects (i) and (ii) described above. These analytical expressions contain dependencies on parameters that represent the physical variables in a multi-infusion set-up, such as: intended flow rates of the various pumps, compliances of the syringes, resistances of the various tubes, height differences within the multi-infusion set-up, etc.
As opposed to conventional numerical simulations, our objective therefore is to obtain explicit analytical expressions for these deviations, in which physical characteristics of the multi-infusion set-up are represented explicitly as variables.
We see the need for a tool for clinicians and medical physicist that provides understanding of the role that various physical parameters (i.e., characteristics of the multiinfusion hardware) play in the emergence of a dosing error. In order to make these roles explicit, we will use our new mathematical approach, presented in the "Method: analytical model, and in-vitro set-up" section of this paper, to explicate these roles in the form of direct mathematical relationships between physical hardware parameters and dosing errors.
For our mathematical model we used the strong symbolic calculation capabilities of the Mathematica package (Mathematica 10, Wolfram ® Inc., USA) to process Laplacetransforms and Z-transforms analytically. This will be explained in "Method: analytical model, and in-vitro set-up" section. Furthermore, we have performed in-vitro measurements in order to verify our mathematical results experimentally.

Method: analytical model, and in-vitro set-up
To start with, we analyse the total chain of causality, from a change in the pump flow rate setting value, up to the moment that a dosing error enters the bloodstream of the patient. As has been mentioned in the introduction, it is essential to recognize that the lumen of the catheter (i.e. the internal volume C inside the catheter, starting from the mixing point M, at which the medications from all syringes come together, up to the distal tip P of the infusion line inside the vasculature of a patient) constitutes a memory in which the effects of previous changes in pump flow rate settings are stored. In the explanation of our method, we focus on a single-lumen catheter. This however can be easily extended to a multi-lumen catheter. In order to represent this memory in the mathematical model, the internal volume of the catheter C between the points M and P is divided into tiny voxels A k , in which the index k runs from k = 0 at point P to k = N at point M (see Fig. 1), in which N is a very large number. Let L denote the length of the catheter. Hence, the length of a single voxel equals γ, with γ = L/N. In our mathematical model, we use a general set-up containing any number of infusion pumps. For simplicity, however, we start with three pumps, in which each of these pumps contains a solution of a different medication, in which we used colours ("Red", "Green", and "Blue", or R, G, and B, respectively) to denote the three different solutions. It is important to note that the R, G, and B denote the three different solutions as stored in their respective syringes, not the medications themselves. Furthermore, in this paper, the dosing errors will be expressed as volumes of the R, G, and B solutions, instead of dosages of the medications themselves. In Fig. 1a, the situation directly after t = 0 is depicted, i.e. the point in time at which a change in pump flow rate setting value of one of the pumps takes place. As a result, the first small voxel at k=N near the mixing point M is now being filled with a droplet ξ featuring a new mixing ratio between the solutions R, G, and B, resulting from the new pump flow rate setting values at t = 0. The rest of the voxels inside the line C, however, contain a fluid mixture that still has the old mixing ratio corresponding to the steady-state situation before the change in pump flow rate setting value. In Fig. 1b, the situation at t = t delay is depicted, i.e. the point in time at which this specific droplet ξ , that has been situated inside voxel A k=N at t = 0, has now reached the distal tip of the catheter at point P at t = t delay . Hence, at t = t delay , the entire line C has now been filled Example of a general multi-infusion set-up. In this example, the number of pumps is three. The method, however, allows for any number of pumps. a The catheter C contains a mixture of the fluids "R", "G", and "B", in which "R", "G", and "B" are the contents of syringe 1, syringe 2, and syringe 3, respectively. We consider this mixture to be completely homogeneous, despite the "stratified" rendering of the three colours in catheter C in this figure. b Situation at t = t delay ; since t delay is by definition the time needed to travel through the catheter C, the contents ξ that was inside the voxel k = N (near the mixing point M) at time t = 0 in a, has now reached the very tip P of the catheter at time t = t delay , and is entering the blood stream of the patient  16:18 with the new mixing ratio. The contents of voxel A k , for any value of k, is described by the 3D-vector ā k , in which in which The ratio a (R) k represents the volume fraction of fluid (R) inside voxel k, i.e. the volume of solution (R) inside voxel k divided by the total volume inside voxel k. The voxel for which k = 0 is the very last voxel at point P inside the catheter, i.e. the distal tip of the catheter. This specific voxel at k = 0 releases its contents ā k=0 directly into the blood stream of the patient. Now consider the situation in which the flow rate setting of the "Green" pump is changed suddenly at t = 0, whereas the flow rate settings of the "Red" and "Blue" pumps are never altered. If the "old" total flow rate inside the catheter C before t = 0 has been constant and equal to u old for all t < 0, and the "new" total flow rate is also constant and equals u new for all t > 0, then the partial flow rate u (R) patient of the "Red" fluid leaving the catheter at point P (and thus entering the patient) at t = 0 equals: in which the "old" mixing ratio refers to the mixing ratio produced by the pump rate settings of t < 0.
As is explained in Fig. 2, the presence of the "old" mixing ratio in Eq. (3) can give rise to effects that may be unexpected to physicians, in which a temporary increase in the partial flow rate of (R) entering the blood stream occurs, although the setting value of the flow rate of the pump of (R) has never been changed. We will refer to this unwanted temporary increase in the u (R) patient as the "push-out" effect, which has sometimes been called "dead volume" effect in the literature [7]. In Fig. 2a, the situation before t = 0 is depicted. In Fig. 2b, a temporary, unplanned, and undesirable increase in u (R) patient occurs (see figure legend for explanation). After t = t delay , the value of u (R) patient returns to its original value (Fig. 2c).
Furthermore, due to the mechanical compliance of the syringes, however (as mentioned at point (ii) in the introduction), the actual flow rate u cath (t) in the catheter does not follow the changes in pump flow rate setting values immediately. As a result, the actual flow rate u cath (t) is not exactly equal to the sum of the pump flow rate setting values u (R) pump (t), u (G) pump (t), and u (B) pump (t), due to the transient mechanical compliance effects that cause extra fluid to be stored in, or released from, the syringes directly after changes in pump flow rate setting values: = u new · fraction of (R) in the old mixing ratio for t > 0 and t < t delay Therefore, the general expression from which the "delay time" t delay needs to be established is: in which u cath depends on the compliance characteristics of the syringes as well.

New method for incorporating the memory effect of the catheter into an analytical model
In this paper, we will follow the following strategy to calculate dosing errors analytically (see Fig. 3): First, standard techniques [17] (involving the Laplace Transform) are used to calculate pressures and total flow rates (see Fig. 3a), in which the differences between the "Red", Fig. 2 Schematic, illustrating the "push-out" effect, i.e.: the temporary increase in the outflow of the "red" fluid into the patient due to an increase in speed of the "blue" pump, as is explained below. a Two infusion pumps, filled with a "blue" and a "red" solution, respectively, produce the same flow rate, which is constant in time.
The two flows are merged at the "mixing point" M, and subsequently travel towards point P, which represents the distal tip of the catheter where the fluid enters the blood stream of the patient. Since the flow rates of both pumps are equal, the catheter from point M to point P contains equal amounts of red and blue fluid, in the form of a mixture. This mixture is travelling with a flow rate that is twice the flow rate of each pump. A constant amount r of red fluid is entering the patient during each unit time interval (illustrated by a canister at point P for each unit time interval). b The same set-up, but now, at t = 0, the flow rate of the "blue" pump is suddenly increased by a factor 5. As a result, after t = 0, the amount of red fluid entering the patient at point P now temporarily equals 3r per unit time. Therefore, in b, the "memory effect" is visible: the distal part of the catheter (near P) still contains the "old" mixing ratio corresponding to the situation before t = 0. This "old" mixture is now being pushed out at the new speed; hence we dubbed this temporary increase in output of red fluid the "push-out" effect. In the literature, this effect is sometimes referred to as "dead volume" effect, or "catheter memory effect". c After t = t delay , the amount of red fluid entering the blood stream per unit time equals r again, which equals the intended dose, just like the situation before t = 0 in a "Blue" and "Green" solutions are disregarded; i.e. in this "color-blind" calculation, only the total flow rate inside the catheter u cath (t) is calculated. The mechanical compliance of the catheter is very small with respect to the mechanical compliance of the syringes [18], and therefore the mechanical compliance of the catheter is neglected in our model. As a result, the flow rate entering the tube at time t and the flow rate leaving the tube (entering the blood stream) at the same time t are both equal to u cath (t), disregarding the different constitutions that the fluids entering, and leaving, respectively, may have in terms of the partial fluids R, G, and B. Secondly, our new analytical method is introduced (see Fig. 3b), which enables incorporation of the "memory effect" of the catheter into the model. In order to make the analytical approach possible throughout the entire calculation up to the end of point (d) in the figure, a formulation in the Z-domain (using the Z-transform, which is a discrete variant of the Laplace transform) is introduced. This analytical method uses the results form the standard (Laplace-based) method [see (a)] as an input. This input has the form of general expressions for the various total flow rates as function of time. The Z-transform formulation in our method [see (b)] therefore does not replace the Laplace-based method from (a), but is used after it.
Third, the Z-domain formulation from (b) enables an easy, analytical, incorporation of the Poiseuille mixing effect (c) into the model, because the convolution-like nature of the Poiseuille mixing effect is reduced to a mere multiplication in the Z-domain. The mathematical details of the Poiseuille mixing effect are described in "Poiseuille mixing effect" section.
Finally (d), we derive expressions for the volume, and the first and second moment, of the dosing error as function of hardware parameters. These moments can be calculated directly within the Z-domain. From the first and second moment, key characteristics concerning timing and duration of the dosing error are derived. Throughout the entire process, we have used the symbolic calculation capabilities of the Mathematica package (Mathematica 10, Wolfram ® Inc., USA) to process Laplacetransforms and Z-transforms analytically, as can be seen in "Results" section.

Starting point: Laplace transform and Kirchhoff's laws
In this paper, we focus on the multi-infusion of aqueous solutions, in which all solutions have a viscosity close to the viscosity of water. Therefore, in our mathematical model, all fluids are assumed to have approximately the same viscosity. For the calculation of the pressures and total flow rates, we use an electric circuit model that is analogous to the infusion set-up, in which current sources (the infusion pumps), resistances (infusion lines and catheter), and capacitors (mechanical compliance of the syringes) are present, and in which Kirchhoff 's laws are applied on the voltages (pressures) and currents (flow rates) in the Laplace domain. The output U cath (s) of the calculation is the Laplace transform of the total flow rate u cath (t) inside the catheter. This u cath (t) equals the total flow rate entering the tube at time t as well as the total flow rate leaving the tube (entering the blood stream) at the same time t. Let Ŵ denote the complete set of these changes in pump flow rate setting values, together with the physical hardware parameters of the set-up, which are present in the form of explicit parameters in the analytical expressions for U cath (s). Using the Mathematica package (Mathematica 10, Wolfram ® Inc, USA), we were able to perform an inverse Laplace transform in order to retrieve an analytical expression for u cath (t, Ŵ) from the U cath (s, Ŵ), in which the above mentioned hardware parameters Ŵ are present as explicit variables. This result is rendered in Eq. (33) in "Results" section.

Key parameters of dosing errors
The dosing errors that are produced after changes in pump flow rate setting values in a multi-infusion set-up are of a temporary nature, i.e. "dead volume" and "mechanical compliance" (points (i) and (ii), respectively, as mentioned in the introduction) give rise to only a temporary deviation of the actual dose rates (entering the blood stream) from the intended ones (see e.g. Fig. 2). As a result, each of these dosing errors has the shape of a "bolus".
Let β (R) patient (t) denote such a bolus-shapes dosing error, i.e. let β (R) patient (t) denote the deviation from the intended partial flow rate [in this case of solution (R)], as a function of time, in which Defining t = 0 as the moment that a change in pump flow rate setting value takes place, and defining β (R) patient (t) as the dosing error due to (only) this specific change in pump flow rate setting value at t = 0, then we have by definition: in agreement with the bolus-shaped nature of the dosing error β (R) patient (t).
There are two basic mechanisms that contribute to the dosing errors, and the effect of each of these two mechanisms is restricted to a particular time interval: During the time interval 0 < t < t delay , the dosing error entering the blood stream is caused by the push-out effect, as explained in Fig. 2. During the time interval t > t delay , however, the dosing error entering the blood stream is caused by the effects of mechanical compliance, because the "push-out effect" is by definition restricted to the time interval 0 < t < t delay . The effects of mechanical compliance are particularly strong directly after changing the pump flow rate settings, and thus the effects of mechanical compliance (in the form of deviating mixtures) are entering the catheter at point M at t = 0, and hence entering the patient at point P directly after t = t delay .
In the following, we will use t POIS delay instead of t delay , because, as will be demonstrated below, the Poiseuille profile of the flow affects the t delay , resulting in a modified delay time that we will refer to as t POIS delay . Furthermore, we will consider a situation with two pumps ("Red" and "Green") only, and will change the flow rate setting of the "Green" pump only (at t = 0), leaving the flow rate setting of the "Red" pump (denoted as u (R) pump ) constant. This does not affect the applicability of our method in more general cases. The flow rate setting of the "Green" pump changes from u (G) delay is due to the push-out effect, and hence its calculation is relatively easy. The β (R) patient (t) for t > t POIS delay however is due to the mechanical compliance, and hence its calculation is much more complex. Therefore, we split the β (R) patient (t) into two parts according to the two corresponding distinct time intervals: in which τ = t − t POIS delay and ψ (R) patient (τ ) represents the effects of mechanical compliance for t > t POIS delay , and in which the dimensionless ψ (R) patient (τ ) is bolus-shaped of its own: In the following, we concentrate on finding analytical expressions for (the moments of ) ψ (R) patient (τ ). To this end, the discretized version ψ is defined as:

Contents of the catheter without poiseuille mixing effect
In this subsection, we derive a generic expression for the Z-Transform of the contents {a (R) k } of the catheter, without incorporating the Poiseuille mixing effect. The Poiseuille mixing effect will be incorporated in the model later on.
First, let tot (t) denote the distance that the fluid inside the catheter has traveled as a result of the total flow rate u cath , i.e.:  16:18 Furthermore, let u (R) M denote the actual partial flow rate of the Red fluid entering the catheter at the mixing point M.
For some voxel k inside the catheter, the value of the ratio a (R) k reflects the u (R) M (t * ) divided by the total flow rate u cath (t * ), at the time t * that the contents of that particular voxel entered the catheter at M, i.e.: in which γ = L/N. In order to facilitate a Z-transform, we rewrite Eq. (12) in the following form: The value of the Dirac delta δ(γ k − tot (t)) is zero everywhere, except around t = t * ; hence evaluation of the integral in Eq. (13) shows the equivalence of Eqs. (12) and (13): in which ε and η are both very close to zero.
Using the displacement rule from Z-transform theory [17], Eq. (13) now yields the following expression A(z) in the Z-domain: Here it needs to be mentioned that, since the Z-transform is merely a discrete form of the Laplace transform, it would be, theoretically, possible to replace this Z-transform by yet another Laplace transform. However, for clarity, we have chosen to use the Z-transform, because it connects to the "voxel-based" description of the contents inside the catheter in an intuitive way, and because it creates a clear distinction from the standard Laplace method used to calculate the pressures and flows that enter the mixing point M.
Similarly, we may define the basic input for the calculation of the dosing error of fluid (R) as the difference between u (R) M (i.e. the actual partial flow rate of the Red fluid entering the catheter at the mixing point M), and the intended u (R) pump (i.e. the intended partial flow rate of the Red fluid entering the catheter at the mixing point M):  16:18 This approach now enables easy incorporation of the Poiseuille mixing effect into our model using a simple multiplication in the Z-domain, as is explained below.

Poiseuille mixing effect
The basic input for calculation of the Poiseuille mixing effect is the input from the Red line at M in the form of which corresponds to A diff (z) in the z-domain [see Eq. (17)]. The mixing effect due to the Poiseuille profile causes a spreading out of the dosing error in time, in which the first arrival of the dosing error occurs sooner then it would without Poiseuille mixing effect. In a Poiseuile flow, the velocity along the centerline of the tube equals two times the average velocity. This implies that, after a change in pump flow rate settings at t = 0 , the tip of the parabolic flow profile reaches point P at the end of the catheter already at t = t delay /2. Hence, we define t POIS delay as: t POIS delay = t delay /2. Because the Poiseuille mixing effect is a phenomenon that is created within the catheter during the very passage of the liquid through the catheter from point M to point P, an accurate description of the full extent of this effect can only be given in the form of the contents of the very last voxel at point P, just before this contents is released into bloodstream of the patient. This is illustrated in Fig. 4, in which the top diagram (Fig. 4a) shows the Poiseuille profile along the length of the catheter from point M to point P, in which the tip of the red, innermost, parabola has just finished its journey through the catheter, and is just to be released into the blood stream. Each parabola in Fig. 4a constitutes a boundary between two volumes of liquid that once were subsequent, undistorted, "regular" voxels when they entered the catheter at M. These initially flat boundaries between voxels near M are being stretched and distorted into the parabola-shaped boundaries during their travel to P.
In "Appendix", it is derived that the Poiseuille flow profile causes a mixing effect that features a simple linear relation between concentration and distance along the line from M to P (see Fig. 4b), which is consistent with earlier work by e.g. Taylor [19] and Hutton and Thornberry [20].
Using this simple linear relationship, the constitution of the final voxel of the catheter at P can be defined in terms of the original, undistorted, voxels that once entered the catheter at M, as is illustrated in Fig. 4c. For instance, the contents of the "voxel * ", indicated by the symbol * , which is voxel nr i counting from point P, and which is limited by two downward sloping lines and the thick blue horizontal line segment in Fig. 4c, contributes to the final voxel at P with weight factor w i (indicated by the small vertical thick blue line segment). The dimensionless weight factor w i equals the voxel length γ times the first derivative of the expression x L+x , which is depicted in Fig. 4b, and which is established on the basis of geometrical reasoning using the triangles involved:  16:18 which yields Furthermore, the contents of each "voxel" i, with i running from 0 to infinity, equals that of voxel b (R) k−i , in which b (R) refers to the green curve in Fig. 4d, which is a "stretched-out" version of the original voxels a (R)diff j that entered the catheter at M, i.e.: which reflects the fact that the velocity along the centerline of the tube equals two times the average velocity. As a result, the summation of all contributions b (R) k−i with weight factors w i yields a (discrete) convolution, which is visible in the first term of the following equation, in which the fraction of the Red solution inside the last voxel at point P is (20)  The second term in Eq. (22) does not depend on the compliances and resistances in the multi-infusion set-up, but represents the remnant fluid from the old mixture due to the poiseuille flow profile.
When calculating the first and second moment of ψ (R)patient k (or t central and the variance σ 2 from Fig. 5, respectively), these moments are simply the sum of the moments of the two distinct terms in Eq. (22). The moments of the second term can be calculated directly in the spatial domain (using an upper limit to the amount of fluid). For the first term, however, the Z-transform is needed, as will be explained below. Therefore, we split the ψ (R)patient k according to:   in which B(z) is the Z-transform of b j , and W(z) is the Z-transform of w j , which yields: Since L/γ = N is a large number, this constitutes a rapidly converging series for z ≈ 1.

Calculation of moments of the dosing error distribution after T delay
We will now use the "theorems of moments" from Z-transform theory to provide expressions for key characteristics of the ψ (R) patient (τ ) that enters the patient. The zero'th moment yields: and hence the total volume of the dosing error ψ (R) patient (τ ) equals: The first moment yields: and hence Furthermore, using the second moment, the following general expression for σ is obtained, in which the typical "duration" or "width" of a dosing error is characterized as 2σ (see Fig. 5): In "Results" section, Eqs. (27)-(31) will be used to produce analytical expressions for various situations. Furthermore, also in "Results" section, Eqs. (27)-(31) will be evaluated for situations without the Poiseuille mixing effect, in which case the � (R)POIS (z) in these equations is replaced by merely A(z).

Spectrometric in-vitro set-up for feasibility tests
A schematic of the in-vitro set-up is rendered in Fig. 6. This allowed a spectrometer QE65000 (Ocean Optics ® , Dunedin, FL, USA), also connected to the same flowcell, to continuously measure an absorption spectrum of the dye mixture. Concentrations of each dye were acquired from this absorption spectrum [4]. After the flow cell, the cumulative flow rate was measured using three M12p flowmeters (Bronkhorst ® , Ruurlo, The Netherlands). Sample time of all the measurements was 1 s.

General result from laplace transform as a starting point
Performing the first step in the scheme rendered in Fig. 3, a general expression for u has been derived in the Laplace domain, using a known technique, i.e. Laplace Transform in combination with Kirchhoff laws, on the basis of the electric analog depicted in Fig. 7: , and u (G) downstep equals the change in pump flow rate setting of the Green pump at t = 0, and C 1 and C 2 are the mechanical compliances of the "green" and "red" syringe, respectively, and R 1 and R 2 the resistances of the "green" and "red" feeding lines, respectively. Meanwhile, the pump flow rate setting of the Red pump, u (R) pump remains unchanged all the time. Using the inverse Laplace transform function of Mathematica, the following general expression was obtained in the time domain: in which pump , u (G) old and u (G) downstep are rendered in Table 1.

Typical example: resulting dosing error during syringe exchange
We now consider the specific case in which the "Green" syringe is exchanged, during which the green pump is stopped and green line has been clamped (with a Kocher), but the red line is not clamped, and the settings of the red pump remain unchanged. Let T restart denote the time between clamping the green line and the reopening the green line, i.e. T restart denotes the time duration of the entire procedure of exchanging the green syringe. During the syringe exchange time interval, the red line still oozes red fluid in undiluted form into the first voxels of the catheter directly after the mixing point M. We now consider two scenario's within the general scenario of the syringe exchange with the green line being clamped off during the entire exchange time interval. In the first scenario, the red line is not clamped off, but the red pump is switched off (pump flow rate setting value is zero; still connected to the red line) simultaneously at the beginning of the syringe exchange procedure, and restarted to its original pump flow rate setting value at the very end of the syringe exchange procedure. In the second scenario, no actions are performed in relation to the red pump or red line at all. The first scenario entails that accumulation of undiluted red fluid directly beyond the mixing point is caused only by the compliance of the red syringe, in combination with the drop in pressure caused by the clamping of the green line. This is indicated by the text "compliance" under the second brace in Eq. 36. In the second scenario, however, the regular pumping action of the red pump (at flow rate u (R) pump ) adds to the accumulation of undiluted red fluid beyond the mixing point as well (indicated by the text "red pump on" under the first underbrace in Eq. 36). As a result, in the second scenario, the total dosing error Q (R) dosingerror (T exchange ) is the sum of the "compliance" and the "red pump on" effects. See Fig. 9. After t delay , this total dosing error is entering the blood stream of the patient within a very short time interval, due to the fact that the green pump has a high flow rate.
In Eq. (36), the second term ("compliance") has been calculated on the basis of integrating the right side of Eq. (33) over time.

Analytical results for reducing the flow rate of the green pump, without Poiseuille mixing effect
We now consider the more general case in which the pump flow rate setting of green pump, i.e. of the fast pump, is lowered with an amount u (G) downstep , without clamping any line. The volume Q (R) dosingerror does not depend on the specific distribution of the dosing error over time, and hence yields the same result as in Eq. (35).
Here we present the general result of the calculation of the first moment of A(z) without the Poiseuille flow effect:  16:18 Applying this result to the specific case described in Eq. (32), i.e. lowering the pump flow rate setting of the green pump, and substituting Eq. (33) into Eq. (37), yields the following result: In many clinical situations we have C 1 = C 2 = C and R 1 = R 2 = R, in which case Eq.
(38) reduces to: Applying Eq. (31) from the "Method: analytical model, and in-vitro set-up" section to the present situation, without the Poiseuille flow effect reads: Evaluation of Eq. (40) yields a very long and unwieldy expression. However, as has been noted before, in many clinical situations we have C 1 = C 2 = C and R 1 = R 2 = R, and, furthermore, for most catheters we have R ≪ R cath . Applying these assumptions, we obtain the following result: This result will be compared with the findings from the in-vitro experiments in "Results of in-vitro experiments, and comparison with theoretical predictions" section. If, in Eq. (41), the u (G) downstep would be very small with respect to u final , then this expression for σ approaches the familiar σ ≈ 2CR cath .

Analytical results for reducing the flow rate of the green pump, incorporating the Poiseuille mixing effect
Evaluation of Eq. (31) from "Method: analytical model, and in-vitro set-up" section to the situation of reducing the pump flow rate setting of the green pump (still without any clamping) once more, but now substituting � (R)POIS (z) [instead of A(z)] in order to incorporate the Poiseuille flow effect, yields a very long and unwieldy expression. Therefore, again, we apply the assumptions C 1 = C 2 = C and R 1 = R 2 = R and R ≪ R cath , which yields the following result: As can be seen in Eq. (42), the strength of the contribution of the Poiseuille mixing effect to the width σ POIS depends on the internal volume of the catheter V cath with respect to the time CR cath in combination with the flow rates (u (G) downstep ) and u final ; or, more precisely, on the ratio between 3(V cath ) 2 and . If the catheter has a large internal volume but the time CR cath is short and the flow rates are low, then the σ POIS in Eq. (42) reduces to: If, however, in Eq. (42), the u (G) downstep would be very small with respect to u final , and the time CR cath would be very large and the V cath would be small, then this expression for σ approaches σ POIS ≈ 4CR cath , which is two times larger than the σ ≈ 2CR cath that was calculated before when omitting the Poiseuille mixing effect. This factor two is a result from the fact that, in a Poiseuille flow, the velocity at the centerline of the catheter equals two times the average velocity.

Results of in-vitro experiments, and comparison with theoretical predictions
Three types of in-vitro experiments have been performed, in order to compare theoretical results with actual measurements. These three types of experiments are:

Discussion
In this paper, explicit expressions have been derived for the total volume (Q), central time point (t central ) and "width" or "duration" (2σ) of dosing errors as function of hardware parameters such as mechanical compliance of syringes, resistance of the catheter, length of the catheter, for two general cases of a change in pump flow rate setting value in a multi-infusion set-up, as well as for a typical case of syringe change-over. Consistency of the resulting analytical expressions has been examined for limiting cases, and, more importantly, three types of in-vitro measurements have been performed to obtain a first experimental test of the validity of the theoretical results derived in this paper.
The experimental measurement of the t delay (Fig. 11) showed a standard deviation between theory and experiment of 16.2 s, which is 14% of the total range of variations in t delay . The experiment measuring the t central (Fig. 12) showed a reasonable agreement between theory and measurements, in which standard deviation between theory and experiment was 8.4 s, which is 13% of the total range of variations in t central .

General findings from the expressions derived
In many of the resulting expressions as presented in "Results" section, the relative contribution of various factors affecting the dosing errors, such as the poiseuille mixing effect, can be discerned clearly in the structure of these expressions. The characteristic "time constant" R cath C is clearly recognisable in the expressions. Furthermore, the ratio between the u (G) downstep (the size of the change in pump flow rate setting value) and the u final (the stabilized final flow rate) appears as an important factor in determining the nature of the dosing error bolus, particulary how the dosing error will spread out in time (the "width" or "duration" (2σ) of the dosing error bolus). The contribution of the Poiseuille mixing effect is visible in the resulting expressions. Of particular clinical importance is the fact that the Poiseuille flow profile, once fully developed, causes a significant reduction of the time that is needed for a newly administered medication to reach the patient. This may come as an unexpected effect for the clinician, as may also the fact that the value of 2σ (the "spread out") is increased at the same time. The results in this paper may help to determine the magnitude of this "spread-out" effect as function of the hardware parameters (most notably, the characteristic R cath C time), and e.g. the length and resistance of the catheter.

Limitations of the method, and possible extensions
A number of limitations can be identified in our method in its present form; most apparent is the assumption that when a time interval of duration t POIS delay has lapsed after having changed a pump flow rate setting value, the actual flow rate has already stabilized and reached the value u final . This needs not to be true in a general case. However, using the same reasoning as presented in this paper, our method can be extended to include nonstable flow rates at t = t POIS delay as well. Another limitation may arise by the fact that in our method we did not examine other elements (other than just syringes, infusion lines, catheters, and pumps) that may be present in an infusion set-up, such as a non-return valve (preventing flow back from the catheter into a line towards a pump), or filters. The non-return valves may be incorporated into the method by restricting flow rates in the catheter to positive values or zero, whereas filters may be modeled on basis of a resistance and a mechanical compliance, as indicated in the literature (e.g. [21]). Another complicating factor may be the use of high viscosity fluids in the infusion set-up. In our method, as it stands now, all fluids in the system were assumed to be of the same viscosity. Incorporation of deviating (high) viscosities into our method would entail an extension of the laminar flow profile used in this paper, because differences in viscosity of mixing liquids may produce destabilization of laminar flow. An easy, but useful, extension of our method is the incorporation of the possibility that the height at which the syringes are positioned in a set-up near the patient, is altered during the period that a patient receives medication using the multi-infusion set-up. Since a change in height may be modeled as the addition of an extra pressure source within the model during the process of infusion, this results in an extra bolus of dosing error, and, hence, in the addition of some extra terms and factors in the equations in "Results" section. Finally, the parabolic flow profile of a laminar flow does not come into existence instantaneously at the mixing point; it has been calculated [20] how long, c.q. what distance, it takes for the flow profile to approach the parabolic shape. All of these complicating factors need to be incorporated into the model to make it more realistic. As far as we can anticipate now, we do not expect these extra factors to be incompatible with our general approach outlined in this paper; however, further research is needed.

Potential use of the results in clinical practice
Healthcare professionals working with infusion technology in critical care have expressed the desire for a real-time tool that visualizes the multi-infusion drug therapy, e.g. continuously calculates predictions to indicate when the drugs will be entering the blood stream of the patient and in what dose. Such a tool may also visualize the causal consequences of an intervention (i.e., change in a pump flow rate setting value), before a clinician decides to proceed with such an intervention. In order to develop such a tool in the future, a fast and generic model will be necessary, combining all the relevant physical effects. The desirability of such a visualization tool, and the mathematical modeling that is a pre-requisite for the development of such a tool in the future, is what prompted us to go beyond the state-of-the-art and to develop the fully analytical method described in this paper. We envision three types of developments in which the results from this paper may be useful: (i) an interactive tool, running synchronized with the multi-infusion system on a smartphone device or on a bed-side display, e.g., next to the vital signs display monitors. Such a tool could then be used in several cases, such as inotropic titration, or to visualize the effects of a syringe changeover, or even the consequence of changing the height of a pump during infusion. (ii) Another application of the method presented in this paper could be actual computer control of an infusion system. It has been shown that a computer-controlled pump with "knowledge" about the dead volume and the mixing effect within the dead volume can be useful in preventing overshoot [16]. Moreover, a control system with a feedback approach has also been attempted, where the mean arterial blood pressure was used to control the administration of a fast-acting vasodilator [22,23]. In both cases, however, increasingly complex situations and infusion setups were encountered where only the incorporation of all the interdependent physical effects, as described in this paper, would provide a sufficiently accurate prediction in order to make computer control feasible. (iii) The model can potentially also be used as an analysis and design tool, prior to investing in potential new infusion hardware. It is known that flow characteristics are influenced by valves [24,25], syringes [26], infusion lines and catheters [27][28][29], and filters [21]. By using the method from this paper, benefits of these components can be compared against potential tradeoffs.

Conclusions
We have developed a new method that combines mechanical compressibility (compliance) effects with poiseuille flow and push-out effects ("dead volume") in one single mathematical framework for calculating dosing errors in multi-infusion set-ups.
In contrast to existing numerical methods, our method produces explicit expressions that indicate the mathematical dependencies of the dosing errors on hardware parameters and pump flow rate settings.
The results from the in-vitro experiments show a reasonable to good agreement between measurements and theoretical results.
The relative contribution of various factors affecting the dosing errors, such as the poiseuille mixing effect, resistance and internal volume of the catheter, mechanical compliance of the syringes and the various pump flow rate settings, can now be discerned clearly in the structure of the expressions generated by our method.
This enables insight and predictability in a large range of possible situations involving many variables and dependencies, which is potentially very useful for e.g. the development of a fast, bed-side tool "calculator" that provides the clinician with a precise prediction of dosing errors and delay times interactively for many scenario's. The interactive nature of such a device has now been made feasible by the fact that, using our method, explicit expressions are available for these situations, as opposed to conventional timeconsuming numerical simulations. Other potential applications of our method involve analysis and design tools for new infusion hardware, and interactive devices connected to the infusion hardware that counteract impending dosing errors using predictive calculations.

Abbreviations
ICU: intensive care unit; OR: operation room; MABP: mean arterial blood pressure.

List of symbols
V cath internal volume of the catheter t delay time needed to travel through the catheter u cath volumetric flow rate inside the catheter M mixing point where fluid enters the catheter P tip of catheter where fluid enters the blood stream L length of the catheter N large dimensionless number, representing the number of voxels inside the catheter γ L/N, i.e. length of a single voxel inside the catheter t = 0 point in time at which a change in pump flow rate setting value takes place a (R) k fraction of the volume of voxel k that contains "red" solution u (R) patient (t) partial flow rate of the "red" solution that enters the blood stream at time t u (R) pump , u (G) pump pump flow rate setting of the "red", and "green", pump, respectively Ŵ complete set of flow rate settings and physical hardware parameters β (R) patient (t) dosing error of "red" solution, i.e. deviation from intended partial flow rate, as function of t t POIS delay = in which x is a function of r boundary and represents the distance that the boundary (see Fig. 13) has travelled in the longitudinal direction (along the length of the catheter), for a given radial position r boundary . The relation between x and r in Eq. (47) is depicted in Fig. 13: Conversely, solving r boundary from Eq. (47) for all values x ∈ [0, t u max ], we have: which is depicted in Fig. 14.
In order to calculate the amount of the fluid R that is present inside a thin, disk-shaped, cross-sectional volume S of the catheter at position x, expressed as a fraction f R (x, t) of the total fluid in the thin volume S, in which the disk-shaped volume S is perpendicular to the central longitudinal axis of the catheter at point x, we need to integrate over the cross-sectional area of the catheter at point x, from r = 0 to r = r boundary (x, t), and divide this by the total cross-sectional area πr 2 0 : (44) u avg = r 2 0 �p 8ηL (45) u max = 2u avg