Purpose
The modern permanent prostate brachytherapy procedure, introduced three decades ago [1], features transperineal deposition of radioactive sources under ultrasound-guidance. Superior disease control has been achieved due to accurate enactment of an implant plan that was developed either pre- or intra-operatively with computerized treatment planning system (TPS). While the general principles of brachytherapy treatment planning have been covered elsewhere [2], the current paper reviews the evolution of computerized treatment planning for brachytherapy in the United States. There were of course developments outside the U.S. during this same era, and the interested reader may learn more on these advances through examining the following references [3-6]. Although the authors recognize the contributions of European and Asian workers to early advances in brachytherapy dosimetry, the scope of this paper is limited to the authors’ area of our expertise – developments in the U.S.
Material and methods
This study is based upon a literature review and original interviews, supported by contemporaneous literature, that were obtained from the following physicists and engineers: Lowell Anderson (March 2005), Stephen Balter (June 2009 and February 2013), Philip Heintz (December 2012), Ted Jackson (October 2012), Michael Saylor (August 2013), Marilyn Stovall (February 2013).
Working in the dark
In the earliest days, brachytherapists relied on atlases, tables, and experience to determine source distribution, implant strength, and duration. By the middle of the 20th century, interstitial brachytherapy (more properly, brachyradium [7]) dosimetry relied on two source distribution systems, which were developed in the 1930s: Paterson-Parker and Quimby. The Paterson-Parker System, developed in Manchester’s Holt Radium Institute, specified an inhomogenous spatial distribution of sources to achieve a relatively homogenous dose distribution [8,9]. The system developed by physicist Edith Quimby at New York’s Memorial Hospital stipulated a homogenous distribution of sources to generate an inhomogenous dose distribution [10,11]. Both systems utilized pre-planning (i.e., preceding clinical implantation) to determine the strength and geometric distribution of implanted sources needed to deliver a specified dose (typically the minimum dose) to a target.
The procedure used at New York’s Memorial Hospital during the 1950s and 1960s is described as an example. The brachytherapist was called to the operating theater to implant a tumor that had just been deemed unresectable (a common occurrence in an era without advanced imaging). He had, at hand, a caliper and a supply of 222Rn seeds. The total strength of implanted sources was calculated by the following process: the three dimensions of the target were measured with a ruler and the volume calculated by using the appropriate formula for the target shape. The maximally tolerated dose for the target volume was determined by consulting a Memorial chart (larger volumes were found to tolerate lower doses). Next, the Paterson-Parker or Quimby tables were used to determine the total milligram-hours of radium (i.e., mg-h 226Ra, an obsolete unit of source strength) necessary to achieve the dose for that volume; the result was then multiplied by the appropriate Paterson-Parker elongation-correction factor. Since most Memorial implants of that era used 222Rn seeds, the mg-h 226Ra value had to be divided by the 222Rn mean life (i.e., 132 h) to determine the total implant strength. This value was divided by the current strength of each seed to determine the number of seeds to be implanted. Finally, a decision had to be made as to how to spatially distribute the seeds [12].
This process proved too cumbersome (especially as the patient’s chest or abdomen had already been opened), so the following shortcut was implemented: the three measured dimensions were averaged. To determine the 222Rn source strength to implant, this average dimension was multiplied by 10 (an empirically derived factor, specific to 222Rn); a factor of 5 was used when 125I seeds were introduced [12]. Thereafter, a Quimby-style (homogenous distribution of sources) implant was performed, which was appropriately fast and simple. Physicist Lowell Anderson further simplified the process through the introduction of nomograms for 125I and 103Pd [13,14].
While the use of atlases, tables, and shortcuts systematized the process, the implants were often seriously flawed. It was difficult (even for expert brachytherapists) to reproduce the ‘ideal’ source geometry specified by the systems, and the achieved ‘minimum’ target dose typically fell well below the prescribed dose [15]. Larger tumors were frequently underdosed, a problem that was exacerbated by the introduction of 125I. Initially, only the inverse-square law was used in calculations; this was not an issue for 226Rn, 222Rn, 198Au, or 192Ir due to their high-energy γ-ray emissions, but was inadequate when calculating the attenuation of low-energy photons emitted by radionuclides such as 125I [16]. Underdosing was not recognized because post-implant dosimetry was not initially performed. Manual calculation of dose at more than a few points was tedious, especially for a seed implant where the position of more than 50 sources had to be identified and their dose distributions had to be plotted and summed.
What was needed was a system to identify a large number of sources, then calculate, sum, and display a spatial representation of the dose distribution over the entire volume of interest (not just the periphery). Computers have the capacity to calculate dose distribution rapidly enough, so that timely adjustments could be made post-operatively (for afterloaded implants), or even intra-operatively (for permanent implantation). This historical review explores the introduction of computer dosimetry to brachytherapy.
Introduction of computerized dosimetry at Memorial Hospital
Kia-Chi Tsien (1916-1993) introduced a TPS for multi-field teletherapy using ‘tabulating machines’ in 1955 [17,18]. The depth-dose characteristics, specified as a function of field size and focus-to-skin distance, for all 60Co, orthovoltage, and megavoltage equipment in use at Memorial Hospital were stored in a library of punch cards (Fig. 1). The appropriate punch cards for the selected fields were entered into the tabulating machine, which summed their contributions and plotted the isodose lines relative to an external contour of the body. This system quadrupled the efficiency of dosimetrists, allowing selection of the most favorable of several plans in the same time that it would have taken to run a single plan manually. Richard Nelson and Mary Lou Meurk extended this work to perform the first computerized brachytherapy dose calculations in 1958 [19]. The positions of implanted sources were localized through stereo-shift radiographs and assigned a location at the nearest point on a three-dimensional Cartesian coordinate lattice with 5 mm spacing (Fig. 2). The location of each source was transferred to a punch card, and the tabulating machine summed the contribution of all sources to plot the dose distribution. Output was represented as a matrix of points with 1 cm spacing; isodose lines were drawn by hand. By 1961, the system was programed in the FORTRAN code for running on a time-shared model 7090 IBM (International Business Machines, Corp., Armonk, NY, USA) computer system, which was “sufficiently fast that its results have been used to modify interstitial implants before and during treatment” [20]. While the system was feasible for 226Ra needles and 192Ir ribbons in afterloading tubes, localization of individual seeds was “possible, but not always practicable” as the program accommodated an “equally spaced array of seeds” [20].
Stephen Balter (born in 1940) revised the system shortly after joining the Memorial Physics staff in 1963. He had learned how to program the UNIVAC (UNIVersal Automatic Computer) by Eckert-Mauchly Computer, Corp. (Philadelphia, PA, USA) while he was an undergraduate student at the Polytechnic Institute of Brooklyn (B.S. in Physics, 1961), and had access to Columbia University’s computer laboratory during graduate work (M.S. in Radiologic Physics, 1963). He learned the FORTRAN II computer programming language from a manual while commuting to work by subway. It was his intention to write a program for Memorial’s model CDC-160A (Control Data Corp., Minneapolis, MN, USA) 4k 12-bit minicomputer to spatially display the dose distribution in any plane for a large seed implant, to calculate average and minimum doses, and to generate dose-volume histograms [15,21]. The program had to be broken into modules that ran sequentially because computer memory could not store the entire program. First, seed location was determined by the stereo-shift method. This was a tedious task because implants often contained more than 50 seeds (Fig. 3). After the source locations were entered, dose was calculated to points on a matrix. Dose calculation points in the vicinity of the target volume were typically 1 cm apart, so there were 1,000 dose points for a 10 cm3 volume. It took about 0.1 s to calculate dose to a point from a single source, 100 s to calculate the dose to all dose points from a single source, and over an hour to calculate the sum of contributions from all sources to all dose points. A printout of dose to points on any plane was generated, but isodose lines were still drawn by hand (Fig. 4). The situation improved when the Memorial Physics Department installed a stereographic reference frame (Fig. 5) to maintain fixed geometry for stereoshift or orthogonal radiographs, a digitizer to enter seed location, a program to automatically digitize seed location, and a model 1800 IBM 16k computer that increased calculation speeds.
Although the ‘Memorial Implant Dosimetry Application System’ (‘MIDAS’) was freely shared, few institutions had the computing power to utilize it. Beginning in 1967, the Memorial Physics Department offered a computational service, whereby other institutions submitted implant source strength and spatial distributions. Dose distributions were calculated using a model 360/67 IBM central processor on the General Electric (GE, Corp., Fairfield, CT, USA) commercial timeshare service, and results were returned to the participating institution by teletype [22]. Initially, six institutions utilized the service [23], but more than 60 hospitals participated by 1972 [24].
Introduction of computer dosimetry at M. D. Anderson Hospital
Similar events transpired concurrently in Houston. Unlike at New York’s Memorial Hospital where 222Rn seeds were commonly implanted, the radiotherapists at M. D. Anderson Hospital preferred 226Ra needles and 198Au seeds. Upon graduation from Baylor University in 1951 with a B.S. degree in Mathematics and Zoology, Marilyn Stovall (born in 1931) answered a newspaper advertisement for employment at Anderson. Robert Shalek (1923-2007), acting Head of Physics, offered her a dosimetrist position with on-the-job training if she would commit to remaining there for at least a year. She accepted his offer and has been there for more than six decades.
Her initial responsibilities included machine calibration, treatment planning, and checking dose delivery using a Sievert chamber in the radiation field. As custodian of the 226Ra sources, she was present at the time of implantation and also when the sources were removed (i.e., at any hour of the day or night). After several years, her position was advanced to ‘Applied Mathematician’, and she was asked to write computer code for brachytherapy dose calculations. She learned the FORTRAN computer programming language at the University of Houston’s night school, and wrote a program with the assistance of a professor at Texas A&M; punch cards were transported between Houston and College Station, TX on passenger bus. The program was ready in 1960, but the only usable computer, a RAMAC-model IBM, was in the business office. The program was first used to calculate post-implant dosimetry for radium needle implants using orthogonal x-ray images, then still commonly referred to as ‘roentgenograms’. Source localization was done by transverse tomography. The computer determined the square of the distance to points on a 1 mm grid, assigned doses based on a tabular lookup, and summed the doses [25]; isodose curves were drawn manually. It took an hour for each plane to be calculated, at which point a new set of cards would be fed into the computer. As dose points were spaced at 0.5 cm (there were 8 times as many dose calculation points per volume as had been utilized at Memorial), the entire process initially took 10 hours. As the computer was used for billing during the day, Stovall set up a cot in the business office and ran post-implant dosimetry at night. This situation continued until the Anderson Physics Department acquired its own faster computer, a model 1620 IBM, which could perform the same task within minutes.
In a retrospective analysis of implants, Stovall had determined that most recurrences and complications (e.g., tissue necrosis) occurred in regions of under- or overdosage [26], respectively. Timely dose calculation, with rapid depiction by plotter or cathode ray tube, allowed compensation for less than ‘ideal’ geometry by varying the timing of needle removal for interstitial implants and through source selection for afterloading implants [27].
Similar to Memorial’s MIDAS dosimetry program, Anderson’s RADCOMP (RADiation COMPutation) dosimetry program, was in the public domain [28]. As few institutions had the computing capacity to utilize it, the Anderson Department of Physics also had a thriving multi-institutional computational service.
Advances in software algorithms
After the introduction of computerized brachytherapy treatment planning, advancements beyond the point-source approximation were introduced to account for the elongated physical distribution of the radionuclide [29,30]. Cooperating with Shalek at the Anderson Department of Physics, Meisberger et al. addressed the issue of radiation attenuation and scatter in the medium through dose correction factors that were expressed as simple polynomial functions [31]. While initially used only for high-energy photon-emitting sources such as 198Au, 192Ir, 137Cs, 226Ra, and 60Co in the late 1960s, this work was later extended to low-energy photon-emitting sources such as 125I, where the radiological influence of medium is more pronounced [32]. Another important influence on clinical dose calculations in this era was Monte Carlo methods for radiation transport simulations, where the first full 3D model of a brachytherapy source was performed in 1971 [33]. Results from these types of Monte Carlo simulations were simplified for integration into existing planning algorithms and development of future systems.
Development of proprietary software
Peter Wootton (1924-2004), departed Anderson to join the Physics faculty at the University of Washington. One of his doctoral students, Philip Heintz (born in 1943), had learned computer programming as an undergraduate engineering student at California Polytechnic State University from 1961 to 1965. Together with physicist Douglas Jones, they used elements of the Memorial and Anderson public domain programs to generate the TPS at the University of Washington. Heintz took the program with him when he went into private practice in California in 1971.
Heintz joined Ralph Worsnop in 1985 to develop MARx Plan, a program written in the PASCAL computer programming language to be run on a personal computer (PC). The partnership with Worsnop dissolved in 1987 and Heintz had the program re-written in the Microsoft C (Microsoft, Corp., Redmond, WA, USA) computer programming language to be run on the IBM AT platform. His new company SSGI (Small Systems Group, Inc., Chico, CA, USA) introduced the Prowess 2000 TPS software in 1988. In addition to teletherapy planning, it had modules for brachytherapy with 226Ra and 137Cs tubes, 192Ir and 125I seeds, and 182Ta wire. In the early 1990s, the program was rewritten for the Windows (Microsoft, Corp., Redmond, WA, USA) operating system and was renamed Prowess 3000.
Another commercial brachytherapy TPS was the ROCS program (Radiation Oncology Computer Systems, Inc., Carlsbad, CA, USA) introduced in 1988 by physicist, founder, and president David Neblett. It was written in the BASIC computer programming language for the PC and included an ultrasound-guided prostate brachytherapy package.
William Saylor was trained as a medical physicist at the Naval Postgraduate School in Monterey, CA in 1962. In 1974, he introduced Therpac (Pinecliff Associates, Wilmington, NC, USA), the first complete TPS written for the model HP 9830A (Hewlett-Packard, Corp., Palo Alto, CA, USA) microcomputer. His son, Michael Saylor, joined him in 1981 and rewrote the program in Microsoft C computer programming language for the PC; Therpac-PC was introduced in 1986. In 1993, a transperineal ultrasound-guided implant (TUI) module (Therpac Plus 6.5 B3Dtui) was introduced. Within a few years, the module could incorporate ultrasound images for pre-planning and CT images for post-implant dosimetry. William Saylor retired in 1996, and sold the company to Multimedia Medical Systems (MMS), Inc. (Charlottesville, VA, USA), a telemedicine company specializing in video-conferencing and image-sharing. Michael Saylor rewrote the TUI module with graphical-user interface for the Windows operating system, introduced in 1997. MMS’s chief engineer, Ted Jackson (Biomedical Engineering M.Sc. in 1990 and Ph.D. in 1994 from the University of Virginia) collaborated with the physicians and physicists at the Mount Sinai Hospital in New York to develop a user-defined dose-optimization module, and the facility for interactive intra-operative planning with real-time display of an implanted dose distribution. Two years later, the MMS program was acquired by Varian Medical Systems, Inc. (Palo Alto, CA, USA) and renamed ‘VariSeedTM’; it has been the dominant program in use in the U.S. for interstitial brachytherapy treatment planning for over a decade.
Conclusions
Computerized brachytherapy treatment planning software was initially developed in the Physics Departments of Memorial Hospital (by Nelson, Meurk, and Balter) and M. D. Anderson Hospital (by Stovall and Shalek). These systems were public-domain, and available to smaller hospitals through multi-institutional computational services. Commercial brachytherapy treatment planning programs, designed to run on microcomputers computers, were introduced in the late 1980s and early 1990s. These systems brought computer brachytherapy dosimetry into the clinic and surgical theatre.
Disclosure
Authors report no conflict of interest.
References
1. Holm HH, Juul N, Pedersen JF et al. Transperineal Iodine-125 seed implantation in prostatic cancer guided by transrectal ultrasonography. J Urol 1983; 130: 283-286.
2. Rivard MJ, Venselaar JLM, Beaulieu L. The evolution of brachytherapy treatment planning. Med Phys 2009; 36: 2136-2153.
3. Rosenwald JC. Repérage des sources radioactives en curiethérapie. Utilisation d’un ordinateur pour calculer leur position. Ann Phys Biol en Méd 1970; 3: 139-150.
4. Worthley BW, Nicholls RL. Computer evaluation of Sievert’s integral by recursive formula. Phys Med Biol 1972; 17: 854-865.
5. Chassagne D, Gerbaulet A, Dutreix A et al. Utilisation pratique de la dosimétrie par ordinateur en curiethérapie gynécologique. J Radiol Electrol 1977; 58: 387-393.
6. van der Laarse R. Computerized radiation treatment planning. Ph.D. thesis, Antoni van Leeuwenhoek Hospital, Amsterdam, The Netherlands (1981). Available through the collection of the International Nuclear Information System of the International Atomic Energy Agency, Vienna, Austria, vol. 13(5), reference 13656920; http://www.iaea.org/inis/collection/NCLCollectionStore/_Public/13/656/13656920.pdf
7. Forssell G. La lutte social contre le cancer. J Radiologie 1931; 15: 621-634.
8. Paterson R, Parker HM. A dosage system for gamma-ray therapy. Brit J Radiol 1934; 7: 592-612.
9. Paterson R, Parker HM. A dosage system for interstitial radium therapy. Brit J Radiol 1938; 11: 252-266.
10. Quimby EH. Physical factors in interstitial radium therapy. Am J Roentgenol 1935; 33: 306-316.
11. Quimby EH. Dosage table for linear radium sources. Radiology 1944; 43: 572-577.
12. Henschke UK, Cevc P. Dimension averaging a simple method for dosimetry of interstitial implants. Rad Biol Ther 1968; 9: 287-298.
13. Anderson LL. Spacing nomographs for interstitial implants of I-125 seeds. Med Phys 1976; 3: 48-51.
14. Anderson LL, Moni JV, Harrison LB. A nomograph for permanent implants of palladium-103 seeds. Int J Radiat Oncol Biol Phys 1993; 27: 129-135.
15. Balter S, Freed BR, Ragazzoni GD et al. An extension of the Memorial system for implant dosimetry. Radiology 1966; 87: 475-482.
16. Anderson LL, Aubrey RF. Computerized dosimetry for I-125 prostate implants. In: Hilaris BS, Batata MA (eds.). Brachytherapy Oncology – 1983. New York, NY: Memorial Sloan-Kettering, 1983; 57-63.
17. Tsien KC. The application of automatic computing machines to radiation treatment planning. Brit J Radiol 1955; 28: 432-439.
18. Tsien KC. A study of basic external radiation treatment techniques with the aid of automatic computing machines. Brit
J Radiol 1958; 31: 32-40.
19. Nelson RF, Meurk ML. The use of automatic computing machines for implant dosimetry. Radiology 1958; 70: 90.
20. Laughlin JS, Siler WM, Holodny EI et al. A dose description system for interstitial radiation therapy. Am J Roentgenol Radium Therap Nuclear Med 1963; 89: 470-490.
21. Holt G, Hilaris B, Balter S et al. Experience with computerized implant dosimetry. Am J Roentgenol 1968; 102: 688-693.
22. Holt JG, Balter S, Laughlin JS et al. Experience with a dose distribution computation service. Ann NY Acad Sci 1969; 161: 344-447.
23. Balter S. Memorial implant dosimetry application system. Radiology 1969; 92: 161-162.
24. Randall G, Balter S, Holt JG et al. The Memorials implant dosimetry automated system. Computer Prog Biomed 1972; 2: 137-152.
25. Shalek RJ, Stovall MA. The calculation of isodose distributions in interstitial implantations by a computer. Radiology 1961; 76: 119-120.
26. Fletcher GH, Stovall M. A study of the explicit distribution of radiation in interstitial implantations: II. Correlation with clinical results in squamous-cell carcinomas of the anterior two-thirds of tongue and floor of mouth. Radiology 1962; 78: 766-782.
27. Shalek RJ, Stovall M. The computation of dosage in interstitial and intracavitary radiation therapy. J Chron Dis 1966; 19: 519-522.
28. Stovall M, Shalek RJ. The M. D. Anderson method for the computation of isodose curves around interstitial and intracavitary radiation sources: III. Roentgenograms for input data and the relation of isodose calculations to the Paterson-Parker system. Am J Roentgenol 1968; 102: 677-687.
29. Powers WE, Bogardus Jr CR, White W et al. Computer estimation of dosage of interstitial and intracavitary implants. Radiology 1965; 85: 135-142.
30. Shalek RJ, Stovall M. The M. D. Anderson method for the computation of isodose curves around interstitial and intracavitary radiation sources: I. Dose from linear sources. Am J Roentgenol 1968; 102: 662-672.
31. Meisberger MM, Keller RJ, Shalek RJ. The effective attenuation in water of the gamma rays of gold-198, iridium-192, cesium-137, radium-226 and cobalt-60. Radiology 1968; 90: 953-957.
32. Thomadsen BR, Williamson JF, Rivard MJ et al. Anniversary Paper: Past and current issues, and trends in brachytherapy physics. Med Phys 2008; 35: 4708-4723.
33. Krishnaswamy V. Calculation of the dose distribution about californium-252 needles in tissue. Radiology 1971; 98: 155-160.
