Conic tangency equations and Apollonius problems in - RISC

14 Pages·2002·0.27 MB·English

MathematicsandComputersinSimulation61(2003)101–114
Conic tangency equations and Apollonius problems in
biochemistry and pharmacology
Robert H. Lewisa,∗, Stephen Bridgettb
aDepartmentofMathematics,FordhamUniversity,Bronx,NY10458,USA
bQueensUniversity,Belfast,Ireland
Received11April2002;accepted25April2002
Abstract
TheApolloniusCircleProblemdatestoGreekantiquity,circa250b.c.Giventhreecirclesintheplane,findor
constructacircletangenttoallthree.Thiswasgeneralizedbyreplacingsomecircleswithstraightlines.Viéte[Canon
mathematicusseuAdtriangulacumadpendicibus,Lutetiae:ApudIoannemMettayer,Mathematicistypographum
regium, sub signo D. Ioannis, regione Collegij Laodicensis, p. 1579] solved the problem using circle inversions
before1580.Twogenerationslater,Descartesconsideredaspecialcaseinwhichallfourcirclesaremutuallytangent
toeachother(i.e.pairwise).Inthispaper,weconsiderthegeneralcaseintwoandthreedimensions,andfurther
generalizationswithellipsoidsandlines.Webelieve,wearethefirsttoexplicitlyfindthepolynomialequationsfor
theparametersofthesolutionsphereinthesegeneralizedcases.Doingsoisquiteachallengeforthebestcomputer
algebrasystems.Wereportlatersomecomparativetimesusingvariouscomputeralgebrasystemsonsomeofthese
problems.Wealsoconsiderconictangencyequationsforgeneralconicsintwoandthreedimensions.
Apolloniusproblemsareofinterestintheirownright.However,themotivationforthisworkcameoriginallyfrom
medicalresearch,specificallytheproblemofcomputingthemedialaxisofthespacearoundamolecule:obtaining
thepositionandradiusofaspherewhichtouchesfourknownspheresorellipsoids.
©2002IMACS.PublishedbyElsevierScienceB.V.Allrightsreserved.
Keywords:Conictangencyequations;Apolloniusproblems;Medialaxis;Polynomialsystem
1. Introduction
PerhapsthebestintroductiontotheApolloniusCircleProblemisinCourantandRobbins([7],pp.125
and161).Theydescribethecircleinversiontechnique,amongotherapproachestothegeneralandspecial
problems.ManypeopleworkedontheApolloniusCirclequestionsinthe20thcentury,includingBoyd
[3], Coxeter [6], Kasner and Supnick [11], and Pedoe [15]. Frederick Soddy, a Nobel Prize winner in
∗ Correspondingauthor.
E-mailaddress:[email protected](R.H.Lewis).
0378-4754/02/$–seefrontmatter©2002IMACS.PublishedbyElsevierScienceB.V.Allrightsreserved.
PII:S0378-4754(02)00122-2
102 R.H.Lewis,S.Bridgett/MathematicsandComputersinSimulation61(2003)101–114
Fig.1.Specialcase,twodimensions,whereallcirclesaremutuallytangent;solutioncircle:innermost.
chemistryin1921,expressedasolutiontothespecialcaseasatheoremintheformofapoem,“TheKiss
Precise,” which was published in the journal Nature [18]. Soddy proved that for four mutually tangent
circlesthecurvaturesarerelatedby
2(κ2+κ2+κ2+κ2) = (κ +κ +κ +κ )2
1 2 3 4 1 2 3 4
(whichwasknowntoDescartes)andsimilarlyforn+2mutuallytangentn-spheresin(n+1)-space:
(cid:1) (cid:3) (cid:1) (cid:3)
(cid:2)n+2 (cid:2)n+2 2
n κ2 = κ
i i
i=1 i=1
Recently,Lagariasetal.[13]havewrittenapaperdescribingpackingsofsuchcircles/spheresinhyper-
bolic n-space and other geometries. Roanes-Lozano and Roanes-Macias[17] have a different approach
tosomeofthequestionsweaskhere(seeFigs.1and2).
Fig.2.Generalcase,twodimensions,solutioncircleistangenttotheoriginalthree,whicharearbitrarilyplaced.Ingeneralthere
areeightsolutions,basedoneachoftheoriginalsbeinginsideoroutsidethesolution.
R.H.Lewis,S.Bridgett/MathematicsandComputersinSimulation61(2003)101–114 103
Fig.3.Left:medialaxis(lines)foraverysimpletwo-dimensionalmolecule(grayatoms)withonemedialdiscshown(dashed
circle).Right:thediameteroftheenclosing(dashed)circlearoundthismoleculequicklyconfirmsthatitistoobigtofitintothe
pocketinthemoleculeontheleft.
2. Biochemicalmotivation
Apollonius problems are of interest in their own right. However, the motivation for this work came
originally from medical research, specifically the problem of computing the medial axis of the space
around a molecule: obtaining the position and radius of a sphere which touches four known spheres or
ellipsoids(Fig.3).
Alllifeformsdependoninteractionsbetweenmolecules.X-raycrystallographyandnuclearmagnetic
resonance(NMR)yieldthethree-dimensionalstructureofmanylargeandimportantproteins.Thecoor-
dinates of atoms within these molecules are now available from public databases and can be visualized
oncomputerdisplays.
Withinthesemoleculestherearereceptorsites,whereacorrectlyshapedandcorrectlychargedhormone
ordrug(oftenreferredtoasa“ligand”)canattach,distortingtheoverallshapeofthereceptormolecule
to cause some important action within the cell. This is similar to the concept of the correct key fitting
intoalock.However,inbiochemistry,the“lock”oftenchangesshapeslightlytoaccommodatethekey.
Moreover, there are often several slightly differently shaped keys that fit the lock, some binding to the
receptorstronglyandothersbindingweakly.
Forinstance,theshapeoftheblood’shaemoglobinmoleculemeansanoxygenmoleculebindsweakly
toit,sotheoxygencanbereleasedwhereneeded,whereascarbonmonoxide,beingaslightlydifferent
shape,bindsverystronglytothesamesiteonhaemoglobin.Thiscanbefatal.
2.1. Automateddockingalgorithms
Identifyingwhichnaturallyoccurringligandmoleculesfitintowhichreceptorisaveryimportanttask
in biochemistry and pharmacology. Moreover, in rational drug design the aim is often to design drugs
whichbetterfitthereceptorthanthenaturalhormone,soblockingtheactionofthenaturalhormone.The
designoftheACEinhibitorusedtoreducebloodpressureisagoodexampleofthis.
Theuseofcomputeralgorithmstoselectapotentialdrugfromthousandsofknownligandsforaparticu-
larreceptorisreceivingmuchresearchinterestbyseveralgroups.Someaimsforsuchanalgorithminclude
104 R.H.Lewis,S.Bridgett/MathematicsandComputersinSimulation61(2003)101–114
• helpingtoidentifypotentialbindingsitesinproteins,
• estimatingifaligandcanreachthebindingsite,
• predictingtheligand’smostsuitableorientationinthesite,
• determininghowgoodafittheligandwouldbe,comparedtootherknownligands.
Software such as DOCK 4.0 [8] is being used by pharmaceutical companies. For larger molecules,
methodshavebeeninvestigatedsuchasAckermannetal.[1]whichuseanumericalbest-fitestimate.
2.2. Themedialaxis
We present here another approach, the “medial axis”. The medial axis is defined as the locus of the
center of a maximal disc (in two-dimensional), or sphere (in three-dimensional) as it rolls around the
interiororexteriorofanobject.Itiseffectivelyaskeletonoftheoriginalobject,andhasbeenfoundhelpful
incomputeraidedengineeringformeshing,simplifyinganddimensionallyreducingcomponents.
Inthemedialaxisofamoleculethemedialedgesorsurfacesareequidistantfromtheatomswhichthe
spheretouches,andtheradiusofthemedialspherevariestobemaximal(Fig.3).Thisisdifferentfrom
theConellysurfaces,whereaprobesphereofconstantsizeisused.Astheimagesbelowdemonstrate,the
medialaxisapproachleadstotheApolloniustypequestionsofspherestouchingotherspheres,ellipsoids,
orlines.
3. Mathematicalapproach
WeareinterestedinthegeneralApolloniusproblem,asolutioncircle(orsphere)tobetangenttothe
originalones,whicharearbitrarilyplaced.Wewantto“know”thesolutionsymbolicallyandalgebraically.
Thatis,fortheradiusandforeachcoordinateofthecenter,wewantanequationexpressingitintermsof
thesymbolicparametersofthethree(four)givencircles,ellipses,orlines(spheres,ellipsoids).Ideally,
theequationwillbeoflowdegreesothatifnumericalvaluesarepluggedinforthesymbolicparameters
of the given shapes, a relatively simple one variable equation will result. We will see later that in some
caseswecanattainthisideal,butinotherswecompromisebyassumingthatthelinesorellipsoidsarein
certainimportantbutspecialorientations.
Webeginwithasystemofpolynomialequations
f (x,y,z,a,b,...) = 0
1
f (x,y,z,a,b,...) = 0
2
f (x,y,z,a,b,...) = 0
3
...
expressing the intersection and tangency conditions. Suppose that x, y, and z are the desired variables
of the solution circle (sphere) and that a,b,... are the parameters of the given figures. For each of x,
y, z in turn we want to derive from the system {fi} a single equation, the resultant, containing only it
and the parameters a,b,... To accomplish this elimination, we use either Gröbner Bases or methods
from the theory of resultants. Gröbner Bases are fairly well known (see for instance [2]), but here is
a brief description: from a set of polynomials, for example, the {fi} above, we generate an ideal by
R.H.Lewis,S.Bridgett/MathematicsandComputersinSimulation61(2003)101–114 105
taking all sums of products of other polynomials with the elements in the set. Then analagous to the
basis of a vector space, a certain set of polynomials that are “better” generators than the original {fi}
canbeproduced.IfthisGröbnerBasisischosencarefully,itwillcontainthedesiredresultant.Resultant
methodsarelesswellknown,especiallytheapparentmethodofchoice,theBezout–Cayley–Dixon–KSY
method[4,12,14],whichwedescribeinsubsequentsections.
We shall compare solutions done with computer algebra systems Maple 6, Mathematica, CoCoA 4.0
[5],andFermat2.6.4[9].
4. TheCayley–Dixon–Bezout–KSYresultantmethod
Hereisabriefdescription.Moredetailsarein[12].
To decide if there is a common root of n polynomial equations in n−1 variables x,y,z,... and k
parametersa,b,...
f (x,y,z,...,a,b,...) = 0
1
f (x,y,z,...,a,b,...) = 0
2
f (x,y,z,...,a,b,...) = 0
3
...
• Create the Cayley–Dixon matrix, n × n, by substituting some new variables t1,...,tn−1 into the
equationsinacertainway.
• Compute cd = determinant of the Cayley–Dixon matrix (a function of the new variables, variables,
andparameters).
• Form a second matrix by extracting the coefficients of cd relative to the variables and new variables.
Thismatrixcanbelarge;itssizedependsonthedegreesofthepolynomialsfi.
• Ideally, let dx = the determinant of the second matrix. If the system has a common solution, then
dx = 0.dxinvolvesonlytheparameters.
• Problem: the second matrix need not be square, or might have det = 0 identically. Then the method
appearstofail.
• However,wemaycontinue[12,4]:findanymaximalranksubmatrix;letksybeitsdeterminant.Exis-
tenceofacommonsolutionimpliesksy = 0.
• ksy = 0 is a sufficient equation, but one must be aware of spurious factors, i.e. the true resultant is
usually a small factor of ksy. Indeed, finding the small factor without actually computing ksy is very
desirable.
5. Two-dimensionalApolloniusresults
Problem1. Threecirclesgiven,findafourthtangenttoallthree.
Firststep:GivenonecircleAwithcenter(ax,ay)andradiusarfindequation(s)thatasecond“solution”
circleS,(sx,sy,sr)mustsatisfytobetangent.Wecouldbeginbyletting(x,y)bethepointofintersection.
Then(x,y)isonthefirstcircleiff(x−ax)2+(y−ay)2 = ar2,andonthesecondiff(x−sx)2+(y−sy)2 =
106 R.H.Lewis,S.Bridgett/MathematicsandComputersinSimulation61(2003)101–114
Fig.4.Twotangentcircles.
sr2 (see Fig. 4). Using derivatives, tangency gives a third equation, so we may eliminate x,y. This first
stepgivestheequationthatthesixparametersmustsatisfyfortangency.
However, the desired first step equation is geometrically obvious without going through the process
inthepreviousparagraph:fromFig.4,thepointoftangencyisonthelineconnectingthecenters.Each
centerisonthecirclewiththeother’scenterandthesumofthetworadii.Thus,theequationwewant,
isobviously
(sx−ax)2+(sy−ay)2 = (sr+ar)2
So,giventhreecirclesA,B,andC,thesolutioncircleparameterssx,sy,srmustsatisfythreeequations:
(sx−ax)2+(sy−ay)2 = (sr+ar)2
(sx−bx)2+(sy−by)2 = (sr+br)2
(sx−cx)2+(sy−cy)2 = (sr+cr)2
These are the equations for the case of interest to us, when each of the original circles is outside the
solutioncircle.Otherpossibilitiesinvolvechangingsomeplusestominusesabove.
Secondstep(thesolution):Thereareseveralwaystoproceed.Aspointedoutin[7],onemayexpand
thepreviousequationsandsubtracttoeliminatesquareterms,thensolveforsomevariablesandsubstitute
back. Then too, one may simplify the equations by assuming that one of the circles, say A, is centered
at the origin, as this is simply a translation, and furthermore that ar = 1, as this is just a change of
scale.However,forconsistencyandcomparisonoflatermoredifficultcases,letusproceedwiththefully
symbolic,Dixon–KSYmethodwithoutsubstitutinganyconstants.Weusedthecomputeralgebrasystem,
Fermat 2.6.4 [9] running on a Macintosh Blue and White G3 at 400MHz. Think of the above as three
equationsinthe“variables”sxandsrandthe“parameters”syplustheninea,b,cparameters.Eliminate
sx,srandobtaintheresultantinsy,i.e.theequationsymustsatisfyintermsoftheoriginaldata.Ittakes
about 3MB of RAM, 0.6s. The answer, an irreducible polynomial, has 593 terms, is degree two in sy.
Similarresultsobtainforsxandsr.
Problem2. Threeellipsesgiven,findacircletangenttoallthree.
R.H.Lewis,S.Bridgett/MathematicsandComputersinSimulation61(2003)101–114 107
Fig.5.Alleightsolutioncircles(blackandgreen)tangenttothreegivenellipses(red).
Fig.8.Sixsolutionconicstangenttofivegivenellipses.
108 R.H.Lewis,S.Bridgett/MathematicsandComputersinSimulation61(2003)101–114
Again,thequestionisofinterestinitsownrightmathematically.Biochemically,itmaybeusefultomodel
atomsbyellipsesinsteadofcircles.Weassumetheellipsesareallparalleltotheaxes(seeSection7.1).
WritetheequationofellipseAandcircleS,thenthetangencyequations:
a2(x −ax)2+a2(y −ay)2 = a2a2
2 1 1 2
(x −sx)2+(y −sy)2 = sr2
a2(y −ay)(x −sx)−a2(x −ax)(y −sy) = 0
1 2
Firststep:JustasinProblem1,wemusteliminatex,y fromtheabove.Butitisnotsoeasynow!The
pointoftangencyneednotlieonthelineconnectingthecenters.
Using again Fermat/Dixon/Mac G3, the elimination is completed using 2.1s and about 3MB RAM.
Theirreducibleanswerhas696termsintheparametersax,ay,a ,a ,sx,sy,sr.Recallthatthisfirststep
1 2
wastrivialinthethreecirclescaseandfiteasilyononeline.
Secondstep:Forafullysymbolicsolution,wewouldnowproducethreesuch696termpolynomials,
oneforeachofthreegivenellipsesA,B,C.Thenwewouldgetanequationfor,say,srbyeliminatingsx
andsy.However,alreadyintheDixon–KSYalgorithmattheearlystepofcomputingcd,thedeterminant
of the Cayley–Dixon matrix, Fermat exhausted all 700MB of RAM on the machine. Therefore, instead
of a fully symbolic solution, we made up numerical examples of three ellipses, with all parameters
integer.Eachofthethreeequationsnowhadonlythethreevariablessr,sx,sy,between36and92terms.
We applied Dixon–KSY to get the equation for sr. The second matrix was 88×88, with each entry a
polynomialovertheintegersinsronly.Therankwas72.Thedeterminantcomputation(byonemethod)
took92min,andthefinalanswerhad169terms.
However, more realistic examples will probably have parameters that are decimals or “floats”, not
integers,withsomeerrortolerance.So,alternatively,weforgoedtheDixon–KSYmethodintheprevious
paragraph. We took the three equations in sr,sx,sy with 36–92 terms and solved numerically with a
multivariate Newton’s method. This converges in a few seconds. One example is shown in Fig. 5, with
alleightsolutioncircles.
6. Three-dimensionalApolloniusresults
Problem3. Fourspheresgiven,findaspheretangenttoall.
Analogoustothecirclesintwo-dimensionalspace:
(sx−ax)2+(sy−ay)2+(sz−az)2 = (sr+ar)2
(sx−bx)2+(sy−by)2+(sz−bz)2 = (sr+br)2
(sx−cx)2+(sy−cy)2+(sz−cz)2 = (sr+cr)2
(sx−dx)2+(sy−dy)2+(sz−dz)2 = (sr+dr)2
Byelimination,getanequationforsrintermsofthea,b,c,d variables,etc.Theanswerhas18,366
termsandisquadraticinsr.
R.H.Lewis,S.Bridgett/MathematicsandComputersinSimulation61(2003)101–114 109
Several people used their favorite CAS to work on this problem. Cejtin [10] used a Gröbner Basis
library written in Scheme. He used an elimination-ordering. To get the equation for sr took just over
1h and about 50MB of RAM. Doing it for the special case where ax, ay and az are all = 0 took 46s
and used a bit under 8MB of RAM. Solving for sx took 1.75h CPU time and used 60MB RAM. In the
ax = ay = az = 0specialcaseittook300sCPUtimeandusedunder13MBofRAM.Allofthesetimes
areonan400MHzPentiumIImachine.
Coauthor Stephen Bridgett ran Mathematica with the simplified version of the equations (i.e. ax =
ay = az = 0) with the Solve command. Computing sy took 4.6h and 41MB RAM. It could not do the
fullproblem.
Israel[16]ranMaple6with(Gröebner)onaSunSPARCSolaris.Maplegaveupafter42min,using
647MBRAM.Onthesimplifiedproblemax = ay = az = 0itcompletedin14minusing574MBRAM.
CoauthorLewistriedthefullproblemwithCoCoA4.0onthesameMacintoshasabove.UsingGröbner
basesandtheElimcommand,thefullproblemwassolvedforsrin95minand29MBRAM.
Coauthor Lewis tried the full problem with Fermat/Dixon/Mac G3. Computing sr took 1.7s and
5MB of RAM. Parameters sx, sy, and sz took essentially the same amounts of time and space. This
is with Fermat 2.6.4, released in January 2002. These times are approximately a 6-fold improvement
over the earlier version of Fermat, and a 10-fold improvement in space. Note that even the times
from an earlier version of Fermat (10s) are two orders of magnitude faster than any competitor
above.
Problem4. Fourellipsoidsgiven,findaspheretangenttoall.
Similar to two-dimensional case of three ellipses. Define the ellipsoid by semi-axes a ,a ,a , center
1 2 3
(ax,ay,az). (So, the ellipsoid is parallel to the axes. But see Section 7.1.) Let (x,y,z) be the point of
common tangency with a sphere sx, sy, sz, sr. As before, the first step is to get two equations saying
(x,y,z)isontheellipsoidandsphere,andtwomorefrompartialderivatives.Trytoeliminate(x,y,z).
Asbefore,wecouldconsiderfirstplugginginax = ay = az = 0(thereducedcase).
Bridgett and Cejtin both reported complete failure after many hours, even with the reduced case.
Fermat/Dixonsolvesthereducedcaserathereasily,130susing113MBofRAM.Asachallenge,Lewis
had Fermat/Dixon solve the full case. It took about 400MB of RAM and 438min. Determinant cd had
1.8millionterms,andthefinalanswer80,372terms.
Recallthatthisisalljuststepone.
Next, as with the ellipse case in two dimensions, we made up an example and ran a multivariate
Newton’smethod.Thisfinishedrathereasilyandyields,forexample,Fig.6.
Problem5. Fourlinesgiven,findaspheretangenttoall.
A line b is defined by a point on the line (bx,by,bz) and a vector (through the origin) parallel to the
line(bxn,byn,bzn)(seeFig.7).
Firststep:Fourequationsthatapoint(x,y,z)mustsatisfyifitliesonthespheresx,sy,sz,r andthe
linebandisthepointoftangencyofthelinetothesphere:
(x −sx)2+(y −sy)2+(z−sz)2−r2 = 0
bxn(x −sx)+byn(y −sy)+bzn(z−sz) = 0
110 R.H.Lewis,S.Bridgett/MathematicsandComputersinSimulation61(2003)101–114
Fig.6.Solutionspheretangenttofourgivenellipsoids.
Fig.7.Parameterizationoflinetangenttosphereinthree-dimensionalspace.