After administration and distribution, a drug has therapeutic effect if the molecules have affinity and selectivity for its pharmacological target. But in the same time it is very important to realise an optimal concentration interval because concentrations above or below this interval can produce toxic manifestations or no therapeutic effects. The way by which a drug is released from a particular formulation can have a remarkable effect on its efficacy or toxic effects. In this context, the aim is to reduce to minimum the risks of administration by maintaining drug level within a desired range. Following conventional drug administration is possible to avoid high variation of plasma concentration if it used multidose therapy. The interval between doses is calculated using a pharmacokinetic parameter named half-life time of the drug. In the last years, the interdisciplinary research that combines chemistry, pharmacology, and molecular biology released new pharmaceuticals forms which provide a specific quantity of a therapeutic substance for a prolonged period of time to a target area within the body.

One of these forms are the drug delivery systems (DDS), based on biocompatible polymers. Depending on DDS and application type, studies reveled that several phenomena occur simultaneously or concurrently. These phenomena, mentioned in the order of their appearance, are as follows:

water diffusion inside DDS due to water concentration gradient between release environment and DDS;

In most of the theories developed so far, diffusion is considered the dominant phenomena, in the approximation that the effect of all other phenomena is negligible.

If the diffusion takes place only on x-axis, Fick's first law can be written in the form(1)J=−D∂C∂x,where J is the diffusion flux, C the drug concentration, ∂C/∂x the drug concentration gradient, and D the diffusion coefficient. Assuming that the ensemble, DDS and release environment, is isolated, the concentration variation in time ∂C/∂t is numerically equal with the flux difference ∂J/∂x:(2)∂C∂t=−∂J∂x.

By combining these equations, the diffusion equation, that is, Fick's second law, results:(3)∂C∂t=∂∂xD∂C∂x.

The generalized diffusion equation, for the case of diffusion along all three axes x, y, and z can be written as [2](4)∂C∂t=∂∂xD∂C∂x+∂∂yD∂C∂y+∂∂zD∂C∂z.

Since water and drug diffusion take place simultaneously, it can be further generalized to(5)∂Cm∂t=∂∂xDm∂Cm∂x+∂∂yDm∂Cm∂y+∂∂zDm∂Cm∂z,where C_m is the concentration and D_m is the diffusion coefficient for water (m = 1) and drug (m = 2).

In the particular case of a cylindrical DDS, the diffusion equation can be written under the form(6)∂Cm∂t=1r∂∂rrDm∂Cm∂r+∂∂θDmr∂Cm∂θ+∂∂zrDm∂Cm∂z,where C_m and D_m have the same significance as above, r is the radial coordinate, z is the axial coordinate, θ is the angular coordinate, and t is time [3].

Crank [4] solved these equations, assuming constant diffusion coefficient and known initial and boundary conditions, offering an extensive number analytical solution, for different geometries [4].

But, in reality, and, implicitly in the case of DDS, the diffusion coefficient is influenced by time, position, and solute concentration, in which case the diffusion equations can not be solved, due to the high number of variable dependence.

To overcome this aspect, having in view the need to predict in an easier manner the drug release, empirical and semiempirical equations were used: the Higuchi [5] is as follows:(7)MtM∞=KH·t1/2and the Korsmeyer-Peppas [6] is as follows:(8)MtM∞=KP·tn,where M(t) is the released drug amount at t, M_∞ is the released drug amount after very large time intervals, tending to infinity, equal, in most cases, with the drug amount loaded into polymeric matrices, K_H is the Higuchi diffusion constant, K_P is the Korsmeyer-Peppas constant, an indicator of release rate, and n is an exponent that depends on the shape of the polymeric matrix and can indicate the drug release mechanism. The Korsmeyer-Peppas is actually a generalization of the Higuchi equation; Higuchi equation is considered applicable in the first stage of the release (up to 20%), while Korsmeyer-Peppas equation is considered applicable in the first stage of the release up to 60%, until the equilibrium plateau is reached.

These models were confirmed by a plethora of experimental data, for all types of polymeric matrices and drugs, proving their validity and confirming the existence of

burst effect phase (I) characterized by a high drug release rate (large drug amount released in a very short time), determined by high concentration gradients;

In addition to these, for long time intervals, a fourth phase has been identified, namely, degradation phase (IV), in which polymer fragments derived from the polymer matrix bonds to released drug molecules, and, consequently, a decrease of the released drug amount is observed (Figure 1) [7, 8].

The main drawback of the above equations is that they have been demonstrated from Fick's second law, considering diffusion dominant, ignoring all other phenomena, and, moreover, assuming the release medium homogeneous and isotropic in relation to diffusion.

The use of these approximations, necessary to reduce the number of variables from the equations system that characterizes the system evolution, in order to determine its solutions, is not justified in reality. Actually, it indicates the “failure” of the classical models [1–6] and therefore the “incapacity” of some mathematical procedures based on the assumption of dynamic variables differentiability, in the analysis of complex phenomena involved in drug delivery.

