HUMAN SKIN BURNS INDUCED BY DEFIBRILLATOR DEFAULT CURRENT: A MATHEMATICAL SIMULATION MODEL JAMESH. K. Yoo Deparfmem of Bioengineering, The University of Texas Health Center at San Antonio, 7703 Floyd Curl Drive, San Antonio, Texas 78284 (USA)
THOMASJ. WHITE Deparrment of Mathematics, San Antonio College, San Antonio, Texas (USA)
and DAVIDL. STONER Biomedical Service, Wilford Hall USAF Medical Center, San Antonio, Texas (USA)
(Received: 18 May, 1976)
SUMMARY
A preliminary model was preparedfor the development of superficial burns which can result from the interaction of jbulty dejibrillator and other electrical operating room equipment. A thermal model previously used to predict retina1 coagulation was adapted to give temperature rise history of the skin. Skin resistaizce is assumed to decrease rapidly due to the instability of a small volume near the surface. Heating due to dissipation of the dejibrillator capacitor energy in this skin volume is shown to be high enough to cause burns. A sensitivity analysis is included to isolate the dominant input parameters. Recommendations are made for future research.
SOMMAIRE
Un modgle prkliminaire simule les brtilures superficielles quipeuvent Otre causiespar une d$brillation dkfectueuse ou d’autres manipulations 6lectriques. Un modkle thermique, pr&ctdemment utilise pour predire la coagulation rktinienne, a 6t6 adapt& pour retracer l’histoire de l’tY&vationde temperature de la peau. La rksistance de la peau est supposie dkcroitre rapidement ir cause de l’instabilit6 d’un petit volume situ& prPs de la surface. L’&hau#ement dii ci la dissipation d’tnergie d’un dbjibrillateur au niveau de la peau est sufisamment important pour provoquer des brtilures. Une analyse fine est incluse dans ce travail, ajin de mettre en evidence les principaux parametres d’entrte. Des recommandations sont faites pour des recherches futures. 109 Int. J. Bio-Medical Computing (8) (1977)~-10 Applied Science Publishers Ltd, England, 1977
Printed in Great Britain
110
JAMES H. K. YOO, THOMAS J. WHITE, DAVID L. STONER
PARTIAL LIST OF SYMBOLS
AJ(J): AP: BF,.: BF, : C: C(J):
DAMAGE: ~FR: G: H: K: L: POT: PRM: PRS: ;. RB: RE(I): RES( I, J): RED: RIM: R,: RMIN: RMAX: RN: RP: RR: RS: RSKIN: S(r, z, 2): SR:
Annular ring area (cm2) Defined as the area for RP or rc (2.52 - RIM2) (cm2) Blood flow rate in the r direction (cmjsec) Blood flow rate in the z direction (cmjsec) Circuit capacitance (PFarad) Current through the Jth cylindrical annulus Fractional skin damage Overall approximate blood change rate (set- ‘) The total conductance of the instability volume (ohm-‘) Coefficient of surface heat transfer (watt cmm2 “C- ‘) Thermal conductivity (watt cm- r “C- ‘) Circuit inductance (henry) The applied defibrillator potential of the form, I’, exp (- rt) sin (WC) (volt) A . RMAX( ‘Cm’) used in reaction rate for eduction of specific resistance Voltage drop across total skin resistance (volt) Power expended in the instability volume (watt) Initial charge on C (Coulomb) Internal circuit resistance of the defibrillator pulse generator (ohm) The lumped body resistance (ohm) Resistance of the Zth annular ring segment around the skinbreakdown instability centre (ohm) Resistance of elemental layer I of annular ring J of the instability volume (ohm) The initial reduction factor applied to the specific resistances in the instability volume Instability radius (cm) Load resistance (ohm) The assumed minimum thermo-chemical reaction rate (set-‘) The assumed maximum thermo-chemical reaction rate (set- ‘) The radius ofthe r boundary for the modelvolume. It is assumed to be 2.5 cm The remaining skin resistance to current flow from the electrode, outside the instability area. Assumed to be 170 ohm in the model The reaction rate for reduction of specific resistance (set-‘) Resistance of the instability volume (ohm) Lumped resistance of skin to current flow (ohm) Source strength (heat generation source) confined to the instability volume (watt cme3) Specific resistance (ohm-cm)
HUMAN
T: TMAX: t: ZM:
AZ: PC:
Ar:
SKIN BURNS
111
minus the original skin Temperature rise, i.e. temperature temperature (“C) Maximum temperature rise (“C) Time elapsed since the voltage application (set) Half the total depth of the model volume, assumed to be 2.5cm Grid spacing in the uniform part of the axial grid (cm) Volumetric specific heat (joule cm- 3 ‘C- ‘) Grid spacing in the uniform part of the radial grid (cm)
INTRODUCTION
The expanding use of complex electronic instrumentation in the Operating Room and Intensive Care Unit is increasing the risk to the patient of severe burns and electrocution. Previous papers (Yoo and Stoner, 1975; Stoner et al., 1976) described skin burns resulting from the interaction of a faulty defibrillator and an electrosurgical unit with a grounded dispersive plate. In this case, a short developed between the defibrillator switching circuitry and the discharge circuit. The short grounded one side of the discharge capacitor, thereby allowing the defibrillator to discharge through any grounded points on the patient. Thus, the grounded electrosurgical dispersive plate provided a second discharge path for the defibrillator. The burns resulting from this second discharge path were 3-10 mm in diameter and occurred on the medial and posterior aspects of the right popliteal space and the upper two-thirds of the lower leg, as shown in Fig. 1 (Stoner et al., 1976). The number of burns correlated with the number of defibrillation attempts and these burns were caused by current arcing from the patient’s leg to the grounded electrosurgical dispersive plate. In an effort to more fully understand the burn process, a preliminary model for the development of superficial burns has been devised. Insight into the physiology of the burn process can now be gained in two ways. First, by describing physiological changes in skin histology with physical changes in engineering components, one can develop a better understanding of the effect of surgical procedures on the physiological properties of the skin. Secondly, modelling allows one to simulate a multitude of burn conditions as functions of a variety of parameters which could not realistically be duplicated experimentally. Adjustment of specific resistances in the skin volume corresponding to the observed burns is used to achieve predicted temperature rises of 1l-48 “C within the first 5 to 18 msec of the damped oscillatory defibrillator pulse. While damage may occur because of the passage of electrical current through the skin layers due to mechanisms other than heating, this range of temperature increments corresponds to those known to cause burns in skin and retina following absorption of incident radiation. While it is clear that skin burns are not a primary concern during open heart
112
JAMES H. K. YOO, THOMAS J. WHITE, DAVID L. STONER
Fig. I.
Photograph
of patient’s
burns immediately
after the incident.
surgery, reduction of additional patient insult is a valid concern. The thermalelectrical model described here could lead to a better understanding of the mechanism of electrical burns when extended to include a description of the mechanism for the initial reduction in skin resistance necessary to allow for high current flow through the outer layers. A detailed sensitivity analysis covering both skin and defibrillator parameters is included. This analysis shows that the behaviour of skin resistance with temperature rise and the initial reduction due to breakdown at the electrode-skin interface and in the horny and basal skin layers are controlling factors for adequate damage prediction.
PHYSICAL, THERMAL AND ELECTRICAL
MODEL
A typical electrode-skin interface equivalent circuit is given in Fig. 2 (Strong, 1970). Typical skin impedance value for a 1 cm2 plate is 10 kR + jo (10 pF), where o is the angular frequency (Strong, 1970). In this paper, however, the capacitive loading is
HUMAN SKIN BURNS
113
ignored assuming that the initial punch-through current surge across the skin-electrode interface diminishes the capacitive effect to negligible magnitude at the specific instability point. No interface contact impedance is included in this model. ELECTRODE INTERFACE
Fig. 2.
IMPEDANCE
Simplified typical electrode-skin interface equivalent circuit.
This simplification is also motivated by calculational convenience during the mathematical model analysis. Research on the electrical properties of skin is not sufficient to allow adequate modelling of human skin impedance. In addition, enormous individualvariations of human skin properties and unknown thermal and electrical characteristics of the skin-electrode contact dictate this simplification. Therefore, in this paper, only gross estimates for the skin properties could be used. Since no models of skin damage due to passage of electrical current seem to be available, skin properties used are those which have been measured in connection with studies of damage due to exposure to incident radiation in the ultraviolet, visible and infrared parts of the spectrum. The work of other authors may be consulted for extensive literature on this type of injury (Stolwijk and Hardy, 1965; Hardy et al., 1965; Lipkin and Hardy, 1954; Weaver and Stoll, 1969). The burns described in this paper are confined to the first few millimetres of skin depth so that the parameters used with radiant energy injury models should be accurate below this depth. The electrical and thermal properties must change radically with the development of the initial instability and consequent very high current flow. No model is as yet available for this breakdown in skin resistance andwe have therefore simulated this process by an initial reduction in specific resistance and have thereafter continued the reduction of specific resistance according to an assumed rate of change dependent on temperature rise. The main purpose of the analysis of predicted temporal and spatial temperature profiles in this paper is thus limited to (1) illustration of the development of skin
JAMES H. K. YOO, THOMAS J. WHITE, DAVID L. STONER
114
burns due to absorption of electrical energy in the skin resistance and (2) exploration of better ways to use electrical surgical aids to minimise patient injury. Table 1 summarises the physical, thermal and electrical properties of the skin-electrode model. Layer depths are taken from the work of Stolwijk and Hardy (1965) and Hardy et al. (1965). Skin thermal parameters, volumetric specific heat, pc, and conductivity, K, are taken from the work of Stolwijk and Hardy (1965) Hardy et al. (1965), Lipkin and Hardy (1954), Weaver and Stoll(1969) and Draper and Boag (197 1). Specific resistances, SR, are estimates based on the work of Geddes and Baker (1967, 1968). TABLE 1 SKIN-ELECTRODE THERMAL AND Layer
ELECTRICAL
Depth
K x lo3
(watt cm-
(cm) Bedding Electrode Horny Basal Vascular Fatty Muscle
PARAMETERS
2.07 0.41 0.01 0.01 0.07 1.11 1.32
0.600 1.59 4.35 4.40 4-31 2.12 4.83
‘C-’ )
5.00 2.50 2.09 3.98 4.19 2.10 4.10
SR
(ohm-cm)
10: 1:: 7000 1500
The electrode is modelled as a thin aluminium sheet ( w 50 cc)glued to a wood backing about 4mm thick. Detailed consideration of the thermal resistances at the skin-electrode boundary would include those due to skin-conductive grease, greasealuminium foil, aluminium foil-glue, and glue-wood backing. All of these have been combined into one thermal resistance estimated as 1 cm2 “C watt- ‘, There are no thermal resistances between internal layer surfaces. The electrode is supported by 2 cm thick bedding and there is 10 cm2 “C- ’watt- 1 thermal resistance between the bedding and the electrode wood backing. Thermal parameters for the electrode and bedding are also given in Table 1 and are based in part on the work of Lipkin and Hardy (1954). Figure 3 shows the physical arrangement of the skin model layers and Fig. 4 the
electrode bedding
I
backing I horny Fig. 3.
vascular layer
I basal
fatty layer
Physical model of the skin layers.
I
muscle
115
HUMAN SKIN BURNS
corresponding electrical resistance model. Simulation of the initial breakdown m skin resistance is confined to a cylinder of 3 mm depth and radius RIM. RIM is varied in the computer runs from 0.66 mm to 3.5 mm. In this volume, referred to hereafter as the instability volume, temperature rise, T,willvary with both depth and radius as a function of time. Therefore, in order to more accurately simulate further reduction of skin layer specific resistances in thisvolume due to temperature rise, the
Fig. 4.
Electrical
resistive model.
model cylinder is divided into 10 annular cylinders, the innermost with radius Ar/2 and the remaining 9 with radius Ar so that Ar = RIMi9.5. Each cylindrical annulus is further divided into 16 layers: 2 each in the horny and basal skin layers, 9 in the vascular and 3 in the fatty skin layers. Each of these 160 elemental volumes comprising the instability volume has resistance RES(Z, J) computed as:
RES(Z,J) = SR(Z,J)AZ(I)/AJ(J)
(1)
where SR(Z,J)is specific resistance at depth Z(Z), AZ(I) is resistance layer depth, AJ(J) is annular ring area and I and J are axial and radial computer program indices, respectively. The resistance of annular ring J is then: 16
RE(J) = (1 I=1
SR(I, J> AZ(I)
Y
NJ),
J=
1,2,...,10
(2)
The total conductance of the instability volume is: 10
G= c
[l/W41
(31
J=l
In Fig. 4 RP represents the remaining skin resistance to current flow from the electrode, assumed to be 5 cm in diameter. RP is thus the resistance of a cylinder of inner radius RIM and outer radius 2.5 cm and is in parallel with the resistance of the instability volume. Specific resistance in RP is assumed not to be affected by temperature rise. RP may be calculated from:
116
JAMES
H. K. YOO, THOMAS
J. WHITE,
DAVID
L. STONER
RP = (O*OlSRl + O.OlSR2 + 0*07SR3 + 0.21SR4)jAP
(4)
where SRI, SR2, SR3 and SR4 are initial specific resistances for the horny, basal, vascular and fatty layers, respectively. AP = 7r(2*.5’- RIM*) is the area for RP. For the model used here RP is about 170 ohm. Since electrical resistance between the electrode and the horny layer is zero, the remaining resistance in the electrical path is the lumped body resistance RB (Fig. 4). RB is modelled as a cylinder 5 cm in radius, 50 cm long, with specific resistance of 250 ohm-cm. For these values: RB = 1590hm RB is not affected by temperature
rise in the skin. The applied defibrillator potential, POT in Fig. 4, has the form: POT = V,, exp ( - rt) sin (or) (volt)
(5)
This form is a lumped constant form of the voltage appearing across the load in a series RCL circuit. The three constants have been simplified as: r = (R + RJ2L ~~{&(!s~)‘-’ V, = QR,/LCw
where R, is load resistance; R is internal circuit resistance; L is circuit inductance; C is circuit capacitance; Q is initial charge on C, and w is real for underdamping. A reasonable form of reaction rate for reduction of specific resistance with increasing temperature is suggested by Milhorn (1966) and Masterson (1975) as the solution of: dRR = A. RR(RMAX
dT
where RR(T) is rate with maximum RMAX, The solution, which includes T = 0, is:
- RR)
A is a constant and T temperature rise.
RMAX
RR=
l+
RMAX
- RMZN
RMZN
(7) exp(-PRM.
T)
where RMZN, the minimum rate, occurs at T = 0, RMAX is the asymptotic rate at T = co and PRM = A. RMAX. Each of the resistances RES(Z, J) in the unstable volume is reduced with increasing temperature at this rate by integration of: dSR. dt
-=
-RR.SR
(8)
117
HUMANSKINBURNS
or: (9)
SR = SROexp( - ~ORRP-(*)ldl)
where SR, is initial reduced specific resistance and T(t) is temperature rise as a function of time for the elemental resistance. COMPUTATIONALDETAILSOFTHEMATHEMATICALMODEL
Figure 5 shows the physical model used in the machine computations. The equations to be solved for temperature rise, T, are:
(10) T(r, z, 0) = 0
(11)
R(RN, z, t) = T(r, 0, t) = T(r, 2ZM, t) = 0
(12)
where RN (2.5cm) is the radius of the r boundary, ZM (2.5cm) is half the total depth, and S(r, z, t) is source strength in watt cm -3 due to current flow confined to the instability volume. The method of solution is described in Mainster et al. (1970). instability
electrode
Fig. 5.
Schematic
three-dimensional
computational
model
The total volume used for heat conduction is about 100cm3 and is considered adequate for simulation of an infinite medium when thevolume of primary interest is less than lo-* cm3. The grid used has 3 1 points in the r direction and 37 points in the z direction with 20 uniform grid spaces in r and 16 uniform axial grid spaces centred at the middle of the z axis. The first point in the horny layer has index 17. AZin the uniform portion of the grid is 50~. Ar in the uniform r grid is RIMl9.5 when the instability contains 10 grid points in the radial direction. Grid stretching factors are 1.7 for the z grid and 1.55 for r when RIM = 1 mm. (Varies with RIM). (a) Thermal resistance and discontinuous conductivity Treatment of discontinuities in conductivity between layers and layer boundary
118
JAMES H. K. YOO, THOMAS J. WHITE, DAVID L. STONER
thermal resistance is achieved by modification of the approximating difference equation coefficients for z (Mainster et al., 1970). Figure 6 shows temperature rise at four grid points near a surface of separation. Heat flux must be continuous in the z direction at the boundary: +I,
= --Kg1
(13)
where n indicates normal derivative and indices 1 and 2 refer to opposite sides of the surface. Ti_l, Ti, Ti + 1 and Ti+2 are temperature rises at grid points zi_ I, zi, zi+ 1 and Zi+z, respectively. iij;and iii+, are temperature rises on either side of the surface
i i+l i+2 Kit-+-, Ki+t
i-l Fig. 6.
Media boundary with Newton cooling.
of separation of the two media of different conductivities, Ki and Ki+ 1. We assume Newton cooling at the surface with thermal resistance l/H (cm2 “C watt- ‘). Then: and :
-&+1
‘ii+1- K+1 (zi+l
_
zi)
=
H(‘ii+1- T>=A+1
(15)
2 Fij- Ti
= H(‘ii-K(zi + 1 - zi> 2
T+,> =fi
(14)
These two equations can be solved for the normal heat fluxes with the result that:
-fi
=A+1 =
KiKi+,(Ti+l - ri)
(Ki+Ki+l)(i’+;-z’)+5+
(16)
HUMANSKINBURNS
119
In the difference equations use of this result at z index i gives: 2Ki zi
+
I
-
zi
-
1
Ki+,)(Zi+l
-Zi)
2
Ki+I
while at z index i+ 1:
Ti - Ti_,
Ti+, - Ti +
+&KC+,
-
Zi-Zi-1
H
I -(Ki+
Kit,)(Zit,
-Zi>
Ki
These difference equations (eqns. (17) and (18)) are used at the bedding-electrode backing boundary and at the electrode-horny layer boundary. When there is perfect thermal contact as between skin layers, l/H = 0, T = T+ 1 and the term with H drops out of the equations. Where the grid is uniform and there is no thermal resistance ((l/H) = 0) and no conductivity difference (Ki = Kit 1), eqns. (17) and (18) reduce to: __ (d)Z2
,+,-2Ti+Ti_,} iT,
(19)
and :
as expected for centred second differences.
(b) Bloodflow in vascular layer Although the defibrillator pulse is short, an estimate of the effect of blood flow in the vascular layer has been included in the interest of possible further applications of the program. Since the co-ordinate system is cylindrical, simulation of blood flow cooling is restricted to axial and radial flow. The fluxes in the r and z directions when allowance is made for fluid motion are: f,=
-K&cT.BF,
(21)
and :
f,= -Kg
+ pcT. BF,
(22)
where BF,. and BF, are flow rates in the r and z directions, respectively, in cm/set. The terms involving derivatives with respect to r in eqn. (10) become:
120
JAMES
-
H. K. YOO, THOMAS
fi(rf,)=!A&rKg (
d!?+”
ar
- ii(rpcT.
J. WHITE,
DAVID
L. STONER
BF,)
>
(K!? ar) (“F “irF)TP+)?$ -pc
___C+__L
(2%
and :
BF, has been assumed to be zero and flow in the radial direction has been used to approximate blood flow parallel to the skin layers. BF, is assumed to go to zero like r as r + 0 so that: BF,+, dBF = constant r ar In view of the overall approximation, BF,/r is assumed to be this same constant and eqn. (10) is modified to: K --pcr.BFR r
$+
Kg+:
-pcBFR.T=pcg-S(r,z,l)
(25)
where BFR in units of set-’ is to be interpreted as an overall approximate blood change rate. In all computations, BFR = 5 set- ‘. Partial difference coefficients are modified from those reported by Mainster et al. (1970) to include the additional terms in eqn. (25).
S(r, z, t) The voltage drop across the total skin resistance is:
(c) Source strength
PRS = RP. POT] [RP + RB(l
+ G. RP)]
(26)
where G is given in eqn. (3), RP in eqn. (4) and POTin eqn. (5). The current through the elemental resistances with index J is: C(J) = PRSIRE(J)
(27)
and : S(Z, J, t> = RES(Z, 4[C(412/
[Az(Z)AJ(J)1
cw
is the source strength in watt cm -3 for elemental resistance RES(Z, J). Source strengthvaries with time through thevariation of POT, RES(Z, .I), RE(J) and C(J). Equations (7) and (9) are used with averagevalue of temperature rise at grid point I, J during time step k to reduce the specific resistance in eqn. (1). The RES(Z, 4 from eqn. (1) then changes the source strength (eqn. (28)) for use in time step k + 1
HUMAN SKIN BURNS
121
so that the source term closely follows the temperature history in the instability volume according to the parameters chosen for the reaction rate of eqn. (7). (d) Initial instability Since, the appearance of the burn is associated with random locations on the electrode and since the mechanism of initial breakdown in the resistance of the outer skin layers is not understood, the initial instability has been simulated by reduction of all specific resistances in the assumed instability volume prior to the start of calculation of temperature rise by the factor RED. As indicated above, the instability volume is made to correspond to the observed size and shape of the defibrillator leakage current burns, a few millimetres in radius and depth. The breakdown in resistance used to estimate the damaging current flow is confined to this volume. (e) Damage prediction
A skin damage prediction calculation is made a part of the computer program using a rate also defined by eqn. (7). Skin fractional damage is defined by DAMACE=l.-exp(-[/R[T(t)]dt)
(29)
Parameter values used are RMAX = 1000 set- ‘, RMZN = 0.01 set- ‘, and A = 0.0009 set “C - ‘. (f) Digital computing machine requirements
All runs were made on the San Antonio College IBM 370-158 computer by Mr Robert Lay of the Data Processing Department Staff. In Fortran H, storage required is 228K bytes for compilation. CPU running time averages about 20 set for 34 time steps or about l/2 msec/time step/grid point for the 1147 (37 x 3 1) point grid used.
RESULTS
Table 2 summarises the results of use of the program and model described in the previous sections. Run No. 1 was considered as the standard and succeeding runs of Table 2 show only those parameters which were changed. ‘V,, r, and o refer to defibrillator potential as shown in eqn. (5). Entries in the column labelled SR are factors bywhich the standard initial specific resistancesof Table 1 were multiplied. RED is the initial reduction factor applied to the specific resistances in the instability volume. RA4AX is maximum specific resistance reduction rate and PRM = A . RMAX(refer to eqn. (7)). TMAX is peak temperature rise for the entire grid. Actual predicted peak
122
JAMES H. K. YOO, THOMAS J. WHITE, DAVID L STONER
TABLE 2 PEAK TEMPERATURE
Run
No. 1
:
V, (volt) 4632 5500 2000
SR (ohm-cm)
(se:- ‘) (se:-‘) 260
4
150
ii ii
400
RISE AND PARAMETER
1000
RED
1.0
500 1500
VARIATION
RMAX (set-1)
0.5
PRM (“Cl)
1400
1000
0.2
9 10
5.0 0.75
i: 13 14
0.333
1:
TMAX (“Cat t msec)
(k%l)
19 (6) 4 (11)
21 ::
ii (18)
1044
ii i;
:‘: 21
I;] 13 (8)
14i 33 ::
2000 500
1.03 0.93 1.371
29 18 (5) (9) ;;
2:
500 100
0.713 0.815
21 (9)
22 31
temperature would be obtained by adding normal skin temperature 30-35°C to TMAX. The figures in parentheses after TMAX are times of occurrence in milliseconds of peak temperature rise. RS is the resistance of the instability volume at the time of peak temperature rise. Additional runs were made with varying instability radius. TMAX remained at 19 “C at 6 msec but RS, as expected, was different from the standard run, as shown in Table 3. TABLE 3 INSTABILITY VOLUME
RADIUS
(RIM),
AND
GRID
(RS’),
RESISTANCE SPACING
OF THE RADIAL
OF THE
IN THE
INSTABILITY
UNIFORM
PART
GRID (Ar)
RIM (cm)
RS at TMAX @ohm)
RS at t = 0 @ohm)
0.0665
49
240
0.09975 0.1425 0.3515
22 11 1.7
106 52 8.5
& 70
105 150 370
Figure 7 shows time dependence of PS (power in the instability), RS (instability resistance), T (temperature rise), and DAMAGE (fractional skin damage) for the standard run (No. 1). For the same run, Fig. 8 shows temperature rise with skin depth at 1, 6, and 45 msec. Between 0.75 and 1 mm and beyond 2.13 mm the profiles are estimates necessitated by the coarseness of the grid beyond O-75mm. Figure 9 shows the time variation of the electrical resistance profile on the z axis (r = 0), RES (I, 1). Numbers on the graph are times in milliseconds. For those runs
123
HUMANSKINBURNS
30. 6-
60RS
1.0. 2%
5-
50-
0.6- 20- 4-
40.
0.2-
IO-
s-
I-
10-4
10-3
10-Z
TIME bed
Fig. 7. /--
Instability history.
INSTABILITY -bj
DEPTH (mm)
Fig. 8.
Temporal and spatial temperature profiles
124
JAMES
H. K. YOO, THOMAS
1. WHITE,
DAVID
L. STONER
with peak temperature rise of about 20 “C and instability radius of 0.1 cm, energy expended in the instability is about 0.1 joule. Energy loss in the total skin resistance is about 30 joules and since the body resistance is about the same as total skin resistance, energy loss in the body resistance is also about 30 joules. When the instability radius is varied, energy per unit volume remains constant at about 15 joules/cm3. For run 2 (T = 4”C), instability energy is 0.02 joule while for run 4 (T = 48 “C), instability energy is 1 joule. Horny
J
IO7
-Basol
1 Vosculor
Fatty
0 =: 2 105:
-6
E -
104-
0,6,45 0,6,45 - 0,6,45 0,6,45
-
-45
IO3 0
2
I
i
DEPTH (mm) Fig. 9. Axial elemental resistances at 0,6 and 45 msec for the standard run. The horny and basal layers occupy grid points 17 through 20 and thevascular layer, points 21 through 29. Only 1.86 mm of the fatty layer is shown. The plot extends to the depth of the skin peripheral resistance used in the program, 2.74 mm.
In Fig. 7 it is seen that even after power in the instability has become negligible, the temperature is still high because of the time required for heat conduction and the assumed thermal resistance between the electrode and the horny layer. For this same reason, temperature rise in the instability volume does not follow the resistance history in this same volume, as may be seen by comparing Figs. 8 and 9. Although the resistance of the vascular layer is larger at 45 msec than the resistance of the instability volume, temperature rise is low due to the larger vascular volume. Conversely, even though fatty volume is again much larger, its resistance is so high that a possibly significant temperature rise occurs, 5-6°C.
HUMAN SKIN BURNS
125
For the standard run (No. l), the maximum applied potential occurs at 2 msec and results in maximum powers in the skin resistor and instability of 15 kwatt and 49 watt, respectively. Runs 1 and 12 through 16 involve specific resistance reduction rate parameters which result in the RR versus Tcurves shown in Fig. 10. The numbers on the curves are the corresponding run numbers. The table in Fig. 10 summarises the different parameters used.
16
0
5
IO
RUN NO.
RMIN
RMAX
PRM
13 12 15 14 18
.Ol .Ol .Ol .Ol .Ol
looo. 500. 2ooo. SW. 1000. loo.
1. .93 1.03 ,815 1.371 ,713
15
20
T Fig. 10.
Specific resistance
reduction
rates.
RATE - 0.9 RMAX AT T= 14 14 14 18 10 13
I 25
126
JAMES
H. K. YOO, THOMAS
I. WHITE,
DAVID
L. STONER
The parameter sensitivity results summarised in Table 2 lead to several general conclusions regarding the model’s applicability. (1) V,: Temperature rise is sensitive to the amplitude of the applied voltage. For the three runs with different VO: 1.8 TMAX=
&
(
(30)
1
The exponent 1.8 is near the expected value of about 2 for power dissipation in a resistive load for short times. The difference is due to heat conduction and the varying instability resistance. (2) z, CO:Source strength is proportional to power in the instability volume. For short times, temperature rise is also proportional to power if the instability resistance does not change appreciably due to current flow. Power is in turn proportional to the square of the applied voltage and we may therefore estimate the effect of varying r and o by comparing the integral: (31) for runs 4 through 7 with that for the standard run (No. 1). Table 4 shows the result of using eqn. (31) to predict maximum temperature rise. TABLE 4 MAXIMUM
TEMPERATURE
RISE
EQUATION
(TMAX)
PREDICTED
BY
(31)
Run
19 x M/M,,, (“C)
TMAX (“C)
1
19 34 11 :i
19 48
22 044
ti 19
:‘: 21
4 5 6 I
(kzl)
The second column is the predicted TMAX to be compared with the computer generated TMAXin the third column, The third and fourth columns are taken from Table 2. For all except run 4, the prediction is reasonably good, indicating that changes in TMAX are proportional to the changes in the integral in eqn. (3 1) caused by changes in z or o. For run 4, the small value of T extends the current pulse long enough for a significant reduction in the instability resistance to occur during the pulse and the computer generated TMAX is somewhat higher. More computer runs woul’d be required to completely describe the effect of o and w for z smaller than that of the standard run. (3) SR: Runs 8 and 9 were started with initialvalues of specific resistance reduced
HUMAN SKIN BURNS
127
by a factor of 5 and increased by a factor of 5, respectively, prior to the application of the instability reduction factor, RED. The power in the instability resistance is: PS =
(PO2 RS[l + RB(l/RS + l/RP)]2
(32)
When both RP and RS are scaled by a factor a, representing a change in all specific resistances, the instability power becomes: PS’ =
ci(POT)2 RS[a + RB(l/RS
+ l/RP)12
(33)
As a function of LX,Ps’ is maximised for: M= RB(l/RS
+ l/RP)
(34)
For run 1, RS = 52,000 ohm, RP = 170 ohm, RB = 159 ohm and a = 1. From eqn. (34), tl = 0.94 for these resistancevalues, showing that power in the insta,bility is nearly optimum for the standard run for maximising temperature rise when all other parameters remain unchanged. For runs 8 and 9, PS’was calculated from eqn. (33) with 01= 0.2 and 5. The ratios PS’/PS for these runs are 0.58 and O-53. When multiplied by the 19°C temperature rise for run 1, the results are 11 and 10°C respectively, in good agreement with the computer generated values of 12 and 11 “C. Changes in specific resistance from those of the standard run thus reduce TMAX. (4) RED: This factor is applied to the specific resistances in the instability just prior to the first time step and for the standard run is l/2. As a function of RS, PS in eqn. (32) is maximised for: RB. RP RB+RP
(35)
RP(POq2 4RB(RP + RB)
(36)
RS = and is: ps =
For RB = 159 ohm and RP = 170 ohm, RS = 82 ohm which corresponds to RED = l/645. Such’s largevalue of the reduction factor is considered neither realistic nor necessary for simulation of the instability. Runs 10 and 11 were made with RED = 314 and l/3, respectively. Since RS $ RP and RB, PS - l/RS for RED near l/2 and use of eqn. (32) gives predicted TMAX of 13 “C and 28 “C, respectively, for the two runs. For RED near i/2, TMAX is well estimated by: TMAX - 9.5/RED
(37) (5) RR: Runs 1 and 12 through 16 show very little variation in maximum temperature rise in spite of very large variation in the parameters of the specific resistance reduction rate, RR. In all the runs of Table 2,60 per cent of the instability
128
JAMES
H. K. YOO, THOMAS
J. WHITE,
DAVID
L. STONER
resistance is in the horny and basal layers, 40 per cent in the fatty layer and only 0.3 per cent in the vascular layer. Source strength is much higher in the horny and basal layers than in the fatty layer because of the smaller volume associated with this resistance. TMAX always occurs at a point near the base of the horny layer, about 75 p below the surface. The effect of resistance reduction due to temperature rise is therefore confined primarily to the outermost skin layers and has the effect of limiting temperature rise. This limiting effect is slightly less for those runs involving smaller values of RMAX (runs 13, 15 and 16) and slightly more for RMAX larger (run 12). For all runs except No. 16, RS is further reduced at time of TMAX to 40 per cent of its initial reducedvalue due to application of rate RR. For run 16, reduction is only to 70 per cent. (6) RIM: Increase of instability radius RIM has the simultaneous effect of an initial reduction RED and an elemental resistance volume increase both equal to (RZM)2/(0.1)2. The source strength therefore remains constant and the instability volume is so small relative to total model volume that TMAX is not affected. (7) Skin damage: Use of eqn. (29) to estimate skin damage results in a DAMAGE value of 1.0 for all runs except 2 and 9, for which DAMAGE is 0.009 and O-95, respectively. These results are for the centre of the instability near the skin surface. No data on electrical resistance heating reaction rates are available, although rates for radiative and conductive heating rates have been used by several authors. (See Weaver and Stoll(1969) for a typical application.) For DAMAGE = 1 - l/e, eqn. (29) requires 1 msec at T = 15 “C, 5.5 msec at 10 “C, and 670 msec at 5 “C steady temperature. (8) Thermal Parameters: While the effect of changing all conductivities or all volumetric specific heats by a constant factor is easily seen from eqn. (lo), the effect of differential changes between layers can be seen only by use of the computer program. No attempt has been made to vary thermal parameters as a function of temperature rise (Weaver and Stoll, 1969), although this is easily accomplished in the same manner as that used for specific resistance.
CONCLUSIONS
The simplified model of burn development presented here makes use of reasonable estimates of skin electrical characteristics and skinelectrode interface breakdown phenomena in a skin temperature rise computer program to achieve prediction of skin burns in agreement with observations and with previous estimates of skin damage criteria. The discussion above shows that predicted temperature rise is most sensitive to the initial charge on the defibrillator capacitor, to z (damping of defibrillator pulse) and to the initial reduction factor, RED. Since 7 is affected by the skin resistance and the internal resistance and inductance of the defibrillator discharge path, and since
HUMAN SKIN BURNS
129
these, in turn, affect the energy delivered to the heart muscle, a complete analysis of delivered waveforms is beyond the scope of the present paper. The initial reduction factor also affects 7. For those runs where skin electrical power can be simply estimated (no large variations in skin resistance): T=
Skin power x time Skin volume x volumetric specific heat
a good estimate of temperature rise. The possibility of skin damage due to the use of defective equipment as described here can obviously be reduced by lowering T, V, and SR, i.e. increased delivery time for the defibrillator pulse, reduced initial charge on the defibrillator capacitor and elimination, if possible, of the electrical resistance of the horny and basal skin layers from the discharge path. Any such reductions must, of course, be consistent with delivery of an effective defibrillation pulse to the heart. Areas which need more research and testing are modelling of the instability development, measurement of skin layer specific resistances, measurement of temperature rise in the skin layers, relation of temperature rise history to skin damage and measurement of skin layer thermal parameters as a function of temperature rise. The computational approach used here is sufficiently general that any variation of input parameters and model characteristics can be accommodated. gives
REFERENCES
CARSLAW, H. S. and JAEGER,J. C., Con&rim of heat in solids. (2nd ed.), Oxford at the Clarendon Press, 1959. DRAPER,J. W. and BOAG, J. W., The calculation of skin temperature distributions in thermography. Phys. Med. Biol., M(2) (1971) p. 201. GEDDES,L. A,. and BAKER,L. E., The specific resistance of biological material-A compendium of data for the biomedical engineer and physiologist. Med. Biol. Eng., 5 (1967) pp. 271-93. GEDDES,L. A. and BAKER,L. E., Principles of applied biomedical instrumentation, John Wiley, New York, 1968, PP. 155 and 164. HARDY. J. D.. STOLWIJK.J. A. J.. HAMMEL.H. T. and MURGATROYD. D. Skin temuerature and cuianeoui pain during warm water immeision. J. Appl. Physiol. 20 (i965) p. 1014. ‘ LIPKIN, M. and HARDY,J. D., Measurement of some thermal properties of human tissues. J. Appi. Physiol., 7 (1954) p. 212. MAINSTER,M. A., WHITE, T. J., TIPS, H. H. and WILSON, P. W., Transient thermal behaviour in biological systems, Bull. Math. Biophysics, 32 (1970) p. 303. MASTERSON, T. S. Private communication. Institute of Surgical Research, Burn Clinic, Brooke General Hospital, Fort Sam Houston, San Antonio, Texas, 1975. MIL;;;;, H. T., JR., The application of control theory tophysiologicalsystems, Saunders, Philadelphia, STOLWIJK,J. A. J. and HARDY,J. D. Skin and subcutaneous temperature changes during exposure to intense thermal radiation. J. Appl. Physiol., 20 (1965) p. 1006. STONER,D. L., Yoo, J. H. K., FELDTMAN,R. W. and STANFORD,W., Human skin burns induced by defibrillator default current. J. Thorac. and Cardiouas. Surg., 72(l) (1976), pp. 157-61. STRONG, P. Biophysical Measurements, Tektronix, Inc., Beaverton, Oregon, 1970.
130
JAMES
WEAVER, J. A.
and
H. K. YOO, THOMAS
STOLL,A.
J. WHITE,
DAVID
L. STONER
M., Mathematical model of skin exposed to thermal radiation. Aerospace
Medicine (January, 1969). p. 24.
Yoo, J. H. K. and STONER,D. L., Is the defibrillator compatible with electrosurgery units? Proc. 28th ACEMB, (1975) p. 136.
THOMASJ. WHITE Deparrment of Mathematics, San Antonio College, San Antonio, Texas (USA)
and DAVIDL. STONER Biomedical Service, Wilford Hall USAF Medical Center, San Antonio, Texas (USA)
(Received: 18 May, 1976)
SUMMARY
A preliminary model was preparedfor the development of superficial burns which can result from the interaction of jbulty dejibrillator and other electrical operating room equipment. A thermal model previously used to predict retina1 coagulation was adapted to give temperature rise history of the skin. Skin resistaizce is assumed to decrease rapidly due to the instability of a small volume near the surface. Heating due to dissipation of the dejibrillator capacitor energy in this skin volume is shown to be high enough to cause burns. A sensitivity analysis is included to isolate the dominant input parameters. Recommendations are made for future research.
SOMMAIRE
Un modgle prkliminaire simule les brtilures superficielles quipeuvent Otre causiespar une d$brillation dkfectueuse ou d’autres manipulations 6lectriques. Un modkle thermique, pr&ctdemment utilise pour predire la coagulation rktinienne, a 6t6 adapt& pour retracer l’histoire de l’tY&vationde temperature de la peau. La rksistance de la peau est supposie dkcroitre rapidement ir cause de l’instabilit6 d’un petit volume situ& prPs de la surface. L’&hau#ement dii ci la dissipation d’tnergie d’un dbjibrillateur au niveau de la peau est sufisamment important pour provoquer des brtilures. Une analyse fine est incluse dans ce travail, ajin de mettre en evidence les principaux parametres d’entrte. Des recommandations sont faites pour des recherches futures. 109 Int. J. Bio-Medical Computing (8) (1977)~-10 Applied Science Publishers Ltd, England, 1977
Printed in Great Britain
110
JAMES H. K. YOO, THOMAS J. WHITE, DAVID L. STONER
PARTIAL LIST OF SYMBOLS
AJ(J): AP: BF,.: BF, : C: C(J):
DAMAGE: ~FR: G: H: K: L: POT: PRM: PRS: ;. RB: RE(I): RES( I, J): RED: RIM: R,: RMIN: RMAX: RN: RP: RR: RS: RSKIN: S(r, z, 2): SR:
Annular ring area (cm2) Defined as the area for RP or rc (2.52 - RIM2) (cm2) Blood flow rate in the r direction (cmjsec) Blood flow rate in the z direction (cmjsec) Circuit capacitance (PFarad) Current through the Jth cylindrical annulus Fractional skin damage Overall approximate blood change rate (set- ‘) The total conductance of the instability volume (ohm-‘) Coefficient of surface heat transfer (watt cmm2 “C- ‘) Thermal conductivity (watt cm- r “C- ‘) Circuit inductance (henry) The applied defibrillator potential of the form, I’, exp (- rt) sin (WC) (volt) A . RMAX( ‘Cm’) used in reaction rate for eduction of specific resistance Voltage drop across total skin resistance (volt) Power expended in the instability volume (watt) Initial charge on C (Coulomb) Internal circuit resistance of the defibrillator pulse generator (ohm) The lumped body resistance (ohm) Resistance of the Zth annular ring segment around the skinbreakdown instability centre (ohm) Resistance of elemental layer I of annular ring J of the instability volume (ohm) The initial reduction factor applied to the specific resistances in the instability volume Instability radius (cm) Load resistance (ohm) The assumed minimum thermo-chemical reaction rate (set-‘) The assumed maximum thermo-chemical reaction rate (set- ‘) The radius ofthe r boundary for the modelvolume. It is assumed to be 2.5 cm The remaining skin resistance to current flow from the electrode, outside the instability area. Assumed to be 170 ohm in the model The reaction rate for reduction of specific resistance (set-‘) Resistance of the instability volume (ohm) Lumped resistance of skin to current flow (ohm) Source strength (heat generation source) confined to the instability volume (watt cme3) Specific resistance (ohm-cm)
HUMAN
T: TMAX: t: ZM:
AZ: PC:
Ar:
SKIN BURNS
111
minus the original skin Temperature rise, i.e. temperature temperature (“C) Maximum temperature rise (“C) Time elapsed since the voltage application (set) Half the total depth of the model volume, assumed to be 2.5cm Grid spacing in the uniform part of the axial grid (cm) Volumetric specific heat (joule cm- 3 ‘C- ‘) Grid spacing in the uniform part of the radial grid (cm)
INTRODUCTION
The expanding use of complex electronic instrumentation in the Operating Room and Intensive Care Unit is increasing the risk to the patient of severe burns and electrocution. Previous papers (Yoo and Stoner, 1975; Stoner et al., 1976) described skin burns resulting from the interaction of a faulty defibrillator and an electrosurgical unit with a grounded dispersive plate. In this case, a short developed between the defibrillator switching circuitry and the discharge circuit. The short grounded one side of the discharge capacitor, thereby allowing the defibrillator to discharge through any grounded points on the patient. Thus, the grounded electrosurgical dispersive plate provided a second discharge path for the defibrillator. The burns resulting from this second discharge path were 3-10 mm in diameter and occurred on the medial and posterior aspects of the right popliteal space and the upper two-thirds of the lower leg, as shown in Fig. 1 (Stoner et al., 1976). The number of burns correlated with the number of defibrillation attempts and these burns were caused by current arcing from the patient’s leg to the grounded electrosurgical dispersive plate. In an effort to more fully understand the burn process, a preliminary model for the development of superficial burns has been devised. Insight into the physiology of the burn process can now be gained in two ways. First, by describing physiological changes in skin histology with physical changes in engineering components, one can develop a better understanding of the effect of surgical procedures on the physiological properties of the skin. Secondly, modelling allows one to simulate a multitude of burn conditions as functions of a variety of parameters which could not realistically be duplicated experimentally. Adjustment of specific resistances in the skin volume corresponding to the observed burns is used to achieve predicted temperature rises of 1l-48 “C within the first 5 to 18 msec of the damped oscillatory defibrillator pulse. While damage may occur because of the passage of electrical current through the skin layers due to mechanisms other than heating, this range of temperature increments corresponds to those known to cause burns in skin and retina following absorption of incident radiation. While it is clear that skin burns are not a primary concern during open heart
112
JAMES H. K. YOO, THOMAS J. WHITE, DAVID L. STONER
Fig. I.
Photograph
of patient’s
burns immediately
after the incident.
surgery, reduction of additional patient insult is a valid concern. The thermalelectrical model described here could lead to a better understanding of the mechanism of electrical burns when extended to include a description of the mechanism for the initial reduction in skin resistance necessary to allow for high current flow through the outer layers. A detailed sensitivity analysis covering both skin and defibrillator parameters is included. This analysis shows that the behaviour of skin resistance with temperature rise and the initial reduction due to breakdown at the electrode-skin interface and in the horny and basal skin layers are controlling factors for adequate damage prediction.
PHYSICAL, THERMAL AND ELECTRICAL
MODEL
A typical electrode-skin interface equivalent circuit is given in Fig. 2 (Strong, 1970). Typical skin impedance value for a 1 cm2 plate is 10 kR + jo (10 pF), where o is the angular frequency (Strong, 1970). In this paper, however, the capacitive loading is
HUMAN SKIN BURNS
113
ignored assuming that the initial punch-through current surge across the skin-electrode interface diminishes the capacitive effect to negligible magnitude at the specific instability point. No interface contact impedance is included in this model. ELECTRODE INTERFACE
Fig. 2.
IMPEDANCE
Simplified typical electrode-skin interface equivalent circuit.
This simplification is also motivated by calculational convenience during the mathematical model analysis. Research on the electrical properties of skin is not sufficient to allow adequate modelling of human skin impedance. In addition, enormous individualvariations of human skin properties and unknown thermal and electrical characteristics of the skin-electrode contact dictate this simplification. Therefore, in this paper, only gross estimates for the skin properties could be used. Since no models of skin damage due to passage of electrical current seem to be available, skin properties used are those which have been measured in connection with studies of damage due to exposure to incident radiation in the ultraviolet, visible and infrared parts of the spectrum. The work of other authors may be consulted for extensive literature on this type of injury (Stolwijk and Hardy, 1965; Hardy et al., 1965; Lipkin and Hardy, 1954; Weaver and Stoll, 1969). The burns described in this paper are confined to the first few millimetres of skin depth so that the parameters used with radiant energy injury models should be accurate below this depth. The electrical and thermal properties must change radically with the development of the initial instability and consequent very high current flow. No model is as yet available for this breakdown in skin resistance andwe have therefore simulated this process by an initial reduction in specific resistance and have thereafter continued the reduction of specific resistance according to an assumed rate of change dependent on temperature rise. The main purpose of the analysis of predicted temporal and spatial temperature profiles in this paper is thus limited to (1) illustration of the development of skin
JAMES H. K. YOO, THOMAS J. WHITE, DAVID L. STONER
114
burns due to absorption of electrical energy in the skin resistance and (2) exploration of better ways to use electrical surgical aids to minimise patient injury. Table 1 summarises the physical, thermal and electrical properties of the skin-electrode model. Layer depths are taken from the work of Stolwijk and Hardy (1965) and Hardy et al. (1965). Skin thermal parameters, volumetric specific heat, pc, and conductivity, K, are taken from the work of Stolwijk and Hardy (1965) Hardy et al. (1965), Lipkin and Hardy (1954), Weaver and Stoll(1969) and Draper and Boag (197 1). Specific resistances, SR, are estimates based on the work of Geddes and Baker (1967, 1968). TABLE 1 SKIN-ELECTRODE THERMAL AND Layer
ELECTRICAL
Depth
K x lo3
(watt cm-
(cm) Bedding Electrode Horny Basal Vascular Fatty Muscle
PARAMETERS
2.07 0.41 0.01 0.01 0.07 1.11 1.32
0.600 1.59 4.35 4.40 4-31 2.12 4.83
‘C-’ )
5.00 2.50 2.09 3.98 4.19 2.10 4.10
SR
(ohm-cm)
10: 1:: 7000 1500
The electrode is modelled as a thin aluminium sheet ( w 50 cc)glued to a wood backing about 4mm thick. Detailed consideration of the thermal resistances at the skin-electrode boundary would include those due to skin-conductive grease, greasealuminium foil, aluminium foil-glue, and glue-wood backing. All of these have been combined into one thermal resistance estimated as 1 cm2 “C watt- ‘, There are no thermal resistances between internal layer surfaces. The electrode is supported by 2 cm thick bedding and there is 10 cm2 “C- ’watt- 1 thermal resistance between the bedding and the electrode wood backing. Thermal parameters for the electrode and bedding are also given in Table 1 and are based in part on the work of Lipkin and Hardy (1954). Figure 3 shows the physical arrangement of the skin model layers and Fig. 4 the
electrode bedding
I
backing I horny Fig. 3.
vascular layer
I basal
fatty layer
Physical model of the skin layers.
I
muscle
115
HUMAN SKIN BURNS
corresponding electrical resistance model. Simulation of the initial breakdown m skin resistance is confined to a cylinder of 3 mm depth and radius RIM. RIM is varied in the computer runs from 0.66 mm to 3.5 mm. In this volume, referred to hereafter as the instability volume, temperature rise, T,willvary with both depth and radius as a function of time. Therefore, in order to more accurately simulate further reduction of skin layer specific resistances in thisvolume due to temperature rise, the
Fig. 4.
Electrical
resistive model.
model cylinder is divided into 10 annular cylinders, the innermost with radius Ar/2 and the remaining 9 with radius Ar so that Ar = RIMi9.5. Each cylindrical annulus is further divided into 16 layers: 2 each in the horny and basal skin layers, 9 in the vascular and 3 in the fatty skin layers. Each of these 160 elemental volumes comprising the instability volume has resistance RES(Z, J) computed as:
RES(Z,J) = SR(Z,J)AZ(I)/AJ(J)
(1)
where SR(Z,J)is specific resistance at depth Z(Z), AZ(I) is resistance layer depth, AJ(J) is annular ring area and I and J are axial and radial computer program indices, respectively. The resistance of annular ring J is then: 16
RE(J) = (1 I=1
SR(I, J> AZ(I)
Y
NJ),
J=
1,2,...,10
(2)
The total conductance of the instability volume is: 10
G= c
[l/W41
(31
J=l
In Fig. 4 RP represents the remaining skin resistance to current flow from the electrode, assumed to be 5 cm in diameter. RP is thus the resistance of a cylinder of inner radius RIM and outer radius 2.5 cm and is in parallel with the resistance of the instability volume. Specific resistance in RP is assumed not to be affected by temperature rise. RP may be calculated from:
116
JAMES
H. K. YOO, THOMAS
J. WHITE,
DAVID
L. STONER
RP = (O*OlSRl + O.OlSR2 + 0*07SR3 + 0.21SR4)jAP
(4)
where SRI, SR2, SR3 and SR4 are initial specific resistances for the horny, basal, vascular and fatty layers, respectively. AP = 7r(2*.5’- RIM*) is the area for RP. For the model used here RP is about 170 ohm. Since electrical resistance between the electrode and the horny layer is zero, the remaining resistance in the electrical path is the lumped body resistance RB (Fig. 4). RB is modelled as a cylinder 5 cm in radius, 50 cm long, with specific resistance of 250 ohm-cm. For these values: RB = 1590hm RB is not affected by temperature
rise in the skin. The applied defibrillator potential, POT in Fig. 4, has the form: POT = V,, exp ( - rt) sin (or) (volt)
(5)
This form is a lumped constant form of the voltage appearing across the load in a series RCL circuit. The three constants have been simplified as: r = (R + RJ2L ~~{&(!s~)‘-’ V, = QR,/LCw
where R, is load resistance; R is internal circuit resistance; L is circuit inductance; C is circuit capacitance; Q is initial charge on C, and w is real for underdamping. A reasonable form of reaction rate for reduction of specific resistance with increasing temperature is suggested by Milhorn (1966) and Masterson (1975) as the solution of: dRR = A. RR(RMAX
dT
where RR(T) is rate with maximum RMAX, The solution, which includes T = 0, is:
- RR)
A is a constant and T temperature rise.
RMAX
RR=
l+
RMAX
- RMZN
RMZN
(7) exp(-PRM.
T)
where RMZN, the minimum rate, occurs at T = 0, RMAX is the asymptotic rate at T = co and PRM = A. RMAX. Each of the resistances RES(Z, J) in the unstable volume is reduced with increasing temperature at this rate by integration of: dSR. dt
-=
-RR.SR
(8)
117
HUMANSKINBURNS
or: (9)
SR = SROexp( - ~ORRP-(*)ldl)
where SR, is initial reduced specific resistance and T(t) is temperature rise as a function of time for the elemental resistance. COMPUTATIONALDETAILSOFTHEMATHEMATICALMODEL
Figure 5 shows the physical model used in the machine computations. The equations to be solved for temperature rise, T, are:
(10) T(r, z, 0) = 0
(11)
R(RN, z, t) = T(r, 0, t) = T(r, 2ZM, t) = 0
(12)
where RN (2.5cm) is the radius of the r boundary, ZM (2.5cm) is half the total depth, and S(r, z, t) is source strength in watt cm -3 due to current flow confined to the instability volume. The method of solution is described in Mainster et al. (1970). instability
electrode
Fig. 5.
Schematic
three-dimensional
computational
model
The total volume used for heat conduction is about 100cm3 and is considered adequate for simulation of an infinite medium when thevolume of primary interest is less than lo-* cm3. The grid used has 3 1 points in the r direction and 37 points in the z direction with 20 uniform grid spaces in r and 16 uniform axial grid spaces centred at the middle of the z axis. The first point in the horny layer has index 17. AZin the uniform portion of the grid is 50~. Ar in the uniform r grid is RIMl9.5 when the instability contains 10 grid points in the radial direction. Grid stretching factors are 1.7 for the z grid and 1.55 for r when RIM = 1 mm. (Varies with RIM). (a) Thermal resistance and discontinuous conductivity Treatment of discontinuities in conductivity between layers and layer boundary
118
JAMES H. K. YOO, THOMAS J. WHITE, DAVID L. STONER
thermal resistance is achieved by modification of the approximating difference equation coefficients for z (Mainster et al., 1970). Figure 6 shows temperature rise at four grid points near a surface of separation. Heat flux must be continuous in the z direction at the boundary: +I,
= --Kg1
(13)
where n indicates normal derivative and indices 1 and 2 refer to opposite sides of the surface. Ti_l, Ti, Ti + 1 and Ti+2 are temperature rises at grid points zi_ I, zi, zi+ 1 and Zi+z, respectively. iij;and iii+, are temperature rises on either side of the surface
i i+l i+2 Kit-+-, Ki+t
i-l Fig. 6.
Media boundary with Newton cooling.
of separation of the two media of different conductivities, Ki and Ki+ 1. We assume Newton cooling at the surface with thermal resistance l/H (cm2 “C watt- ‘). Then: and :
-&+1
‘ii+1- K+1 (zi+l
_
zi)
=
H(‘ii+1- T>=A+1
(15)
2 Fij- Ti
= H(‘ii-K(zi + 1 - zi> 2
T+,> =fi
(14)
These two equations can be solved for the normal heat fluxes with the result that:
-fi
=A+1 =
KiKi+,(Ti+l - ri)
(Ki+Ki+l)(i’+;-z’)+5+
(16)
HUMANSKINBURNS
119
In the difference equations use of this result at z index i gives: 2Ki zi
+
I
-
zi
-
1
Ki+,)(Zi+l
-Zi)
2
Ki+I
while at z index i+ 1:
Ti - Ti_,
Ti+, - Ti +
+&KC+,
-
Zi-Zi-1
H
I -(Ki+
Kit,)(Zit,
-Zi>
Ki
These difference equations (eqns. (17) and (18)) are used at the bedding-electrode backing boundary and at the electrode-horny layer boundary. When there is perfect thermal contact as between skin layers, l/H = 0, T = T+ 1 and the term with H drops out of the equations. Where the grid is uniform and there is no thermal resistance ((l/H) = 0) and no conductivity difference (Ki = Kit 1), eqns. (17) and (18) reduce to: __ (d)Z2
,+,-2Ti+Ti_,} iT,
(19)
and :
as expected for centred second differences.
(b) Bloodflow in vascular layer Although the defibrillator pulse is short, an estimate of the effect of blood flow in the vascular layer has been included in the interest of possible further applications of the program. Since the co-ordinate system is cylindrical, simulation of blood flow cooling is restricted to axial and radial flow. The fluxes in the r and z directions when allowance is made for fluid motion are: f,=
-K&cT.BF,
(21)
and :
f,= -Kg
+ pcT. BF,
(22)
where BF,. and BF, are flow rates in the r and z directions, respectively, in cm/set. The terms involving derivatives with respect to r in eqn. (10) become:
120
JAMES
-
H. K. YOO, THOMAS
fi(rf,)=!A&rKg (
d!?+”
ar
- ii(rpcT.
J. WHITE,
DAVID
L. STONER
BF,)
>
(K!? ar) (“F “irF)TP+)?$ -pc
___C+__L
(2%
and :
BF, has been assumed to be zero and flow in the radial direction has been used to approximate blood flow parallel to the skin layers. BF, is assumed to go to zero like r as r + 0 so that: BF,+, dBF = constant r ar In view of the overall approximation, BF,/r is assumed to be this same constant and eqn. (10) is modified to: K --pcr.BFR r
$+
Kg+:
-pcBFR.T=pcg-S(r,z,l)
(25)
where BFR in units of set-’ is to be interpreted as an overall approximate blood change rate. In all computations, BFR = 5 set- ‘. Partial difference coefficients are modified from those reported by Mainster et al. (1970) to include the additional terms in eqn. (25).
S(r, z, t) The voltage drop across the total skin resistance is:
(c) Source strength
PRS = RP. POT] [RP + RB(l
+ G. RP)]
(26)
where G is given in eqn. (3), RP in eqn. (4) and POTin eqn. (5). The current through the elemental resistances with index J is: C(J) = PRSIRE(J)
(27)
and : S(Z, J, t> = RES(Z, 4[C(412/
[Az(Z)AJ(J)1
cw
is the source strength in watt cm -3 for elemental resistance RES(Z, J). Source strengthvaries with time through thevariation of POT, RES(Z, .I), RE(J) and C(J). Equations (7) and (9) are used with averagevalue of temperature rise at grid point I, J during time step k to reduce the specific resistance in eqn. (1). The RES(Z, 4 from eqn. (1) then changes the source strength (eqn. (28)) for use in time step k + 1
HUMAN SKIN BURNS
121
so that the source term closely follows the temperature history in the instability volume according to the parameters chosen for the reaction rate of eqn. (7). (d) Initial instability Since, the appearance of the burn is associated with random locations on the electrode and since the mechanism of initial breakdown in the resistance of the outer skin layers is not understood, the initial instability has been simulated by reduction of all specific resistances in the assumed instability volume prior to the start of calculation of temperature rise by the factor RED. As indicated above, the instability volume is made to correspond to the observed size and shape of the defibrillator leakage current burns, a few millimetres in radius and depth. The breakdown in resistance used to estimate the damaging current flow is confined to this volume. (e) Damage prediction
A skin damage prediction calculation is made a part of the computer program using a rate also defined by eqn. (7). Skin fractional damage is defined by DAMACE=l.-exp(-[/R[T(t)]dt)
(29)
Parameter values used are RMAX = 1000 set- ‘, RMZN = 0.01 set- ‘, and A = 0.0009 set “C - ‘. (f) Digital computing machine requirements
All runs were made on the San Antonio College IBM 370-158 computer by Mr Robert Lay of the Data Processing Department Staff. In Fortran H, storage required is 228K bytes for compilation. CPU running time averages about 20 set for 34 time steps or about l/2 msec/time step/grid point for the 1147 (37 x 3 1) point grid used.
RESULTS
Table 2 summarises the results of use of the program and model described in the previous sections. Run No. 1 was considered as the standard and succeeding runs of Table 2 show only those parameters which were changed. ‘V,, r, and o refer to defibrillator potential as shown in eqn. (5). Entries in the column labelled SR are factors bywhich the standard initial specific resistancesof Table 1 were multiplied. RED is the initial reduction factor applied to the specific resistances in the instability volume. RA4AX is maximum specific resistance reduction rate and PRM = A . RMAX(refer to eqn. (7)). TMAX is peak temperature rise for the entire grid. Actual predicted peak
122
JAMES H. K. YOO, THOMAS J. WHITE, DAVID L STONER
TABLE 2 PEAK TEMPERATURE
Run
No. 1
:
V, (volt) 4632 5500 2000
SR (ohm-cm)
(se:- ‘) (se:-‘) 260
4
150
ii ii
400
RISE AND PARAMETER
1000
RED
1.0
500 1500
VARIATION
RMAX (set-1)
0.5
PRM (“Cl)
1400
1000
0.2
9 10
5.0 0.75
i: 13 14
0.333
1:
TMAX (“Cat t msec)
(k%l)
19 (6) 4 (11)
21 ::
ii (18)
1044
ii i;
:‘: 21
I;] 13 (8)
14i 33 ::
2000 500
1.03 0.93 1.371
29 18 (5) (9) ;;
2:
500 100
0.713 0.815
21 (9)
22 31
temperature would be obtained by adding normal skin temperature 30-35°C to TMAX. The figures in parentheses after TMAX are times of occurrence in milliseconds of peak temperature rise. RS is the resistance of the instability volume at the time of peak temperature rise. Additional runs were made with varying instability radius. TMAX remained at 19 “C at 6 msec but RS, as expected, was different from the standard run, as shown in Table 3. TABLE 3 INSTABILITY VOLUME
RADIUS
(RIM),
AND
GRID
(RS’),
RESISTANCE SPACING
OF THE RADIAL
OF THE
IN THE
INSTABILITY
UNIFORM
PART
GRID (Ar)
RIM (cm)
RS at TMAX @ohm)
RS at t = 0 @ohm)
0.0665
49
240
0.09975 0.1425 0.3515
22 11 1.7
106 52 8.5
& 70
105 150 370
Figure 7 shows time dependence of PS (power in the instability), RS (instability resistance), T (temperature rise), and DAMAGE (fractional skin damage) for the standard run (No. 1). For the same run, Fig. 8 shows temperature rise with skin depth at 1, 6, and 45 msec. Between 0.75 and 1 mm and beyond 2.13 mm the profiles are estimates necessitated by the coarseness of the grid beyond O-75mm. Figure 9 shows the time variation of the electrical resistance profile on the z axis (r = 0), RES (I, 1). Numbers on the graph are times in milliseconds. For those runs
123
HUMANSKINBURNS
30. 6-
60RS
1.0. 2%
5-
50-
0.6- 20- 4-
40.
0.2-
IO-
s-
I-
10-4
10-3
10-Z
TIME bed
Fig. 7. /--
Instability history.
INSTABILITY -bj
DEPTH (mm)
Fig. 8.
Temporal and spatial temperature profiles
124
JAMES
H. K. YOO, THOMAS
1. WHITE,
DAVID
L. STONER
with peak temperature rise of about 20 “C and instability radius of 0.1 cm, energy expended in the instability is about 0.1 joule. Energy loss in the total skin resistance is about 30 joules and since the body resistance is about the same as total skin resistance, energy loss in the body resistance is also about 30 joules. When the instability radius is varied, energy per unit volume remains constant at about 15 joules/cm3. For run 2 (T = 4”C), instability energy is 0.02 joule while for run 4 (T = 48 “C), instability energy is 1 joule. Horny
J
IO7
-Basol
1 Vosculor
Fatty
0 =: 2 105:
-6
E -
104-
0,6,45 0,6,45 - 0,6,45 0,6,45
-
-45
IO3 0
2
I
i
DEPTH (mm) Fig. 9. Axial elemental resistances at 0,6 and 45 msec for the standard run. The horny and basal layers occupy grid points 17 through 20 and thevascular layer, points 21 through 29. Only 1.86 mm of the fatty layer is shown. The plot extends to the depth of the skin peripheral resistance used in the program, 2.74 mm.
In Fig. 7 it is seen that even after power in the instability has become negligible, the temperature is still high because of the time required for heat conduction and the assumed thermal resistance between the electrode and the horny layer. For this same reason, temperature rise in the instability volume does not follow the resistance history in this same volume, as may be seen by comparing Figs. 8 and 9. Although the resistance of the vascular layer is larger at 45 msec than the resistance of the instability volume, temperature rise is low due to the larger vascular volume. Conversely, even though fatty volume is again much larger, its resistance is so high that a possibly significant temperature rise occurs, 5-6°C.
HUMAN SKIN BURNS
125
For the standard run (No. l), the maximum applied potential occurs at 2 msec and results in maximum powers in the skin resistor and instability of 15 kwatt and 49 watt, respectively. Runs 1 and 12 through 16 involve specific resistance reduction rate parameters which result in the RR versus Tcurves shown in Fig. 10. The numbers on the curves are the corresponding run numbers. The table in Fig. 10 summarises the different parameters used.
16
0
5
IO
RUN NO.
RMIN
RMAX
PRM
13 12 15 14 18
.Ol .Ol .Ol .Ol .Ol
looo. 500. 2ooo. SW. 1000. loo.
1. .93 1.03 ,815 1.371 ,713
15
20
T Fig. 10.
Specific resistance
reduction
rates.
RATE - 0.9 RMAX AT T= 14 14 14 18 10 13
I 25
126
JAMES
H. K. YOO, THOMAS
I. WHITE,
DAVID
L. STONER
The parameter sensitivity results summarised in Table 2 lead to several general conclusions regarding the model’s applicability. (1) V,: Temperature rise is sensitive to the amplitude of the applied voltage. For the three runs with different VO: 1.8 TMAX=
&
(
(30)
1
The exponent 1.8 is near the expected value of about 2 for power dissipation in a resistive load for short times. The difference is due to heat conduction and the varying instability resistance. (2) z, CO:Source strength is proportional to power in the instability volume. For short times, temperature rise is also proportional to power if the instability resistance does not change appreciably due to current flow. Power is in turn proportional to the square of the applied voltage and we may therefore estimate the effect of varying r and o by comparing the integral: (31) for runs 4 through 7 with that for the standard run (No. 1). Table 4 shows the result of using eqn. (31) to predict maximum temperature rise. TABLE 4 MAXIMUM
TEMPERATURE
RISE
EQUATION
(TMAX)
PREDICTED
BY
(31)
Run
19 x M/M,,, (“C)
TMAX (“C)
1
19 34 11 :i
19 48
22 044
ti 19
:‘: 21
4 5 6 I
(kzl)
The second column is the predicted TMAX to be compared with the computer generated TMAXin the third column, The third and fourth columns are taken from Table 2. For all except run 4, the prediction is reasonably good, indicating that changes in TMAX are proportional to the changes in the integral in eqn. (3 1) caused by changes in z or o. For run 4, the small value of T extends the current pulse long enough for a significant reduction in the instability resistance to occur during the pulse and the computer generated TMAX is somewhat higher. More computer runs woul’d be required to completely describe the effect of o and w for z smaller than that of the standard run. (3) SR: Runs 8 and 9 were started with initialvalues of specific resistance reduced
HUMAN SKIN BURNS
127
by a factor of 5 and increased by a factor of 5, respectively, prior to the application of the instability reduction factor, RED. The power in the instability resistance is: PS =
(PO2 RS[l + RB(l/RS + l/RP)]2
(32)
When both RP and RS are scaled by a factor a, representing a change in all specific resistances, the instability power becomes: PS’ =
ci(POT)2 RS[a + RB(l/RS
+ l/RP)12
(33)
As a function of LX,Ps’ is maximised for: M= RB(l/RS
+ l/RP)
(34)
For run 1, RS = 52,000 ohm, RP = 170 ohm, RB = 159 ohm and a = 1. From eqn. (34), tl = 0.94 for these resistancevalues, showing that power in the insta,bility is nearly optimum for the standard run for maximising temperature rise when all other parameters remain unchanged. For runs 8 and 9, PS’was calculated from eqn. (33) with 01= 0.2 and 5. The ratios PS’/PS for these runs are 0.58 and O-53. When multiplied by the 19°C temperature rise for run 1, the results are 11 and 10°C respectively, in good agreement with the computer generated values of 12 and 11 “C. Changes in specific resistance from those of the standard run thus reduce TMAX. (4) RED: This factor is applied to the specific resistances in the instability just prior to the first time step and for the standard run is l/2. As a function of RS, PS in eqn. (32) is maximised for: RB. RP RB+RP
(35)
RP(POq2 4RB(RP + RB)
(36)
RS = and is: ps =
For RB = 159 ohm and RP = 170 ohm, RS = 82 ohm which corresponds to RED = l/645. Such’s largevalue of the reduction factor is considered neither realistic nor necessary for simulation of the instability. Runs 10 and 11 were made with RED = 314 and l/3, respectively. Since RS $ RP and RB, PS - l/RS for RED near l/2 and use of eqn. (32) gives predicted TMAX of 13 “C and 28 “C, respectively, for the two runs. For RED near i/2, TMAX is well estimated by: TMAX - 9.5/RED
(37) (5) RR: Runs 1 and 12 through 16 show very little variation in maximum temperature rise in spite of very large variation in the parameters of the specific resistance reduction rate, RR. In all the runs of Table 2,60 per cent of the instability
128
JAMES
H. K. YOO, THOMAS
J. WHITE,
DAVID
L. STONER
resistance is in the horny and basal layers, 40 per cent in the fatty layer and only 0.3 per cent in the vascular layer. Source strength is much higher in the horny and basal layers than in the fatty layer because of the smaller volume associated with this resistance. TMAX always occurs at a point near the base of the horny layer, about 75 p below the surface. The effect of resistance reduction due to temperature rise is therefore confined primarily to the outermost skin layers and has the effect of limiting temperature rise. This limiting effect is slightly less for those runs involving smaller values of RMAX (runs 13, 15 and 16) and slightly more for RMAX larger (run 12). For all runs except No. 16, RS is further reduced at time of TMAX to 40 per cent of its initial reducedvalue due to application of rate RR. For run 16, reduction is only to 70 per cent. (6) RIM: Increase of instability radius RIM has the simultaneous effect of an initial reduction RED and an elemental resistance volume increase both equal to (RZM)2/(0.1)2. The source strength therefore remains constant and the instability volume is so small relative to total model volume that TMAX is not affected. (7) Skin damage: Use of eqn. (29) to estimate skin damage results in a DAMAGE value of 1.0 for all runs except 2 and 9, for which DAMAGE is 0.009 and O-95, respectively. These results are for the centre of the instability near the skin surface. No data on electrical resistance heating reaction rates are available, although rates for radiative and conductive heating rates have been used by several authors. (See Weaver and Stoll(1969) for a typical application.) For DAMAGE = 1 - l/e, eqn. (29) requires 1 msec at T = 15 “C, 5.5 msec at 10 “C, and 670 msec at 5 “C steady temperature. (8) Thermal Parameters: While the effect of changing all conductivities or all volumetric specific heats by a constant factor is easily seen from eqn. (lo), the effect of differential changes between layers can be seen only by use of the computer program. No attempt has been made to vary thermal parameters as a function of temperature rise (Weaver and Stoll, 1969), although this is easily accomplished in the same manner as that used for specific resistance.
CONCLUSIONS
The simplified model of burn development presented here makes use of reasonable estimates of skin electrical characteristics and skinelectrode interface breakdown phenomena in a skin temperature rise computer program to achieve prediction of skin burns in agreement with observations and with previous estimates of skin damage criteria. The discussion above shows that predicted temperature rise is most sensitive to the initial charge on the defibrillator capacitor, to z (damping of defibrillator pulse) and to the initial reduction factor, RED. Since 7 is affected by the skin resistance and the internal resistance and inductance of the defibrillator discharge path, and since
HUMAN SKIN BURNS
129
these, in turn, affect the energy delivered to the heart muscle, a complete analysis of delivered waveforms is beyond the scope of the present paper. The initial reduction factor also affects 7. For those runs where skin electrical power can be simply estimated (no large variations in skin resistance): T=
Skin power x time Skin volume x volumetric specific heat
a good estimate of temperature rise. The possibility of skin damage due to the use of defective equipment as described here can obviously be reduced by lowering T, V, and SR, i.e. increased delivery time for the defibrillator pulse, reduced initial charge on the defibrillator capacitor and elimination, if possible, of the electrical resistance of the horny and basal skin layers from the discharge path. Any such reductions must, of course, be consistent with delivery of an effective defibrillation pulse to the heart. Areas which need more research and testing are modelling of the instability development, measurement of skin layer specific resistances, measurement of temperature rise in the skin layers, relation of temperature rise history to skin damage and measurement of skin layer thermal parameters as a function of temperature rise. The computational approach used here is sufficiently general that any variation of input parameters and model characteristics can be accommodated. gives
REFERENCES
CARSLAW, H. S. and JAEGER,J. C., Con&rim of heat in solids. (2nd ed.), Oxford at the Clarendon Press, 1959. DRAPER,J. W. and BOAG, J. W., The calculation of skin temperature distributions in thermography. Phys. Med. Biol., M(2) (1971) p. 201. GEDDES,L. A,. and BAKER,L. E., The specific resistance of biological material-A compendium of data for the biomedical engineer and physiologist. Med. Biol. Eng., 5 (1967) pp. 271-93. GEDDES,L. A. and BAKER,L. E., Principles of applied biomedical instrumentation, John Wiley, New York, 1968, PP. 155 and 164. HARDY. J. D.. STOLWIJK.J. A. J.. HAMMEL.H. T. and MURGATROYD. D. Skin temuerature and cuianeoui pain during warm water immeision. J. Appl. Physiol. 20 (i965) p. 1014. ‘ LIPKIN, M. and HARDY,J. D., Measurement of some thermal properties of human tissues. J. Appi. Physiol., 7 (1954) p. 212. MAINSTER,M. A., WHITE, T. J., TIPS, H. H. and WILSON, P. W., Transient thermal behaviour in biological systems, Bull. Math. Biophysics, 32 (1970) p. 303. MASTERSON, T. S. Private communication. Institute of Surgical Research, Burn Clinic, Brooke General Hospital, Fort Sam Houston, San Antonio, Texas, 1975. MIL;;;;, H. T., JR., The application of control theory tophysiologicalsystems, Saunders, Philadelphia, STOLWIJK,J. A. J. and HARDY,J. D. Skin and subcutaneous temperature changes during exposure to intense thermal radiation. J. Appl. Physiol., 20 (1965) p. 1006. STONER,D. L., Yoo, J. H. K., FELDTMAN,R. W. and STANFORD,W., Human skin burns induced by defibrillator default current. J. Thorac. and Cardiouas. Surg., 72(l) (1976), pp. 157-61. STRONG, P. Biophysical Measurements, Tektronix, Inc., Beaverton, Oregon, 1970.
130
JAMES
WEAVER, J. A.
and
H. K. YOO, THOMAS
STOLL,A.
J. WHITE,
DAVID
L. STONER
M., Mathematical model of skin exposed to thermal radiation. Aerospace
Medicine (January, 1969). p. 24.
Yoo, J. H. K. and STONER,D. L., Is the defibrillator compatible with electrosurgery units? Proc. 28th ACEMB, (1975) p. 136.
