To prove the statement that {a,b,c} is a solution triple iff
(abc)(abc+a+b+c)+(ab+ac+bc) is a square, Philip Gibbs wrote:

       (abc)(abc+a+b+c) + (ab+ac+bc) = (ab+1)(ac+1)(bc+1) - 1

so the "only if" part is easy. The "if" part I know is true because
you can now calculate an integer d from
        
        d = 2abc + a + b + c + 2sqrt((ab+1)(ac+1)(bc+1))

The 4 integers are then a solution to the 4-tuple equation and I
claim that the product of any two of them must be one less
than a square.  To prove this, consider non-negative integers 
{a,b,c} such that (ab+1)(ac+1)(bc+1) is a square, and define

       d = 2abc + a + b + c + 2sqrt((ab+1)(ac+1)(bc+1))

Rearrange,

       ( 2abc + a + b + c - d )^2 = 4(ab+1)(ac+1)(bc+1)         

Expand this and cancel terms to get the "4-tuple equation"

     a^2 + b^2 + c^2 + d^2 - 2ab - 2bc - 2ac - 4abcd - 4 = 0

So we have a non-negative integer solution to the 4-tuple equation. 
The equation can also be written,

           ( a + b - c - d )^2 = 4(ab+1)(cd+1)                        

If we knew that (ab+1) was a square then we could deduce that (cd+1)
was also a square and similarly for (ad+1) and (bd+1) so this is
a solution to the original puzzle.

But now we only know that (ab+1)(ac+1)(bc+1) is a square. It is
necessary to prove that any solution to the 4-tuple equation in
non-negative integers has the property that the product of any two is 
one less than a square.

If we have a solution {a,b,c,d}, 0 <= a <= b <= c <= d, of the 4-tuple
equation in non-negative integers we can generate 4 other solutions in 
integers by substituting any of the four variables by the other root of 
the quadratic equation in that variable. E.g. we can replace d, which 
we will take to be the largest of the four integers, by

              d' = 4abc + 2a + 2b + 2c - d

This is a different solution in integers since 

                  (d'-d)^2 = 16(ab+1)(bc+1)(ac+1) 

which is not zero. But d' could be negative, then the 4-tuple {a,b,c,d'}
would be a solution for which a,b,c are non-negative and d' is negative
but from

         ( a + b - c - d' )^2 = 4(ab+1)(cd'+1) 
        
The right is non-negative => cd >= -1, so either c = 0 or cd' = -1
If cd' = -1 then d'=-1, c=1 then a = b = 0.  If d' != -1 then c=0, 
similarly a=b=0 then d' = -2.

So there are exactly two solutions in integers for which a,b,c are
non-negative and d is negative, namely {0,0,0,-2} and {0,0,1,-1}

Next we can show that when we replace d by d', d' is less than c

         dd'   =   (c - a - b)^2 - 4(ab+1)    <    c^2    

( given that c >= b >= a >= 0)  so d >= c => d' < c.

Now we have done enough to show that given a solution { a,b,c,d } we 
can repeatedly substitude the greatest of the four numbers by its 
other root giving a smaller solution each time until we get down to 
one of the two solutions where one of the four numbers is negative, 
i.e.
                {0,0,0,-2} or {0,0,1,-1}

These two solutions certainly have the property that the product
of any two is one less than a square.

To complete the proof we need to show that { a,b,c,d } has the
property that the product of any two is one less than a square
if the same is true for { a,b,c,d' }. This is easy to get from

         ( a + b - c - d )^2 = 4(ab+1)(cd+1)
        
and the similar equations with a,b,c permuted.

In summary the following results have been proved,

(1) Given any three non-negative integers {a,b,c} such that
(ab+1)(ac+1)(bc+1) is a square we can find a non-negative integer
d such that {a,b,c,d} is a solution if the 4-tuple equation,

         a^2 + b^2 + c^2 + d^2 - 2ab - 2bc - 2ac - 4abcd - 4 = 0
        
(2) Any solution of the 4-tuple equation in non-negative integers has 
the property that the product of any two is one less than a square.

(3) From any one solution in non-negative integers we can generate a 
tree of solutions by repeatedly replacing any of the four numbers by 
the other root of the equation as a quadratic in that variable leaving 
the other three fixed. There are exactly two trees of solutions, the 
odd tree which included the solution {0,0,1,3} and the even tree
which includes {0,0,0,2}.

The key equation for the 4-tuplets is
 
                   (a+b-c-d)^2 = 4(ab+1)(cd+1)
 
If the equation is expanded it is seen to be symmetric in a,b,c and d.
When you complete the square on this equation and factorise the
discriminant you get
 
              (2abc+a+b+c-d)^2 = 4(ab+1)(ac+1)(bc+1)
 
So Peter Montgomery's solution for d gives solutions to this equation.
 
The equation can be used to generate a large tree of solutions since for
any three of the four numbers fixed there are two possible values for d
which give solutions to the equation. There are in fact two independent
trees of solutions in integers which can be generated by successively
swapping roots starting from the solutions 1,3,0,0 or 2,0,0,0.
 
An obvious conjecture is that all 4-tuples of positive integers with 
the property that the product of any two is one less than a square are
solutions of this equation. If this could be proved then a corollary
would be that there are no 5-tuples of positvie integers with this 
property.  This would be a non-trivial result since it would solve an 
infinite class of simultaneous Pelian equations of the type Richard 
Pinch mentioned.