However, this situation can be overcome applying nonstandard mathematical procedures based either on the nondifferentiability of the dynamic variables (remaining still tributary to differential methods) through fractal type theories or on invariance groups (any “unspecified potentiality” has the concrete expression in the existence of an invariance group).

In the present paper, we will explain the last procedure, namely, being given an “actual situation” in the form of semiempirical Higuchi law for the one-axial case; we will express “not actual” situations (potentialities) in the form of a Higuchi type law for the three-axial case, as concrete expression of Lie's group existence.

The Higuchi equation is the hull of a parabolas family of the form(9)M2t=a¯2t+b¯,where a¯ and b¯ are constants dependent both on the external constrains (temperature, pressure, etc.) and on the structure of the polymeric matrices. In the selection of expression (9) the fact that this expression is universal was taken into account, meaning that it is valid for any polymeric matrices, no matter the shape, structure, and so on, and, also, it is valid for any “constrain” type rate, from the “burst type effect” domain to equilibrium plateau. Moreover, it involves an intrinsic isotropy, which does not correspond to reality, the drug release mechanism having a high degree of anisotropy. Still, one can establish a correlation between the one-axial case, associated with the isotropy hypothesis and the three-axial case, associated with anisotropy, as we will show next.

If we center around the parabolas family from (9), then it is clear that it must make its mark on a possible experimental plane geometry (M, t). This geometry can be founded on a parametric group which must make the form from relation (9) invariant. This group can be best revealed if the homogenous coordinates (M, t) are used in the form [9](10)x1x=x2y=x31,where(11)x=t+b¯y=Ma¯case in which (9) becomes(12)x22−x1x3=0.

In this situation, the conic from (12) accepts the canonic parameterization [9]:(13)x1τ2=x2τ=x31,where τ is a real parameter and its invariance group is the three-parameter group generated by the homographic transformation of the τ parameter. If this transformation is written under a more convenient form(14)τ=τ+α11−α2−α3,which highlights the unit transformation for α₁ = α₂ = α₃ = 0, then using (13) the following transformation relations for the parameters x₁ and x₂ result:(15)x1=x1+2α1x2+α12α32x1−2α31−α2x2+1−α22,x2=−α3x1+1−α2−α1α3x2+α11−α2α32x1−2α31−α2x2+1−α22,from which a continuous two variables with three-parameter group can be observed. The Lie algebra [10] is given by the operators(16)L1=2y∂∂x+∂∂y,L2=2x∂∂x+y∂∂y,L3=2xy∂∂x+2y2−x∂∂y,with the commutation relations:(17)L1,L2=L1,L2,L3=L3,L3,L1=−2L2,where inhomogeneous coordinates were taken into account in order to simplify the writing.

As it should be, the conics in (12) appear in this situation as (16) group's invariant varieties with two parameters, and this is why they are invariant only with regard to the first two operators from (16). The issue at hand is not to find the two-parameter invariant varieties families, but to find the three-axial that holds three parameters: the main masses, that is, the eigenvalues of the mass tensor. Now the masses evolution group remains to be solved, which must be isomorphic to the group from (16). In order to highlight it we must note that the main masses are the solution to the secular equation of the respective matrix, which can be written as [9](18)M3+3a1M2+3a2M+a3=0,where 3a₁, 3a₂, and a₃ are the orthogonal invariants of the masses matrix. If the masses state varies from M₁, M₂, and M₃ to M₁′, M₂′, and M₃′ then an algebra theorem [9] shows that, between the secular equations, which have the respective values as roots, a linear relation takes place, generated by the homographic transformation(19)M′=aM+bcM+dwhich gives a three-parameter group but in three variables. By writing the roots of the curve from relation (18) in the Barbilian form [11–13],(20)M′=h+εih¯k1+εik,where ε_i³ = 1, h, h- are quantities conjugated one to the other, k is a one-module complex factor, the transformation from (19) induces the quantities h, h-, and k is the real transformations(21)h′=ah+bch+d,h−′=ah−+bch−+d,k′=ch−+dch+dk,which form a three variables with three-parameter group, that is, the Barbilian group [14].

This group is simple transitive, with the infinitesimal generators given by the operators [14](22)A1=∂∂h+∂∂h−A2=h∂∂h+h−∂∂h−A3=h2∂∂h+h−2∂∂h−+h−h−k∂∂kwhich reveals for the associated Lie algebra a structure that is identical with the one from (17). Therefore the two groups are isomorphic, and operators (16) and (22) are generated by the one and the same algebra (12). Moreover, group (22), being simple transitive, is definitely measurable, its elementary measure being given by the differential three-form [14]:(23)dh∧dh−∧dkh−h−2k.

Using the elementary probability(24)dP=dh∧dh¯∧dkh−h¯2kprobabilities theory can be a priori constructed in the space of field variables h,h-,k.

This function quadratic root can be assimilated to the wave function analogue [15] and it will satisfy an equation of Schrödinger type that defines the fractal space-time geodesics [16].

The issue now at hand is to find the invariant varieties families of group (16) with three parameters, having group (22) associated as a parameters group. In our opinion these functions can provide an answer to the problem of the correlation between the one-axial behavior of the mass-time curve and the mass induced in the complex system by the one-axial experimental “strain.”

These varieties families (joint invariant functions) will be solutions of the Stoka equations [17, 18]:(25)2y∂f∂x+∂f∂y+∂f∂h+∂f∂h−=02x∂f∂x+y∂f∂y+h∂f∂h+h−∂f∂h−=02xy∂f∂x+2y2−x∂f∂y+h2∂f∂h+h−2∂f∂h−+h−h−k∂f∂k=0.

This system admits solutions of the form(26)fα0,k02=const. ,where (27)α0=x−y2h−h¯x−h−h¯y+hh¯k02=k2x−2yh¯+h¯2x−2yh+h2.

It can be observed that the last of these integrals is a one-module complex one. In principle, f can be any function which is continuous and derivable in its variables. It is not yet known what kind of interpretation can a general solution such as (26) have, but some specific integrals values from relation (27) can still be interpreted. Thus, if the experimental one-axial constrain is monotonous, then (27) must fulfill the condition y² = x, fact which leads to the specific value x = 0. In this case, the second relation (27) gives(28)k0=ky−h−y−h,from which we can write y as(29)y=h−k−hk0k−k0.

The result we obtained in this case is important mainly because it shows that y can be identified in a specific case with one of the main masses values. Indeed, if k₀ ≡ (−1, −ε, −ε²) then the situation from (20) is again reached. Therefore, we can state that in both these specific cases the mass in the one-axial constrain can be considered as one of the internal masses eigenvalues. However we can draw more from (29). If this equation is written for k₀ = −1,(30)y=h+h−k1+k,and h, h-, and k are explicitly written with regard to the main mass and also the system of (20) is solved with regard to h, h-, and k, then the following relations can be found:(31)h=−M2M3+εM3M1+ε2M1M2M1+εM2+ε2M3,k=M1+ε2M2+εM3M1+εM2+ε2M3.

These can be related to the above-mentioned parameters by means of relations:(32)k=e−3iξh=13M1+M2+M3+16M1−M22+M2−M32+M1−M321/2sin⁡3ξ−icos⁡3ξ,where(33)tan⁡ξ=2M1−M2−M33M2−M3is the “anisotropy” angle.

If we use (32) in (30), we obtain the following:(34)y=13M1+M2+M3+12M1−M22+M2−M32+M1−M321/2sin⁡ξ.

The first term of (34) corresponds to the average drug mass released on the main directions, while the second term corresponds to the average drug mass released on the “mixed” directions. In other words, the first term reflects the linear evolution of the release system, that is, the dominance of individual effects and isotropy, in the first moments of the release (burst effect phase), while the second term reflects the nonlinear evolution, that is, dominance of collective effects and anisotropy, at higher time moments (swelling and equilibrium phases).

The dominance of the linear effects, which implies the functionality of the relation(35)y⟶M1=a¯t+b¯1/2results from (34) through the annulment of release system anisotropy, that is, imposing the restriction(36)sin⁡ξ⟶0,2M1⟶M2+M3.

On the contrary, “reconsidering” the nonlinear behavior (but still remaining tributary to the linear one) through a proper choice of normalization parameters, the functionality of a relation of the following type can be induced:(37)y=Aτ1/2snBτ1/2;s,where sn is the Jacobi elliptic function of module s [19], τ is the normalized time, and A and B are constants specific for the release system. The module s of the elliptic function sn is a measure of system nonlinearity and, implicitly, a measure of the release degree; thus, one can also explain the release mechanism. Two extreme situations can be made explicit by the elliptic function degenerations as functions of its module values:(38)y=A0τ1/2sin⁡B0τ1/2,for  s⟶0,y=A1τ1/2tanh⁡B1τ1/2,for  s⟶1.

We present the graphical representations of relation (37), in three-dimensional (Figure 2(a)) and contour plot (Figure 2(b)) formats. Table 1 reveals the curves of equally concentration for the released drug.

Accepting the Higuchi semiempirical law for the one dimensional case, it is generalized to the three-axial case based on concept of joint invariant functions at the simultaneous actions of two isomorphic groups SL(2R). In this context, a new release law, independent of the polymer matrix and/or drug types, is deducted and validated by the experimental data from the literature. Thereby, also the behavior of biological structures is explained, with a high nonlinear character, under the action of drug delivery systems [20–27].